As it's rather quiet on the forum, here is a fun program to calculate the Sunrise and Sunset times for any date and location on the Earth's surface, using standard astronomical equations with an 08M2 (and most other PICaxes). It needs less than half of the 08M2's available resources, i.e. Program space, EPROM (table) and RAM (variables), so a more sophisticated User Interface, or code to obtain the date from a RTC or even parse data from a GPS, etc., could be added as required.

I originally devised the Program some 4 years ago from a reference NOAA spreadsheet, but have recently revised it to work at higher resolution for at any Longitude and Latitude. Porting the equations did need some creativity, since PICaxe M2s offer only 16-bit integer maths and no trignometric functions, whilst these calculations typically employ double-precision floating-point maths with lots of trigonometry and other functions.

The calculation includes the "Equation of Time", because Solar noon, and hence sunrise and sunset times, drift through a range of over half an hour compared with GMT, during the course of a year. The program also takes into account normal Leap Years (but not the century exceptions) and the overall site Longitude (not just relative to local time) since these can each add more than one minute of error (variation) to the results. Angles (Longitude, etc.) are input to a hundredth of a degree, because a tenth can represent a significant proportion of a minute in the final results. Times are displayed to a tenth of minute resolution (not accuracy!) to allow the smaller "errors" to be analysed without the uncertainty introduced by rounding to the nearest minute.

M2s do not provide any trignometric functions, but this can be considered an advantage because the required subroutines can be optimised specifically for the application. For Sine and Cosine, the use of "fractional binary" notation is particularly advantageous, where the maximum numerical range of +/- 1.0 (i.e. unity with a sign) and one revolution (i.e. 360 degrees) are represented by the maximum 16-bit binary value (i.e. 65536 in unsigned integer form). This has the great benefit that if adding two angles exceeds a full rotation (i.e. > 360 degrees) then the binary overflow still produces the correct result, with no need for any additional checking.

Sine and Cosine (and their inverses) are calculated by interpolating through an optimised lookup table of 17 points (i.e. a 16 "piecewise linear" curve) covering one quadrant (i.e. 90 degrees). These subroutines, and others for signed and fractional multiplication and division, etc. are not described in detail here, but hopefully may be added to the code snippetts section in due course.

The complete program is too large to show in-line in a single forum post, so the subroutines are listed here, with the main part held over to the next post. This is quite convenient because the subroutines can be useful, independently of the main program. They all use the variables b1 - b7 (at most), whilst the solar computation is allocated the (unused in M2s) system variables s_w1 to s_w6 (or above w13 for X2s). Therefore all the remaining variables, b8 to b27, are available for additional programming.

The subroutines use about 400 bytes of program memory, the solar calculations about the same and the largely optional "user interface" code is again similar. It should be fairly easy to customise the required User Input data from the program comments, but the following may assist:

There are some significant locations which include a half-hour TimeZone offset so the program increments in quarter-hours, but the two predefined constants, TZE and TZW, permit simple conversion of integer hours and negative (western) offsets. Daylight Saving Time is not directly supported but easily added to the TimeZone, always in a positive (i.e. easterly) direction. Longitude may be entered from 0 to 360 degrees or with a negative value for the Western hemisphere; the demonstration values use Longitude Zero (Greenwich) to highlight the Solar Noon variation on the specified day.

The threshold sun elevation for sunrise and sunset is around 0.85 degree, to correspond with the top edge of the sun (not its centre) and the significant atmospheric refraction. The angle should be set negative for the northern hemisphere but positive for the southern, or appropriate angles used for Civil, Nautical and Astronomical twilights.

Code:

