Calculating the approximate phase of the moon on a range of devices

Tuesday, 7th April 2026

Calculating the phase of the moon always struck me as a fun little program for a calculator, and with the current Artemis II mission it seemed like as good a time as any to look at a few example programs for this on a range of different calculators and pocket computers.

Sharp PC-1211

I copied the following program from the book 119 Practical Programs for the TRS-80 Pocket Computer and it forms the basis of most of the subsequent programs:

  100 "MOON"INPUT "DATE M?",M,"D?",D,"Y?",Y: GOSUB "DJ"
  110 M=(J+4.867)/29.53058:M=2*(M-INT M)-1:N=ABS M
  120 USING "##.##": PRINT "MOON LIT ABOUT ";N
  130 Z$="NEW": IF M>0 LET Z$="FULL"
  140 PRINT "HEADED FOR A ";Z$;" MOON.": END
  900 "DJ"J=INT (365.2422Y+30.44*(M-1)+D+1):N=M-2+12*(M<3)
  905 Z=Y-(M<3):E=INT (Z/100):Z=Z-100E
  910 W=INT (2.61N-.2)+D+Z+INT (Z/4)+INT (E/4)-2E
  915 W=W-7*INT (W/7):X=J-7*INT (J/7)
  920 J=J-X+W-7*(X<W)+1721061: RETURN

The program has two parts; the first is to convert a date from its year, month and day components into a Julian day number which it does via the subroutine DJ. It then calculates the phase of the moon by adding an offset (4.867 in this case) and dividing by 29.53058, the length of the synodic month (lunar cycle) in days. The fractional part of the result of this calculation corresponds to the phase of the moon, and by multiplying by two and subtracting one you can get a value roughly corresponding to how illuminated the moon is on a particular date and whether it's waxing (heading for a full moon) or waning (heading for a new moon).

Incidentally, the program listing printed in the book does have a bug: line 130 just checks IF M instead of IF M>0, and so the program will nearly always report a waxing moon (heading for a full moon). This is corrected in the version of the code shown above.

Sharp PC-1500

The next page of the book where I copied the PC-1211 listing from showed a graphical representation of the moon's phase which I thought would be fun to replicate, though as the PC-1211's printer doesn't support graphics it would be a very crude representation indeed. Fortunately, the PC-1500's plotter allows for graphical output, so I added an optional plotting routine to the above program as well as some other niceties like automatically populating the month and day fields with the value from the real-time clock:

   30 T=TIME :Y=2000:M=INT (T/1E4):T=T-M*1E4:D=INT (T/1E2)
   40 L=-50
   50 WAIT 0
   60 CLS : PRINT "Year (";STR$ Y;") ";: INPUT Y
   70 CLS : PRINT "Month (";STR$ M;") ";: INPUT M
   80 CLS : PRINT "Day (";STR$ D;") ";: INPUT D
   90 CLS : WAIT 
  110 J=INT (365.2422*Y+30.44*(M-1)+D+1)
  120 N=M-2+12*(M<3)
  130 Z=Y-(M<3)
  140 E=INT (Z/100)
  150 Z=Z-100*E
  160 W=INT (2.61*N-.2)+D+Z+INT (Z/4)+INT (E/4)-2*E
  170 W=W-7*INT (W/7)
  180 X=J-7*INT (J/7)
  190 J=J-X+W-7*(X<W)+1721061
  210 P=(J+4.867)/29.53058
  220 P=2*(P-INT P)-1
  230 N=ABS P
  240 Q=INT (N*100+.5)
  250 CLS : PRINT "Moon lit about";Q;"%"
  260 Z$="full": IF P<0 LET Z$="new"
  270 CLS : PRINT "Headed for a ";Z$;" moon."
  290 IF PEEK &A000<>&C0 END
  300 WAIT 0:P$="Y": PRINT "Print output (Y/N) ";: INPUT P$
  310 IF P$<>"Y" END
  320 "MPRINT" CLS : PRINT "Latitude (";STR$ L;") ";: INPUT L
  330 IF L<-90 LET L=-90
  340 IF L>90 LET L=90
  360 CLS : PRINT "Printing...": WAIT 
  380 M$=STR$ M: IF M<10 LET M$="0"+M$
  390 D$=STR$ D: IF D<10 LET D$="0"+D$
  400 TEXT : CSIZE 3: LPRINT STR$ Y;"-";M$;"-";D$
  410 GRAPH : GLCURSOR (216/2,-216/2+15): SORGN :R=108: DEGREE 
  420 C=9: FOR A=0 TO 360 STEP 6
  430 LINE -(R*SIN (A),R*COS (A)),C
  440 C=0: NEXT A
  450 V=N*2:C0=9:C1=0
  460 IF P>=0 LET V=2-V:C0=0:C1=9
  470 FOR S=1 TO 2:C=9
  480 FOR I=-R TO R STEP 8
  490 XO=I*COS (-L):YO=I*SIN (-L)
  500 W=√(R*R-I*I)
  510 XN=W*SIN (L):YN=W*COS (L)
  520 X=XO-XN+V*XN:Y=YO-YN+V*YN
  530 IF S=1 LINE -(X,Y),C:C=0
  540 IF S=2 GLCURSOR (XO-XN,YO-YN): LINE -(X,Y),C0: LINE -(XO+XN,YO+YN),C1
  550 NEXT I: NEXT S
  560 GLCURSOR (-216/2,-216/2-25): SORGN 
  710 TEXT : CSIZE 2: LPRINT "Moon is lit about"
  720 LPRINT STR$ Q;"% and headed"
  730 LPRINT "for a ";Z$;" moon."
  740 LF 3: END

