Source Code ภาษา Basic จากหนังสือ Astronomy with your Personal Computer โดย Peter Duffett-Smith ใครที่สนใจการเขียนโปรแกรมและมีสตางค์หน่อยก็ช่วยอุดหนุนเขาผ่านทาง amazon.com หน่อยนะครับ

1 REM
2 REM Handling Program HMOONNF from "Astronomy with your Personal Computer" by Peter Duffett-Smith
3 REM

5 DIM FL(20), ER(20), SW(20)
10 DEF FNI (W) = SGN(W) * INT(ABS(W))

15 DEF FNL (W) = FNI(W) + FNI((SGN(W) - 1!) / 2!)
20 DEF FNM (W) = .01745329252# * W
25 DEF FND (W) = 57.29577951# * W
30 DEF FNV (W) = W - FNL(W / 360!) * 360!
35 DEF FNW (W) = W - FNL(W)

36 INPUT "^P"; X
40 PRINT : PRINT : PRINT "New and Full Moon"
45 PRINT "------------------": PRINT

50 Q$ = "Approximate date (D,M,Y)......"
55 PRINT Q$; : INPUT DY, MN, YR

60 PRINT : GOSUB 6800: X = NF - .5: DJ = NI + .5

65 IF ER(1) = 1 THEN GOTO 145
70 IF X < 0 THEN X = X + 1: DJ = DJ - 1
75 FL(2) = 0: GOSUB 1200
80 Q$ = "New Moon on date ............. "
85 PRINT Q$ + DT$
90 X = X * 24! + .008333: SW(1) = 1: GOSUB 1000
95 Q$ = "...and at time (H,M) ......... "
100 PRINT Q$ + " " + MID$(OP$, 4, 6)

105 X = FF - .5: DJ = FI + .5
110 IF X < 0 THEN X = X + 1: DJ = DJ - 1
115 FL(2) = 0: GOSUB 1200
120 Q$ = "Full Moon on date ............. "
125 PRINT Q$ + DT$
130 X = X * 24! + .008333: SW(1) = 1: GOSUB 1000
135 Q$ = "...and at time (H,M) ......... "
140 PRINT Q$ + " " + MID$(OP$, 4, 6)

145 Q$ = "Again (Y or N) ............... "
150 PRINT : GOSUB 960
155 IF E = 0 THEN STOP
160 PRINT : GOTO 50

957 REM
958 REM Input yes/no routine
959 REM

960 PRINT Q$; : INPUT A$
965 IF A$ = "Y" OR A$ = "y" THEN E = 1: RETURN
970 IF A$ = "N" OR A$ = "n" THEN E = 0: RETURN
975 PRINT "What ? ": GOTO 960

997 REM
998 REM Subroutine MINSEC
999 REM

1000 IF SW(1) = -1 THEN GOTO 1050

1005 SN = SGN(X): XP = ABS(X) + 1.39E-06: XD = INT(XP)
1010 A = (XP - XD) * 60: XM = INT(A): S$ = "+"
1015 XS = INT((A - XM) * 6000!) / 100!
1020 IF SN = -1 THEN S$ = "-"
1025 A$ = RIGHT$((" " + STR$(XD)), 3)
1030 B$ = RIGHT$((" " + STR$(XM)), 3)

1035 D$ = RIGHT$(("0" + STR$(INT((XS - INT(XS)) * 100))), 2)

1040 C$ = RIGHT$((" " + STR$(INT(XS)) + "." + D$), 6)
1045 OP$ = " " + S$ + A$ + B$ + C$: RETURN

1050 SN = 1: S$ = "+"
1055 IF XD < 0 OR XM < 0 OR XS < 0 THEN SN = -1: S$ = "-"
1060 X = ((((ABS(XS) / 60) + ABS(XM)) / 60) + ABS(XD)) * SN
1065 IF NC = 0 THEN NC = 9

1070 OP$ = LEFT$(STR$(X), NC)

1075 IF SN = 1 THEN OP$ = S$ + OP$
1080 OP$ = " " + OP$: RETURN

1097 REM
1098 REM Subroutine JULDAY
1099 REM

1100 IF FL(1) = 1 THEN RETURN