```
; ** Insert main program here **
; SUBROUTINES:
;============
readdeg: ; Convert degrees * 100 (signed) to a fractional binary revolution
read b1, WORD w1 ; Angle range +/- 180 degrees
wa = 29826 ; 65536 * 65536 / 36000 / 4 (max pos value is 32767)
call smult ; Signed multiplication 1 * 1 = 1, (i.e. 32768 * 32768 = 32768)
w1 = w1 + w1 ; Double the value because 65536 represents one revolution
return
smult: ; Signed (fractional word) multiplication: w1 = w1 ** wa
if wa > $7fff then
wa = -wa
if w1 < $8000 then negmult
w1 = -w1 ; Then fall into posmult
posmult:
w1 = w1 ** wa * 2 ; 32768 represents +/-1 so multiply result by 2
return
else
if w1 < $8000 then posmult
w1 = -w1 ; Then fall into negmult
endif
negmult:
call posmult
w1 = -w1
return
sdiv: ; Signed (fractional binary) division, w1/wa := fractional result
b1 = 16 ; Clear sign flag (also initialises double-word division, if used)
if w1 > $7fff then ; It's a negative value
w1 = -w1
dsign = 1 ; Sign of result
endif
if wa > $7fff then
wa = -wa
toggle dsign
endif
if w1 >= wa then ; Modulus > 1 so there is no inverse trig function
w1 = $8000 ; Marker for overflow
else
wa = wa / 128 ; Scale divisor down to ~one byte (or use double-word div)
w1 = w1 + w1 / wa ; Maximise numerator and divide
w1 = w1 * 128 ; Scale back to a fractional word
if dsign = 1 then ; Negative result
w1 = - w1
endif
endif
return
#ifdef hiresd ; High resolution division, not currently used
ddiv: ; Basic double-word division routine WX:W1/WY
for b1 = b1 to 255 step 16
wx = wx + wx + bit31 ; Shift Numerator left (top bit lost,carry from w1)
w1 = w1 + w1 ; Shift Num/Result left (low bit = 0)
if wx => wy then ; Can subtract divisor
wx = wx - wy
inc w1 ; Set lowest bit in result register
endif
next b1 ; Exit with Result=W1,Remainder=WX,Divisor=WY,
return
#endif
symbol WIDTH = $4000 / 16 ; Segment length (quadrant angle / 16 segments)
symbol SINTAB = 30 ; Start address of Sine table
data SINTAB,(0,0,141,12,252,24,45,37,2,49,94,60,39,71,66,81,144,90) ; "Symmetrical" errors
data (1,99,125,106,243,112,83,118,144,122,157,125,117,127,255,127) ; Sine table ($7FFF = +1)
cosine: ; Enter with angle in w1 (2^16 = 360 degrees)
w1 = w1 + $4000 ; Now fall through into sine
sine: ; Enter with angle in w1 (2^16 = 360 degrees)
if w1 => $8000 then ; It's negative
w1 = -w1 ; Invert if 180 - 360 degrees
call sinpos
w1 = -w1 ; Restore negative quadrant
else ; Negative quadrant jumps on to RETURN
sinpos:
if w1 > $3fff then ; Reverse 2nd and 4th quadrants
w1 = $8000 - w1
endif
interpxy: ; Piecewise linear interpolation from X-axis to Y-axis
b1 = w1 / WIDTH * 2 + SINTAB ; Segment (cell) number (*2 bytes) + table offset
read b1,word wx,word wy ; Read LHS & RHS to give baseline & top of cell
b1 = w1 / WIDTH ; Horizontal start of segment
wy = wy - wx * 8 ; Height of cell (scaled)
w1 = -b1 * WIDTH + w1 * 8 ; Horizontal target within cell (scaled)
w1 = w1 ** wy ; Height to target within cell
w1 = w1 + wx max $7fff ; Total Positive sine result
endif
return
arccos: ; Input +1 to -1 in w1, output 0 to +half rev (180 degs)
call arcsin
w1 = $4000 - w1
return
arcsin: ; Input -1 to +1 in w1, output -quarter to +quarter rev (90 degs)
if w1 > $7fff then ; It's a negative value
w1 = -w1
call arcpos
w1 = -w1 ; Then jump to RETURN
else
arcpos:
b1 = SINTAB ; Pointer to first element (cell) in sine lookup table
interpyx: ; Piecewise linear interpolation from Y-axis to X-axis
do ; Piecewise linear (interpolated) lookup table
b1 = b1 + 2
read b1,WORD wy ; RH Side of rectangular cell
loop while wy < w1
b1 = b1 - 2
read b1,WORD wx ; LHS
b1 = b1 - SINTAB / 2 ; Number of skipped cells
wy = wy - wx / 16 ; Cell height (max 32k/10) scaled
w1 = w1 - wx * 16 ; Target height
w1 = w1 / wy max 255 * 4 ; Target width within cell, scaled *4
w1 = b1 * WIDTH + w1 ; Total Horizontal angle (WIDTH is 1024)
endif
return
showtime: ; Show hh:mm.m ; Enter with angle (65536 = 24 hours) in w1
w1 = w1 ** 14400 + 2 ; Time of day in tenths of a minute (rounded)
yy = w1 / 600 ; Hours
mm = w1 // 600 / 10 ; Minutes
dd = w1 // 10 ; Tenths of a minute
sertxd(" ",#yy,":")
sertxd(#mm,".",#dd," ")
return
#ifdef report
showswab: ; Display signed word from EEPROM with two decimal places (& advance bptr)
read bptr, WORD w1 ; Read w1 @bptr * Bug if w1 > 255 used in PE5 *
showsw1: ; Show signed w1
if w1 > 32767 then
sertxd("-") ; Negative symbol
w1 = -w1
endif
b1 = w1 // 100 ; Hundredths
w1 = w1 / 100 ; Units
sertxd(#w1,".") ; Print Integer part
w1 = b1 / 10 ; Tenths digit (then w1 is small enough to avoid PE5 bug)
b1 = b1 // 10 ; Hundredths digit
sertxd(#w1,#b1) ; Print decimal part
bptr = bptr + 2 ; Advance to next datafield
return
#endif
```

Cheers, Alan.