The program also prompts for latitude to rotate the drawing appropriately – the moon appears to fill with light from the right to the left in the Northern hemisphere, from top to bottom at the equator and from the left to the right in the Southern hemisphere.

HP-12C

This calculator is really designed for financial applications but it is keystroke programmable and so a phase of moon calculation program would be a good learning project. Fortunately it can already calculate the number of days difference between two dates, so that would save having to write a program to calculate the Julian day number.

Dates on the HP-12C are represented as decimal values, either MM.DDYYYY or DD.MMYYYY depending on the current calculator mode. To keep things simple, rather than calculate the correct Julian day number the number of days since the first of January 2000 is used as the reference as this can be represented as 1.012 in either mode.

This change of date meant that the offset from the Julian day number (previously 4.867) could no longer be used. As the range of moon phases (0 to 1) are mapped to being from full to full (-100% to +100%) the 0 reference value also needs to be a full moon. I found a list of moon phases for the year 2000 which put a full moon at 04:41 on 21st January, so the offset was set to 20.195 – 20 days after the 1st January is the 21st, and 04:41 is (4+(41/60))/24=0.195 hours into the day.

One other very minor change from previous programs is the use of 29.53059 as this is a closer approximation than 29.53058, but in the grand scheme of things it doesn't make much difference to the accuracy. The program can be run by typing in the desired date (e.g. 7.042026) and pressing the R/S key.

Casio fx-3800P

This is another keystroke programmable calculator, though life is made a little more difficult for us compared to the HP-12C as it doesn't have a built-in function to compute the number of days between two dates. One complication with adapting date computation algorithms to scientific calculators is that they often lack functions to truncate values to integers, and even if they do have a way to round a number (e.g. by switching to a mode that only shows a certain number of fixed decimal places and using a "modify" key to convert the displayed number to the stored one) they sometimes don't allow mode switches mid-program.

Fortunately, the Casio fx-3800P will store mode changes in its programs. Even better, you can switch between different numerical bases and the value will be truncated (not rounded) without raising an error which is the behaviour we want when handling date calculations, and when in this BASE-n mode calculations are carried out using integer division and multiplication which suits the algorithm perfectly.

The program is split into two parts: PROG I converts a date (stored in constant registers K1=year, K2=month, K3=day) into the Julian day number and stores the result in the memory register M. The second part, PROG II, converts the Julian day number stored in M into the approximate phase of the moon.

PROG I is adapted from Julian Day Numbers by Bill Jefferys, and is as follows:

PROG II is somewhat simpler, and follows the same sort of logic as previous programs:

As mentioned above the programs make fairly heavy use of mode switching to truncate values to integers. Program flow control is very limited on these programmable scientific calculators, usually only permitting a jump back to the start of the program based on a certain condition – hence the slightly clumsy month/year adjustment at the start and end of PROG I.

The Julian day number calculation returns a value that is 0.5 smaller than the value returned by the PC-1211 program that was the basis for most of the other programs (e.g. for 7th April 2026 the Bill Jefferys algorithm returns the correct 2461137.5, the algorithm in the PC-1211 program returns 2461138). To compensate for this the offset used to calculate the current phase of the moon is made 1 larger; not 0.5, as through some experimentation a value of 1 produced results that more closely matched a lunar phase calculator I found elsewhere.