1105 IF YR = 0 THEN PRINT "*** no year zero ***": GOTO 1150
1110 M1 = MN: Y1 = YR: FL(1) = 1: B = 0: ER(1) = 0
1115 IF Y1 < 1 THEN Y1 = Y1 + 1
1120 IF MN < 3 THEN M1 = MN + 12: Y1 = Y1 - 1
1125 IF Y1 > 1582 THEN GOTO 1160
1130 IF Y1 < 1582 THEN GOTO 1165
1135 IF Y1 = 1582 AND MN < 10 THEN GOTO 1165
1140 IF Y1 = 1582 AND MN = 10 AND DY < 5 THEN GOTO 1165
1145 IF MN > 10 OR DY >= 15 THEN GOTO 1160
1150 ER(1) = 1: FL(1) = 0
1155 PRINT "*** impossible date": RETURN

1160 A = INT(Y1 / 100): B = 2 - A + INT(A / 4)
1165 C = INT(365.25 * Y1) - 694025
1170 IF Y1 < 0 THEN C = FNI((365.25 * Y1) - .75) - 694025
1175 D = INT(30.6001 * (M1 + 1)): DJ = B + C + D + DY - .5
1180 RETURN

1197 REM
1198 REM Subroutine CALDAY
1199 REM

1200 IF FL(2) = 1 THEN RETURN

1210 D = DJ + .5: I = FNL(D): FD = D - I
1215 IF FD = 1 THEN FD = 0: I = I + 1
1220 IF I <= -115860 THEN GOTO 1235