That said, none of the programs above line up particularly well with any other lunar phase calendar, and if you search through days to find when the new moon, full moon and quarters are based on the values closest to 0%, 50% and 100% you'll often find yourself a day off to one side or the other. A more accurate program would be useful, which brings me to the final and most sophisticated program.

Sharp PC-1251

When leafing through old issues of La revue des Sharpentiers, a French publication about all things Sharp from the 1980s, I found an interesting program for the PC-1261: Les phases de la lune.

This program can produce an accurate calendar of moon phases, showing the dates and times of the new moon, full moon and quarters on a month-by-month basis. The article goes into detail about how it works, and I was keen to try it, but unfortunately I do not own a PC-1261! The closest machine I have is the PC-1251, as that matches the 24-column display and printer. However, there are some troublesome differences; the PC-1251 only has a single-line display, and so the PC-1261's code would need to have anything referring to the second line of the display adjusted or removed. A more significant issue is variable names; the PC-1261 supports two-character variable names, whereas the PC-1251 only supports a single character for its variable names.

I typed in the PC-1261 program and made a list of the two-character variable names it used along with a list of the single-character variable names it does not use. There were too many two-character names to fit in the space left over, so I had to make some further adjustments to reduce variable usage such as reusing the same variable in different places for different purposes or reordering code to avoid needing to store a value in an intermediate variable.

    1 "A": PRINT =LPRINT : GOTO 5
    2 "Z": PRINT =PRINT : GOTO 5
    3 REM PHASES DE LA LUNE*J.HERY D APRES J.MEEUS* EDI.20/11/85
    5 CLEAR : WAIT 100: DEGREE :U=0: DIM M$(12)*9,L$(7)*2: RESTORE 
    7 FOR I=1 TO 12: READ M$(I): NEXT I
   10 FOR I=1 TO 7: READ L$(I): NEXT I
   20 PRINT "**PHASES DE LA LUNE**"
   25 INPUT "ANNEE ? ";Y: INPUT "NO MOIS OU AN ? ";S$:Z=Y
   30 G=1: IF Y<1583 LET G=0
   35 USING "#####": PRINT "AN:";Y: IF S$<>"AN"PRINT "MOIS: ";M$(VAL S$)
   40 PRINT " PH.   DATE    TU.(H.M)": PRINT ":--:----------:--------:": WAIT 
   45 K=INT ((Y-1900)*12.3685)
   50 T=(Y-1899.5)/100
   60 I=2415020+29K
   65 L=.0001178TT-.000000155TTT
   70 L=L+.75933+.53058868K
   75 L=L+.00033*SIN (166.56+132.87T-.009173TT)
   80 L=L-.000837T-.000335TT
   85 N=.08084821133K
   90 N=360*(N-INT N)+359.2242
   95 N=N-.0000333TT
  100 N=N-.00000347TTT
  105 O=.07171366128K
  110 O=360*(O-INT O)+306.0253
  115 O=O+.0107306TT
  120 O=O+.00001236TTT
  125 V=.08519585128K
  130 V=360*(V-INT V)+21.2964
  135 V=V-.0016528TT-.00000239TTT
  140 K=4*(VAL S$-1): IF S$="AN"LET K=0
  145 FOR K=K TO 53
  150 J=I+7K:F=L+.38264717K
  160 P=N+K/4*29.10535608
  165 Q=O+K/4*385.81691806
  170 W=V+K/4*390.67050646
  180 IF U=0 OR U=1 GOSUB 300
  185 IF U=.5 OR U=1.5 GOSUB 340
  190 F=F+.5/1440
  195 J=J+INT F:F=F-INT F
  197 R=J+F+1.5:R=R-7*INT (R/7)+1
  200 GOSUB 400
  205 IF Y<Z GOTO 260
  210 IF S$="AN" OR M=VAL S$ GOTO 220
  215 GOTO 255
  220 IF U=0 PRINT "":P$=" NL"
  230 IF U=.5 LET P$=" PQ"
  235 IF U=1 LET P$=" PL"
  240 IF U=1.5 LET P$=" DQ"
  245 PRINT P$;"  ";L$(R);USING "###";D;M;USING "#####.##";DMS H
  255 IF M>VAL S$ AND S$<>"AN"GOTO 270
  260 U=U+.5: IF U=2 LET U=0
  265 NEXT K
  270 PRINT "": END 
  300 F=F-.4068*SIN Q
  305 F=F+(.1734-.000393T)*SIN P
  310 F=F+.0161*SIN (2Q)-.0004*SIN (3Q)
  315 F=F+.0104*SIN (2W)+.0004*SIN (2W+P)
  320 F=F-.0074*SIN (P-Q)-.0004*SIN (2W-P)
  325 F=F-.0051*SIN (P+Q)-.0006*SIN (2W+Q)
  330 F=F+.0021*SIN (2P)+.0005*SIN (P+2Q)
  335 F=F+.0010*SIN (2W-Q): RETURN 
  340 F=F+(.1721-.0004T)*SIN P+.0021*SIN (2P)
  345 F=F-.6280*SIN Q+.0089*SIN (2Q)
  350 F=F-.0004*SIN (3Q)+.0079*SIN (2W)
  355 F=F-.0119*SIN (P+Q)-.0047*SIN (P-Q)
  360 F=F+.0003*SIN (2W+P)-.0004*SIN (2W-P)
  365 F=F-.0006*SIN (2W+Q)+.0021*SIN (2W-Q)
  370 F=F+.0003*SIN (P+2Q)+.0004*SIN (P-2Q)-.0003*SIN (2P+Q)
  380 F=F+SGN (1-U)*(.0028-.0004*COS P+.0003*COS Q)
  385 RETURN 
  400 F=F+.5
  405 IF F<1 GOTO 415
  410 F=F-1:J=J+1
  415 IF G=1 GOTO 425
  420 A=J: GOTO 435
  425 B=INT ((J/36524.25)-51.12264)
  430 A=J+1+B-INT (B/4)
  435 B=A+1524
  440 C=INT ((B/365.25)-.3343)
  445 D=INT (365.25C)
  450 E=INT ((B-D)/30.61)
  455 D=B-D-INT (30.61E)+F
  460 M=E-1:Y=C-4716
  465 IF E>13.5 LET M=M-12
  470 IF M<2.5 LET Y=Y+1
  475 H=24*(D-INT D):D=INT D
  480 RETURN 
  500 DATA "JANVIER","FEVRIER","MARS","AVRIL"
  510 DATA "MAI","JUIN","JUILIET","AOUT"
  520 DATA "SEPTEMBRE","OCTOBRE","NOVEMBRE","DECEMBRE"
  530 DATA "DI","LU","MA","ME","JE","VE","SA"

Another optimisation is related to the PC-1251's lack of support for long variable names; it supports implicit multiplication in certain situations, for example 2*A can be written as 2A. Whereas the original program used T2 and T3 to store the values of T² and T³ respectively, I could instead use TT and TTT in their place and save having to use up more previous variable names. I could also replace instances of Q+Q+Q with 3Q or W+W with 2W. After making these changes the program was a few bytes smaller and had enough spare single-character variable names free to use for the remaining two-character names; the end result is somewhat harder to read, but it does run on the PC-1251 and matches the output of the PC-1261 original.

Sharp PC-1245

A very similar computer to the Sharp PC-1251 is the PC-1245. This can drive the same 24-column printer, however its display is only 16 characters wide. A bigger problem, however, is the amount of available RAM: the PC-1245 only has 1486 bytes free for a program and dynamically-named variables, and the PC-1251 program is 2112 bytes in length. Slimming the program down to fit was a considerable challenge, but here are some of the changes that were made:

• Dynamically-allocated variables were removed; a month number is shown instead of a name, and weekday names come from indexing into a string instead of storing them in an array.
• Where possible, multiple lines of code were condensed into single lines of code separated by colons.
• Constant variables with long sequences of zeroes at the start were replaced with scientific notation where it saved space (e.g. .000000155 to 155𝐄-9).
• The user-supplied month number is stored in a numeric variable S with a value of 0 to print a year instead of a string with a value of "AN" – this saves a lot of comparisons and use of VAL to convert back to a number where required.
• Parentheses around function arguments were removed where not required.
• Conditions were simplified or removed if possible, for example G=1: IF Y<1583 LET G=0 becomes G=Y>1582.
• Decorative text and comments were condensed or removed entirely.
• Support for outputting to the screen was removed; information would be cut off due to the narrower screen, so making it printer-only felt like an acceptable loss.