1225 A = FNL((I / 36524.25) + .99835726#) + 14
1230 I = I + 1 + A - FNL(A / 4)
1235 B = FNL((I / 365.25) + .802601)
1240 C = I - FNL((365.25 * B) + .750001) + 416
1245 G = FNL(C / 30.6001): MN = G - 1
1250 DY = C - FNL(30.6001 * G) + FD: YR = B + 1899
1255 IF G > 13.5 THEN MN = G - 13
1260 IF MN < 2.5 THEN YR = B + 1900
1265 IF YR < 1 THEN YR = YR - 1
1270 ID = INT(DY): FL(2) = 1
1275 A$ = RIGHT$((" " + STR$(ID)), 2)
1280 B$ = RIGHT$((" " + STR$(MN)), 3)
1285 C$ = RIGHT$((" " + STR$(YR)), 6)
1290 DT$ = " " + A$ + B$ + C$: RETURN

6797 REM
6798 REM Subroutine MOONNF
6799 REM

6800 D0 = DY: MO = MN: YO = YR
6805 IF YO < 0 THEN YO = YO + 1
6810 MN = 1: DY = 0: FL(1) = 0: GOSUB 1100: JO = DJ
6815 IF ER(1) = 1 THEN RETURN
6820 MN = MO: DY = D0: FL(1) = 0: GOSUB 1100
6825 IF ER(1) = 1 THEN RETURN
6830 K = FNL(((YO - 1900! + ((DJ - JO) / 365!)) * 12.3685) + .5)
6835 TN = K / 1236.85: TF = (K + .5) / 1236.85
6840 T = TN: GOSUB 6855: NI = A: NF = B: NB = F
6845 T = TF: K = K + .5: GOSUB 6855: FI = A: FF = B: FB = F
6850 RETURN

6855 T2 = T * T: E = 29.53 * K: C = 166.56 * (132.87 - .009173 * T) * T
6860 C = FNM(C): B = 5.8868E-04 * K + (.0001178 - 1.55E-07 * T) * T2
6865 B = B + .00033 * SIN(C) + .75933: A = K / 12.36886
6870 A1 = 359.2242 + 360 * FNW(A) - (.0000333 + 3.47E-06 * T) * T2
6875 A2 = 306.0253 + 360 * FNW(K / .9330851)
6880 A2 = A2 + (.0107306 + 1.236E-05 * T) * T2: A = K / .9214926
6885 F = 21.2964 + 360! * FNW(A) - (.0016528 + 2.39E-06 * T) * T2
6890 A1 = FNV(A1): A2 = FNV(A2): F = FNV(F)

6895 A1 = FNM(A1): A2 = FNM(A2): F = FNM(F)

6900 DD = (.1734 - .000393 * T) * SIN(A1) + .0021 * SIN(2 * A1)
6905 DD = DD - .4068 * SIN(A2) + .0161 * SIN(2 * A2) - .0004 * SIN(3 * A2)
6910 DD = DD + .0104 * SIN(2 * F) - .0051 * SIN(A1 + A2)
6915 DD = DD - .0074 * SIN(A1 - A2) + .0004 * SIN(2 * F + A1)
6920 DD = DD - .0004 * SIN(2 * F - A1) - .0006 * SIN(2 * F + A2) + .001 * SIN(2 * F - A2)
6925 DD = DD + .0005 * SIN(A1 + 2 * A2): E1 = INT(E): B = B + DD + (E - E1)
6930 B1 = INT(B): A = E1 + B1: B = B - B1
6935 RETURN

Source Code ภาษา Basic โดย Daniel P. Franco

100 '*************************************************************************
200 '* PHASES OF THE MOON *
300 '* *
400 '* Programmer: Daniel P. Franco *
500 '* *
600 '* VERSION 1.0.0 *
700 '* March 8, 1987 *
800 '* [73307,3471] *
900 '* *
1000 '* This program calculates the phase of the moon for a given year *
1100 '* and month. The user inputs the year, the month, and the number of *
1200 '* consecutive months data are required for. Output includes Ephemeris *
1300 '* Time of each phase beginning with the new moon. *
1400 '* *
1500 '*************************************************************************
1600 '*************************************************************************
1700 '* *
1800 '* INPUT SECTION *
1900 '* *
2000 '*************************************************************************
2100 CLS
2200 DEFDBL A-Z
2300 PRINT "Enter Year:"
2400 INPUT YEAR
2500 LEAP = YEAR MOD 4'if leap = 0 then year is a leap year
2600 PRINT "Enter Month:"
2700 INPUT MONTH
2800 PRINT "Output For How Many Months:"
2900 INPUT COUNT
3000 IF LEAP <> 0 THEN 3400 ELSE 4700
3100 '**************************************************************************
3200 '* CALCULATION FOR DECIMAL YEARS *
3300 '**************************************************************************
3400 IF MONTH = 1 THEN YD = .0424375935815675#
3500 IF MONTH = 2 THEN YD = .1232059168497121#
3600 IF MONTH = 3 THEN YD = .2039742401178567#
3700 IF MONTH = 4 THEN YD = .2874804726493282#
3800 IF MONTH = 5 THEN YD = .3709867051807998#
3900 IF MONTH = 6 THEN YD = .4544929377122713#
4000 IF MONTH = 7 THEN YD = .5379991702437429#
4100 IF MONTH = 8 THEN YD = .6228743574068778#
4200 IF MONTH = 9 THEN YD = .7063805899383494#
4300 IF MONTH = 10 THEN YD = .7898868224698209#
4400 IF MONTH = 11 THEN YD = .8733930550012924#
4500 IF MONTH = 12 THEN YD = .956899287532764#
4600 GOTO 6000
4700 IF LEAP = 0 GOTO 4800
4800 IF MONTH = 1 THEN YD = .0424375935815675#
4900 IF MONTH = 2 THEN YD = .1245748714813756#
5000 IF MONTH = 3 THEN YD = .2053431947495202#
5100 IF MONTH = 4 THEN YD = .2888494272809917#
5200 IF MONTH = 5 THEN YD = .3723556598124632#
5300 IF MONTH = 6 THEN YD = .4558618923439348#
5400 IF MONTH = 7 THEN YD = .5393681248754063#
5500 IF MONTH = 8 THEN YD = .6242433120385413#
5600 IF MONTH = 9 THEN YD = .7077495445700128#
5700 IF MONTH = 10 THEN YD = .7912557771014844#
5800 IF MONTH = 11 THEN YD = .8747620096329559#
5900 IF MONTH = 12 THEN YD = .9582682421644275#
6000 K = ((YEAR + YD) - 1900) * 12.3685
6100 K = CINT(K)
6200 COUNT = K + COUNT
6300 T = K / 1236.85
6400 T2 = T ^ 2
6500 T3 = T ^ 3
6600 PI = 3.141592653589793#
6700 R = PI / 180
6800 '**************************************************************************
6900 '* SUN MEAN ANOMALY *
7000 '**************************************************************************
7100 SMA = 359.2242 + (29.10535608# * K) - (.0000333 * T2) - (3.47E-06 * T3)
7200 IF SMA > 360 THEN SMA = SMA / 360: SMA = SMA - FIX(SMA): SMA = SMA * 360
7300 '**************************************************************************
7400 '* MOON MEAN ANOMALY *
7500 '**************************************************************************
7600 MMA = 306.0253 + (385.81691806# * K) + (.0107306 * T2) + (1.236E-05 * T3)
7700 IF MMA > 360 THEN MMA = MMA / 360: MMA = MMA - FIX(MMA): MMA = MMA * 360
7800 '**************************************************************************
7900 '* MOON'S ARGUMENT OF LATITUDE *
8000 '**************************************************************************
8100 F = 21.2964 + (390.67050646# * K) - (.0016528 * T2) - (2.39E-06 * T3)
8200 IF F > 360 THEN F = F / 360: F = F - FIX(F): F = F * 360
8300 '**************************************************************************
8400 '* MEAN PHASE OF THE MOON *
8500 '**************************************************************************
8600 JD = 2415020.75933# + (29.53058868# * K) + (.0001178 * T2) - (1.55E-07 * T3) + (.00033 * SIN((R * 166.56) + (R * 132.87) * T) - ((R * .009173 * T2)))
8700 SMA = SMA * R
8800 MMA = MMA * R
8900 F = F * R
9000 '**************************************************************************
9100 '* TRUE PHASE CORRECTIONS FOR NEW AND FULL MOON *
9200 '**************************************************************************
9300 IF K - FIX(K) = 0 OR K - FIX(K) = .5 OR K - FIX(K) = -.5 THEN 9400 ELSE 11100
9400 JD = JD + ((.1734 - .000393 * T) * SIN(SMA))
9500 JD = JD + (.0021 * SIN(2 * SMA))
9600 JD = JD - (.4068 * SIN(MMA))
9700 JD = JD + (.0161 * SIN(2 * MMA))
9800 JD = JD - (.0004 * SIN(3 * MMA))
9900 JD = JD + (.0104 * SIN(2 * F))
10000 JD = JD - (.0051 * SIN(SMA + MMA))
10100 JD = JD - (.0074 * SIN(SMA - MMA))
10200 JD = JD + (.0004 * SIN((2 * F) + SMA))
10300 JD = JD - (.0004 * SIN((2 * F) - SMA))
10400 JD = JD - (6.000001E-04 * SIN((2 * F) + MMA))
10500 JD = JD + (.001 * SIN((2 * F) - MMA))
10600 JD = JD + .0005 * SIN(SMA + (2 * MMA))
10700 GOTO 14300
10800 '*************************************************************************
10900 '* TRUE PHASE CORRECTIONS FOR FOR FIRST AND LAST QUARTER *
11000 '*************************************************************************
11100 JD = JD + (.1721 - .0004 * T) * SIN(SMA)
11200 JD = JD + .0021 * SIN(2 * SMA)
11300 JD = JD - .628 * SIN(MMA)
11400 JD = JD + .0089 * SIN(2 * MMA)
11500 JD = JD - .0004 * SIN(3 * MMA)
11600 JD = JD + .0079 * SIN(2 * F)
11700 JD = JD - .0119 * SIN(SMA + MMA)
11800 JD = JD - .0047 * SIN(SMA - MMA)
11900 JD = JD + .0003 * SIN(2 * F + SMA)
12000 JD = JD - .0004 * SIN(2 * F - SMA)
12100 JD = JD - 6.000001E-04 * SIN(2 * F + MMA)
12200 JD = JD + .0021 * SIN(2 * F - MMA)
12300 JD = JD + .0003 * SIN(SMA + 2 * MMA)
12400 JD = JD + .0004 * SIN(SMA - 2 * MMA)
12500 JD = JD - .0003 * SIN(2 * SMA - MMA)
12600 '*************************************************************************
12700 '* ADDITIONAL FIRST QUARTER CORRECTION *
12800 '*************************************************************************
12900 IF K >= 0 AND K - FIX(K) = .25 THEN 13100 ELSE 13000
13000 IF K < 0 AND K - FIX(K) = -.75 THEN 13100 ELSE 13600
13100 JD = JD + .0028 - .0004 * COS(SMA) + .0003 * COS(MMA)
13200 GOTO 14300
13300 '*************************************************************************
13400 '* ADDITIONAL LAST QUARTER CORRECTION *
13500 '*************************************************************************
13600 IF K >= 0 AND K - FIX(K) = .75 THEN 13800 ELSE 13700
13700 IF K < 0 AND K - FIX(K) = -.25 THEN 13800 ELSE 14300
13800 JD = JD - .0028 + .0004 * COS(SMA) - .0003 * COS(MMA)
13900 GOTO 14300
14000 '*************************************************************************
14100 '* CALENDAR DATE CALCULATION *
14200 '*************************************************************************
14300 JD = JD + .5
14400 Z = INT(JD)
14500 FRAC = JD - FIX(JD)
14600 IF Z < 2299161! THEN A = Z
14700 IF Z >= 2299161! THEN ALPHA = INT((Z - 1867216.25#) / 36524.25)
14800 IF Z >= 2299161! THEN A = Z + 1 + ALPHA - INT(ALPHA / 4)
14900 B = A + 1524
15000 C = INT((B - 122.1) / 365.25)
15100 D = INT(365.25 * C)
15200 E = INT((B - D) / 30.6001)
15300 DOM = B - D - INT(30.6001 * E) + FRAC
15400 IF E < 13.5 THEN M = E - 1
15500 IF E > 13.5 THEN M = E - 13
15600 IF M > 2.5 THEN Y = C - 4716
15700 IF M < 2.5 THEN Y = C - 4715
15800 DAYINT = INT(DOM)
15900 DAYFRAC = DOM - FIX(DOM)
16000 TOTSEC = DAYFRAC * 86400!
16100 TOTHOURS = (TOTSEC / 60) / 60
16200 HOUR = INT(TOTHOURS)
16300 MINLEFT = TOTHOURS - FIX(TOTHOURS)
16400 TOTMIN = (MINLEFT * 60)
16500 MIN = INT(TOTMIN)
16600 SECLEFT = TOTMIN - FIX(TOTMIN)
16700 SEC = (SECLEFT * 60)
16800 IF K >= 0 AND K - FIX(K) = 0 THEN PHASE$ = "NEW MOON"
16900 IF K >= 0 AND K - FIX(K) = .25 THEN PHASE$ = "FIRST QUARTER"
17000 IF K >= 0 AND K - FIX(K) = .5 THEN PHASE$ = "FULL MOON"
17100 IF K >= 0 AND K - FIX(K) = .75 THEN PHASE$ = "LAST QUARTER"
17200 IF K < 0 AND K - FIX(K) = 0 THEN PHASE$ = "NEW MOON"
17300 IF K < 0 AND K - FIX(K) = -.75 THEN PHASE$ = "FIRST QUARTER"
17400 IF K < 0 AND K - FIX(K) = -.5 THEN PHASE$ = "FULL MOON"
17500 IF K < 0 AND K - FIX(K) = -.25 THEN PHASE$ = "LAST QUARTER"
17600 PRINT USING "#### ## ## ## \ \ ## \ \ ##.## \ \ \ \"; Y, M, DAYINT, HOUR, "Hours", MIN, "Min.", SEC, "Sec.", PHASE$
17700 K = K + .25
17800 IF K = COUNT GOTO 18000
17900 GOTO 6300
18000 END
 

Source Code ภาษาC โดย Dylan MacHattie แปลงจาก Source Code ภาษา Basic ของ Daniel P. Franco

/****************************************************************************
* PHASES OF THE MOON in C *
* *
* Translated from the Basic program, ver 1.0.0 (March *
* 8, 1987) by Daniel P. Franco, to Turbo-C by: *
* Dylan MacHattie *
****************************************************************************/

#include
#include
#include

#define PI 3.141592653589793
#define SOLY 365.2421875 /* Solar year */

void syntax_err(void);

struct dat {
int month;
int year;
} show;

void main(int argc, char *argv[])
{
double d, k, t, t2, t3, r, sma, mma, f, jd, z, frack, frac;
double a, b, e, dom, y, m, tothours, day, totmin, hour, sec, min, junk;
char *phase[4] = {"New moon","First quarter","Full moon","Third quarter"};
int i, leap=0, last, first;

puts(" Moon Phase Calculator");

if ((argc < 3)||(argc > 4)) syntax_err();

show.month = atoi(argv[1]);
if ((show.month<1)||(show.month>12)) syntax_err();
if ((show.year = atoi(argv[2])) > 9999) syntax_err();
if (argc == 4)
{
if ((last=atoi(argv[3])) > 9) syntax_err();
}
else last = 1;

last = (last + show.month - 1)%12 + 1;
first = (show.month + 10)%12 + 1;
if (first==12) y = (double)(show.year-1);
else y = (double)show.year;

if ((!(show.year%400))||
((!(show.year%4))&&(show.year%100))) leap = 1;

/* decimal year calculation */
switch (first)
{
case 1: d = 15.5; break;
case 2: d = 45+.5*leap; break;
case 3: d = 74.5+leap; break;
case 4: d = 105+leap; break;
case 5: d = 135.5+leap; break;
case 6: d = 166+leap; break;
case 7: d = 196.5+leap; break;
case 8: d = 227.5+leap; break;
case 9: d = 258+leap; break;
case 10: d = 228.5+leap; break;
case 11: d = 319+leap; break;
case 12: d = 349.5+leap;
}
k = (y - 1900 + d/SOLY) * 12.3685;
if (modf(k,&junk)<0.5) k = floor(k);
else k = ceil(k);

do
{
t = k/1236.85;
t2 = pow(t,2);
t3 = pow(t,3);
r = PI/180;

/* sun mean anomaly */
sma = r * fmod((359.2242+(29.10535608*k)-(.0000333*t2)-(.00000347*t3)),360);

/* moon mean anomaly */
mma = r * fmod((306.0253+(385.81691806*k)+(.0107306*t2)+(.00001236*t3)),360);

/* moon's argument of latitude */
f = r * fmod((21.2964+(390.67050646*k)-(.0016528*t2)-(.00000239*t3)),360);

/* mean phase of the moon */
jd = 2415020.75933+(29.53058868*k)+(.0001178*t2)-(1.55E-07*t3)
+(.00033*sin((r*166.56)+(r*t*132.87)-(r*t2*.009173)));

/* true phase corrections for new and full moon */
frack = modf(k,&junk);
if ((frack==0)||(frack==0.5)||(frack==-.5))
{
jd=jd+((.1734-.000393*t)*sin(sma))+(.0021*sin(2*sma))-(.4068*sin(mma))
+(.0161*sin(2*mma))-(.0004*sin(3*mma))+(.0104*sin(2*f))-(.0051*sin(sma+mma))
-(.0074*sin(sma-mma))+(.0004*sin((2*f)+sma))-(.0004*sin((2*f)-sma))
-(6.000001E-04*sin((2*f)+mma))+(.001*sin((2*f)-mma))+.0005*sin(sma+(2*mma));
}

/* true phase corrections for first and third quarter */
else
{
jd=jd+(.1721-.0004*t)*sin(sma)+.0021*sin(2*sma)-.628*sin(mma)+.0089*sin(2*mma)
-.0004*sin(3*mma)+.0079*sin(2*f)-.0119*sin(sma+mma)-.0047*sin(sma-mma)
+.0003*sin(2*f+sma)-.0004*sin(2*f-sma)-6.000001E-04*sin(2*f+mma)
+.0021*sin(2*f-mma)+.0003*sin(sma+2*mma)+.0004*sin(sma-2*mma)
-.0003*sin(2*sma-mma);
}

/* additional first quarter correction */
if (((k>=0)&&(frack==0.25))||((k<0)&&(frack==-0.75)))
{
jd=jd+.0028-.0004*cos(sma)+.0003*cos(mma);
}

/* additional third quarter correction */
else
{
if (((k>=0)&&(frack==0.75))||((k<0)&&(frack==-0.25)))
{
jd=jd-.0028+.0004*cos(sma)-.0003*cos(mma);
}
}

/* calendar date calculation */
jd=jd+.5;
z=floor(jd);
frac=modf(jd,&junk);
if (z<2299161) a = z;
else
{
a = z + 1 + floor((z-1867216.25)/(100*SOLY)) - floor((z-1867216.25)/(400*SOLY));
}
b = a + 1524 - floor(365.25 * floor((a + 1401.9)/365.25));
e = floor(b/30.6001);
dom = b - floor(30.6001 * e) + frac;
if (e<13.5) m = e - 1;
else m = e - 13;
if (m>2.5) y = floor((a+1401.9)/365.25) - 4716;
else y = floor((a+1401.9)/365.25) - 4715;
tothours = 24 * modf(dom,&day);
totmin = 60 * modf(tothours,&hour);
sec = 60 * modf(totmin,&min);

i = ((int)((frack+1)*4))%4;
if (((int)m != first)&&((int)m != last))
{
printf(" %2d/%02d/%d %2d:%02d:%02d UT %s\n", (int) m, (int) day,
(int) y, (int) hour, (int) min, (int) sec, phase[i]);
}
k = k+0.25;
} while ((int)m != last);
return;
}
void syntax_err(void)
{
puts(" Command line syntax: moonph MM YYYY [N]");
puts(" where: MM is the first month to display (1=Jan, 12=Dec)");
puts(" YYYY is the starting year");
puts(" N is the number (1 to 9) of months to display (default 1)");
exit(0);
}