The resulting code is much harder to read, but the resulting program is not too much of a compromise from the original in my opinion. It comes to exactly 1486 bytes, which means it completely fills the computer's memory.

   25 "A"CLEAR : INPUT "ANNEE?";Y: INPUT "MOIS?";S
   35 Z=Y:G=Y>1582: USING "#####": LPRINT "AN:",Y: IF S LPRINT "MOIS:",S
   40 LPRINT " PH.   DATE    TU.(H.M)": LPRINT ":--:----------:--------:"
   45 K=INT ((Y-1900)*12.3685):T=(Y-1899.5)/100:I=2415020+29K
   65 L=1178€-7TT-155€-9TTT+.75933+.53058868K
   75 L=L+33€-5*SIN (166.56+132.87T-.009173TT)-837€-6T-335€-6TT
   85 N=.08084821133K:N=360*(N-INT N)+359.2242-333€-7TT-347€-8TTT
  105 O=.07171366128K:O=360*(O-INT O)+306.0253+.0107306TT+1236€-8TTT
  125 V=.08519585128K:V=360*(V-INT V)+21.2964-.0016528TT-239€-8TTT
  145 FOR K=(4S-4)*(S>0) TO 53:J=I+7K
  150 F=L+.38264717K:P=N+K/4*29.10535608:Q=O+K/4*385.81691806:W=V+K/4*390.67050646
  180 IF U=INT U GOSUB 300
  185 IF U<>INT U GOSUB 340
  190 F=F+.5/1440:J=J+INT F:F=F-INT F:R=J+F+1.5:R=INT (R-7*INT (R/7)): GOSUB 400
  205 IF Y<Z GOTO 260
  210 IF S*(M<>S) GOTO 255
  220 IF U=0 LPRINT ""
  230 P$=MID$ ("NLPQPLDQ",4U+1,2)+"  "+MID$ ("DILUMAMEJEVESA",2R+1,2)
  245 LPRINT " ";P$;USING "###";D;M;USING "#####.##";DMS H
  255 IF S*(M>S) GOTO 270
  260 U=((2U+1) AND 3)/2: NEXT K
  270 LPRINT "": END 
  300 F=F-.4068*SIN Q+(.1734-393€-6T)*SIN P+.0161*SIN 2Q-4€-4*SIN 3Q
  315 F=F+.0104*SIN 2W+4€-4*SIN (2W+P)-.0074*SIN (P-Q)-4€-4*SIN (2W-P)
  325 F=F-.0051*SIN (P+Q)-6€-4*SIN (2W+Q)+.0021*SIN 2P+5€-4*SIN (P+2Q)
  335 F=F+.0010*SIN (2W-Q): RETURN 
  340 F=F+(.1721-4€-4T)*SIN P+.0021*SIN 2P-.6280*SIN Q+.0089*SIN 2Q
  350 F=F-4€-4*SIN 3Q+.0079*SIN 2W-.0119*SIN (P+Q)-.0047*SIN (P-Q)
  360 F=F+3€-4*SIN (2W+P)-4€-4*SIN (2W-P)-6€-4*SIN (2W+Q)+.0021*SIN (2W-Q)
  370 F=F+3€-4*SIN (P+2Q)+4€-4*SIN (P-2Q)-3€-4*SIN (2P+Q)
  380 F=F+SGN (1-U)*(.0028-4€-4*COS P+3€-4*COS Q): RETURN 
  400 F=F+.5: IF F>=1 LET F=F-1:J=J+1
  420 A=J: IF G LET B=INT ((J/36524.25)-51.12264):A=J+1+B-INT (B/4)
  435 B=A+1524:C=INT ((B/365.25)-.3343):D=INT 365.25C:E=INT ((B-D)/30.61)
  455 D=B-D-INT 30.61E+F:M=E-1:Y=C-4716
  465 IF E>13.5 LET M=M-12
  470 IF M<2.5 LET Y=Y+1
  475 H=24*(D-INT D):D=INT D: RETURN

I had originally hoped to squeeze the program onto the Sharp PC-1246, but that only has 1278 bytes of program memory so I'd need to shave a further 208 bytes from the program which I don't think I'll be able to pull off without some significant reworking. It should be reasonably easy to split the program into two, and have one program perform the initial setup and calculations and then CHAIN the second half that prints the calendar from tape, but for now I think I've got enough calculators calculating phases of the moon to keep me occupied.