# Sunrise and Sunset calculation using an 08M2 :)

#### AllyCat

##### Senior Member
Hi,

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
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
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``````
The Main Program is to follow in the next post.

Cheers, Alan.

#### AllyCat

##### Senior Member
Hi,

Now the main program:

Code:
``````; Sunrise, Solar Noon and Sunset calculation for PICaxe 08M2, etc.
; Updated: March 2017 to resolve Longitude and angles to 2 decimal places, time to 0.1 minute
; June 2017 for continuous day/night, TZ in hours/4 and to write only changed EPROM data

#picaxe 08m2			; And most others
#no_data			; Comment out this line if EPROM (table) data has changed
#define showdeb			; Show debugging values (~125 bytes)
#define report			; Report Site details, etc. (~180 bytes)

; CONSTANTS :
symbol ZERO = 65536		; Use to implement negative constant declarations
symbol TZE = 4			; Eastern time zone (use 2 or 1 for fractional hours)
symbol TZW = 252		; Western time zone (For DST, REDUCE value in TZONE by 1)

; USER SITE DATA
symbol TODAY = 4			; Required Day of month (1 - 31)
symbol MONTH = 7			; Required Month (January = 1)
symbol YEAR =  17			; Year - 2000 (for leapyear correction)
symbol TZONE = TZE * 1  		; E or W Timezone + DST * 4. for datelined hemisphere (-48 to +60)
symbol SITELAT = 5144  		; Signed Latitude in hundredths of a degree, North positive (0 to 900)
symbol SITELONG = 0		; Hundredths of degrees East (-18000 to +36000)
symbol SUNELEV = ZERO - 83	; NORTH POS! Atmospheric refraction (0.6 deg) & sun radius (0.25 deg)

; USER DATA STORED IN EPROM :
symbol EDat  = 0			; Pointer to start of data in EEPROM
symbol DayP  = EDat		; Pointer to Day in EPROM  (Binary 1-31)
symbol MonP  = EDat + 1 		; Pointer to Month in EPROM (Binary 1-12)
symbol YearP = EDat + 2		; Pointer to Year in EPROM (Binary referenced to year 2000)
symbol ZoneP = EDat + 3		; Timezone in 1/4 hours + half-day (so only positive)
symbol LatP  = EDat + 4		; Latitude of site in EPROM (degrees * 100, South neg)
symbol LongP = EDat + 6		; Longitude of site in EPROM (degrees * 100,East positive)
symbol AltP  = EDat + 8		; Elevation of Sun (deg * 100),Typ.-0.83 degree for rise/set
symbol LastP = EDat + 9		; Last address in list

; SOLAR COMPUTATION VARIABLES :
symbol dd = b4			; Current day			;)
symbol mm = b5			; Current month			;} Used only for initial report
symbol yy = b6			; Current year			;) 	and "Showtime"
symbol zz = b7			; Time zone 			;)
symbol Day = s_w5		; Daynumber within 4 year cycle
symbol GS = s_w1			; Mean anomoly of Sun
symbol LS = s_w2			; Mean Lomgitude of Sun    ?LMST=Local Mean Sidereal Time
symbol Lambda = s_w3		; Eclyptic Longditude of Sun
symbol ET = s_w4			; Used for Equation Of Time calculation
symbol delta = s_w2		; Declination of Sun (needed to calculate Sun Alt)
symbol cost = s_w1		; Product of COSs
symbol sint = s_w6			; Product of SINs
symbol noon = s_w6		; Time (angle) when sun is at highest elevation
symbol tss = s_w1			; Used only for Final sunrise to noon time (angle)

; LOCAL / SUBROUTINE VARIABLES :
symbol tempb = b1			; )_Temporary/Local/parameter byte and word (used explicitly)
symbol tempw = w1		; }
symbol wx = w2			; }_Used in Trig routines
symbol wy = w3			; }
symbol wa = w2			; Auxiliary angle (word) for smult & sdiv (may overlay WX or WY)
symbol dsign = bit10		; Sign flag (in b1) compatible with double-word division routine

init:
b4 = TODAY			; Day number 1 to 31
b5 = MONTH			; Month number 1 to 12
b6 = YEAR			; Year // 100 (only // 4 used)
b7 = TZONE + 48			; Time zone in 1/4-hours, plus a day/2, DST and dateline (0 to 108)
w4 = SITELAT			; Uses ZERO- for Southern hemisphere to work around PE limitations
w5 = SITELONG 			; Uses ZERO- for Western hemisphere
w6 = SUNELEV			; Limit elevation at SunRise/Set or twilight (North positive!)
if w5 > 18000 and w5 <= 36000 then
w5 = w5 - 36000		; Map > 180 degrees back onto negative angle
endif
bptr = 4				; b4 = start of user date and location data
for b1 = DayP to LastP		; Copy user data into EPROM
if b2 <> @bptr then
write b1,@bptr	; Update the EPROM data byte when necessary
endif
inc bptr
next

lookup mm,(0,0,31,59,90,120,151,181,212,243,273,304,335),Day	; First value not used
Day = yy // 4 * 365 + Day + dd			; Day number in 4-year cycle (leap year = 0)
if Day < 61 and mm < 3 then
dec Day					; 1st March in leap year is Day 60 (ref 1 Jan = 0)
endif

#ifdef report
sertxd("Date:DMY:",#dd,"/",#mm,"/",#yy," TZone:")	; Confirm Time data
w1 = zz * 256 ** 6400 - 1200			; Hundredths of an hour, restoring 1/2 day offset
call showsw1				; Report w1, signed with 2 decimal places
sertxd(cr,lf,"Site Lat:")			; Confirm location data
bptr = EDat + 4
call showswab				; Report word @bptr, signed with 2 decimal places
sertxd(" Long:")
call showswab				; Also increments bptr for words
sertxd(" Sun Elev:")
call showswab
#endif

main:
; Convert date to part of an orbital revolution
sertxd(" Day=",#Day,cr,lf)
w1 = 0					; * PE5 simulator bugfix *
read  ZoneP,w1				; Read Time Zone (only +ve so byte sign-extends to a word)
GS = Day * 4 // 1461			; Quarter-hours in 4 years would overflow so limit to one year
GS = GS * 24 - w1 + 48			; Hours: 2^16 * years/hours = 65536/365.2422/24=7.476318
LS = GS ** 56956 + GS			; Multiply by 7.476318 / 4 = 1.86908 to give 65536 per year
;** GS = MOD(357.41 + LS,360) ** 	; Original spreadsheet formulae
GS = 65067 + LS				; For year 2020
;** LS = MOD(280.46 + LS,360) **
LS = 51074 + LS				; For 2020 = 280.46 / 360 * 2^16
w1 = GS					; LS + 2 * Sin(GS)
#ifdef showdeb
sertxd("LS=",#LS," GS=",#GS)
#endif
call sine					; Values sent and returned in w1
;** Lambda = MOD(LS + 2 * Sin(GS),360)		; Eclyptic Longditude of sun **
wa = 349					; 1.915 scaled to 180 degrees
call smult					; Signed multiplication
Lambda = w1				; sinGS coefficient
w1 = GS + GS
call sine
wa = 4 					; 0.02 scaled (unnecessary?)
call smult
Lambda = Lambda + w1
ET = -Lambda				; Used in next formula
Lambda = Lambda + LS
#ifdef showdeb
sertxd(" Lambda=",#Lambda)
#endif
;** E1 = -1.9 * Sin(GS)	 ; Already done within Lambda **
;** E2 = 2.5 * Sin(2 * lambda)  **
;** E3 = 0.053 * Sin(4 * lambda)  **
;** ET = E1 + E2 - E3  **
w1 = Lambda + Lambda			; 2.5 * Sin(2 * lambda)
call sine
wa = 449					; 2.5 / fullscale (180)
call smult
ET = ET + w1
w1 = Lambda * 4
call sine
wa = 9					; 0.053 fullscale 180
call smult
ET = ET - w1				; Equation Of Time (Minutes = -**1440)
#ifdef showdeb
sertxd(" EoT= ",#ET," ")
#endif
;** delt1 = Sin(obliq) **
;** delt2 = Delt1 * Sin(lambda) **
;** delta = Arcsin(delt2) **
;** obliq = 23.43 **
;	w1 = 4267				; Obliqiuity scaled (65536 * 23.43 / 360)
;	call sine					; Fixed constant so can pre-calculate
Delta =  13030   				; Or w1 from previous lines
w1 = Lambda
call sine
wa = Delta
call smult
call arcsin
Delta = w1
#ifdef showdeb
sertxd(" Delta=",#w1," ")
#endif
;** cost = Cos(Lat) * Cos(delta) ; Divisor **
;** sint = Sin(Lat) * Sin(delta) ; Numerator **
b1 = LatP			; Get site Latitude
call readdeg			; Start cos(Lat) * cos(delta)
call cosine
cost = w1			; Store to use in numerator
w1 = delta
call cosine
wa = cost
call smult				; -1 < denominator < +1   (+ or - sign)
cost = w1
b1 = LatP			; Reload because b1 used in cosine routine
call sine				; Now calculate sin(Lat) * sin(delta)
sint = w1
w1 = delta
call sine
wa = sint
call smult
sint = w1
;** sint = Sin(Alt) - sint **
;** cost = sint / cost **
b1 = AltP				; 90 degs -> \$4000
call readdeg			; Get Sun Elevation
call sine
w1 = w1 - sint			; Numerator
sint = w1
#ifdef showdeb
sertxd(cr,lf,"Num=",#sint," Div=",#cost)
#endif
;  Now calculate sint/cost			; No result if |sint| > |cost|
w1 = sint
wa = cost
call sdiv				; Signed division w1 / wa
if w1 = \$8000 then 		; No sunrise and/or sunset (no arcsin/cos for value > 1 )
sertxd(cr,lf)
if dsign = 1 then
sertxd("Sun All-")
else
sertxd("No Sun to")
endif
sertxd("day",cr,lf)
else
#ifdef showdeb
sertxd(" Num/Div=",#w1)
#endif
call	arccos
tss = w1				; Daylight length / 2
#ifdef showdeb
w1 = w1 + w1
sertxd(cr,lf,"Daylight=")
call showtime
#endif

showsun:
b1 = LongP			; Get Site Longitude
call readdeg			; Returned in w1   (uses wa)
read zoneP,wy	   		; Calculate TimeZone angle (24 hours = 65536) offset 1/2 day
wy = wy AND 255			; Bug fix for PE5 simulator
wa = wy / 3 			; Multiply by 65536/96 = 682.6667
wa = wy * 683 - wa
;		Noon = \$8000 - ET - w1 + wa	; Original formula with signed timezone
Noon = wa - ET - w1		; Timezone (wa) includes a half-day offset
w1 = Noon - tss			; Sunrise
sertxd(cr,lf,"SunRise=")
call showtime			; Angle in w1 to hh:mm.m
#ifdef report
w1 = noon
sertxd(cr,lf,"Noon=")
call showtime
#endif
w1 = Noon + tss			; Sunset
sertxd(cr,lf,"SunSet=")
call showtime
endif
end``````
In the above calculation, sunrise and sunset are estimated with reference to the local noon, but the astronomical data values change a little during the day. A rigorous solution is to repeat the computation for the updated times, but a simpler method (for PICaxe) is to interpolate between values calculated for the previous and subsequent midnights. For higher latitudes, this method can also detect the two days a year on which there is only a sunrise or a sunset. A version using this, and other refinements such as a 32 element Sine function, higher-precision division and any corrections that others might discover, may be added in due course.

Usually, sunrise and sunset time are only quoted to the nearest minute (which is good enough for almost all practical purposes) but for testing, the following web link might be of interest:

https://sunrise-sunset.org/search?location=Greenwich,+London,+United+Kingdom

Cheers, Alan.

#### premelec

##### Senior Member
Very good! Now how about tides? ;-0 [for Snoopy's benefit...]

#### tmfkam

##### Senior Member
Wow! Very impressive.

#### julianE

##### Senior Member
Very impressive indeed. I enjoy making scientific type instruments and much of it is dependent on solar activity. I noticed that my magnetic sensor changes with the sun.
I'm in a midst of reading a novel set in Victorian England and came across a concept I was not aware of, Telluric currents, there are a few articles on the web I'm working through, it's a fascinating phenomenon. Very much diurnal, here is a quote from Wikipedia,

"These currents are known to have diurnal characteristics wherein the general direction of flow is towards the sun. Telluric currents continuously move between the sunlit and shadowed sides of the earth, toward the equator on the side of the earth facing the sun (that is, during the day), and toward the poles on the night side of the planet."

I'm thinking of ways of making a picaxe based sensor to monitor the Telluric activity and might incorporate your fine code to track Sunrise/Sunset.

I also have to better my understanding of interfacing op-amps to picaxe which you helped me with a while back with a geophone sensor.

#### newplumber

##### Senior Member
Hi Allycat
I think you solved a mystery again ... which of course will be another project when i find time
but I plan on making a globe that will light up where the sun is shining and using your awesome code i think it can be possible
and I plan on using APA102 strips with of course 4 or less picaxe 20x2's since you great people got me started in I2c
Since I tried to find a translucent globe the only thing i can figure using is a beach ball one with extra support
but thanks alot for sharing the sunrise
I drew a picture of idea and of course my globe looks like i kicked it in front of a train but should get the idea
your credit card maxed out buddy

#### Attachments

• 20.9 KB Views: 30

#### AllyCat

##### Senior Member
Hi,

First a quick "advert" that I've now posted a description of the first linear interpolation subroutine here, as used for the Sine and Cosine functions. The "reverse" interpolation, used for arcsin/cos etc., may follow later but is not high on my list of priorities.

Thanks for the replies. I'm not sure if the request for "tides" was serious, but (at least for some parts of the world) they are very dependent on geographical data so not really predictable by a PICaxe. For example, for our English Channel there can be a High tide at one end and Low at the other, only about 300 miles (500 km) apart. That produces "interesting" tides in the middle with Southampton (not far from where Snoopy has been setting sail) having a reported "Four High tides a day".

The globe idea reminds me of a project I planned many, many years ago. At the time I was receiving "imagery" (photos of a kind) directly from the early experimental weather satellites, which were in relatively low orbits so had a rather limited field of view and radio range. FWIW, I had written my own "Look Angle" tracking data from first principles (3D Pythagoras and Trigonometry) in Basic, running on a Time Shared (dial-up) IBM 360 with Teletype output. The idea was to illuminate the inside of a globe rather like day/night, but with a much smaller "spot" having a radius corresponding to the "horizon" of the satellite (around 2000 miles for those satellites and even less for ones like the ISS now).

The globe seemed easy, because at the time one could buy an illuminated translucent world globe, about 10 inches (25 cms) in diameter for a very modest price at one of our stationery stores, so I bought one. However, I soon discovered two "issues"; the first was that the inside of the globe was printed with "political" (national) regions so when the light was on then they were all that was visible. A clever idea, but I was interested in mountains and deserts, etc., not who "owned" what. The other issue was that the two hemispheres were welded together with only a small hole at the base to accept the lamp (holder). The mains voltage lamp had a long filiament (certainly no white LEDs then) and any "projector" needed quite complex optics which certainly wouldn't fit through the hole in the base. Of course I could have cut the globe into halves, but that rather defeated the "elegance" and certainly the convenience of the concept. So that project got "shelved"

I'm not sure if I even thought how I would "steer" the light, but these days two servos, or more probably steppers, and a LED might be all that's required. With any "continuously" rotating system there may be a need for slip rings or similar to couple the electrical signals, or perhaps just use very flexible wires, switch off the light for a few seconds every 24 hours and sweep the motors back to the start. Or rotate the globe and keep the light stationary? Maybe another application for the tiny stepper motors and driver board idea I wrote about some time ago (yet another project awaiting an update). Those steppers (and others like them) have 20 steps/revolution and Rev Ed sell packs of nice 1:3 ratio gears, so two-degree steps (20 x 3 x 3 = 180 steps per revolution) might be quite easy to engineer.

But that idea really needs a calculation of Sun Elevation and Azimuth for any time/day, which indeed is yet another item on my "to do" list.

Cheers, Alan.

#### newplumber

##### Senior Member
Hi Thanks AllyCat
I do like the 1 led trick then the light would show more perfectly 180 degrees and I would keep the led stationary maybe have the globe with out south pole
I guess for me I was thinking like a 30" beach ball ...fill it up and use epoxy clear on the outside to hold the shape (might have many loud bangs to make a close to perfect shaped circle) Someday I should start a new post when I get time/parts because I would definitely need your help coding it...but so far i think the mechanical parts should work like a servo for angling the led tilt for winter summer and 1 stepper motor for rotation on some crazy slide bearing with a gear or timing belt.
Sometime I will design some unreadable drawing

#### hippy

##### Technical Support
Staff member
if w1 > \$7fff then ; It's a negative value

One thing I found with a recent project having to deal with negative values was the very useful #DEFINE ...

Code:
``````#Define Negative(wvar) wvar >= \$8000
#Define Positive(wvar) wvar <  \$8000``````
Code:
``````If Negative(w1) Then
w1 = -w1
End If``````
Adding "#Define Negate(wvar)" was also handy for keeping some operations on a single line ...

Code:
``#Define Negate(wvar) wvar ^ \$FFFF + 1``
Code:
``````w1 = Negate( -w1 * 2 )

w1 = Negate( Negate(w1) * 2 ) ; works but more memory``````

#### BESQUEUT

##### Senior Member
One thing I found with a recent project having to deal with negative values was the very useful #DEFINE ...
Waouh ! Very useful indeed.
Another thing I often use : it's not always necessary to use only positives numbers :
w3=w2+w1
works well, event if w1 is negative...
also :
w1=-4 is accepted by PE6
so
Code:
``````[color=Purple]w1[/color][color=DarkCyan]=-[/color][color=Navy]4[/color]
[color=Purple]w2[/color][color=DarkCyan]=[/color][color=Navy]10[/color][color=DarkCyan]+[/color][color=Purple]w1[/color]
[color=Blue]sertxd([/color][color=Black]#[/color][color=Purple]w2[/color][color=Blue])[/color]``````
gives 6

Even multiply work with negatives numbers :
Code:
``````[color=Navy]#no_data
#simspeed 5

#macro [/color][color=Black]SerNbr[/color][color=Blue]([/color][color=Black]a[/color][color=Blue])
if [/color][color=Black]a [/color][color=DarkCyan]>= [/color][color=Navy]\$8000 [/color][color=Blue]then
[/color][color=Black]a[/color][color=DarkCyan]=-[/color][color=Black]a : [/color][color=Blue]sertxd ([/color][color=Red]"-"[/color][color=Black],#a[/color][color=Blue])[/color][color=Black]:a[/color][color=DarkCyan]=-[/color][color=Black]a
[/color][color=Blue]else
sertxd ([/color][color=Black]#a[/color][color=Blue])
endif[/color]
[color=Navy]#endmacro[/color]

[color=Purple]w1[/color][color=DarkCyan]=[/color][color=Navy]3[/color][color=Black]: [/color][color=Purple]w2[/color][color=DarkCyan]=[/color][color=Navy]10000 [/color][color=Black]: [/color][color=Blue]gosub [/color][color=Black]mult[/color]
[color=Purple]w1[/color][color=DarkCyan]=-[/color][color=Navy]3[/color][color=Black]: [/color][color=Purple]w2[/color][color=DarkCyan]=[/color][color=Navy]10000 [/color][color=Black]: [/color][color=Blue]gosub [/color][color=Black]mult[/color]
[color=Purple]w1[/color][color=DarkCyan]=[/color][color=Navy]3[/color][color=Black]: [/color][color=Purple]w2[/color][color=DarkCyan]=-[/color][color=Navy]10000 [/color][color=Black]: [/color][color=Blue]gosub [/color][color=Black]mult[/color]
[color=Purple]w1[/color][color=DarkCyan]=-[/color][color=Navy]3[/color][color=Black]: [/color][color=Purple]w2[/color][color=DarkCyan]=-[/color][color=Navy]10000 [/color][color=Black]: [/color][color=Blue]gosub [/color][color=Black]mult[/color]
[color=Blue]end[/color]

Mult:
SerNbr(w1)
sertxd("*")
SerNbr(w2)
sertxd("=")
w3=w1*w2
SerNbr(w3)
sertxd(13,10)
return``````
3*10000=30000
-3*10000=-30000
3*-10000=-30000
-3*-10000=30000

Last edited:

#### BESQUEUT

##### Senior Member
For Divide, we need a little more :
Code:
``````[color=Navy]#no_data
#simspeed 5

#Define [/color][color=Black]Negate[/color][color=Blue]([/color][color=Black]wvar[/color][color=Blue]) [/color][color=Black]wvar [/color][color=DarkCyan]^ [/color][color=Navy]\$FFFF [/color][color=DarkCyan]+ [/color][color=Navy]1

#macro [/color][color=Black]Divide[/color][color=Blue]([/color][color=Black]a,b[/color][color=Blue])
if [/color][color=Black]a [/color][color=DarkCyan]< [/color][color=Navy]\$8000 [/color][color=Blue]then
if [/color][color=Black]b [/color][color=DarkCyan]< [/color][color=Navy]\$8000 [/color][color=Blue]then
[/color][color=Black]a[/color][color=DarkCyan]=[/color][color=Black]a[/color][color=DarkCyan]/[/color][color=Black]b
[/color][color=Blue]else
[/color][color=Black]b[/color][color=DarkCyan]=-[/color][color=Black]b
a[/color][color=DarkCyan]=[/color][color=Black]Negate[/color][color=Blue]([/color][color=Black]a[/color][color=DarkCyan]/[/color][color=Black]b[/color][color=Blue])
endif
else
if [/color][color=Black]b[/color][color=DarkCyan]<[/color][color=Navy]\$8000 [/color][color=Blue]then
[/color][color=Black]a[/color][color=DarkCyan]=[/color][color=Black]Negate[/color][color=Blue]([/color][color=DarkCyan]-[/color][color=Black]a[/color][color=DarkCyan]/[/color][color=Black]b[/color][color=Blue])
else
[/color][color=Black]b[/color][color=DarkCyan]=-[/color][color=Black]b
a[/color][color=DarkCyan]=-[/color][color=Black]a[/color][color=DarkCyan]/[/color][color=Black]b
[/color][color=Blue]endif
endif[/color]
[color=Navy]#endmacro

#macro [/color][color=Black]SerNbr[/color][color=Blue]([/color][color=Black]a[/color][color=Blue])
if [/color][color=Black]a [/color][color=DarkCyan]>= [/color][color=Navy]\$8000 [/color][color=Blue]then
[/color][color=Black]a[/color][color=DarkCyan]=-[/color][color=Black]a : [/color][color=Blue]sertxd ([/color][color=Red]"-"[/color][color=Black],#a[/color][color=Blue])[/color][color=Black]:a[/color][color=DarkCyan]=-[/color][color=Black]a
[/color][color=Blue]else
sertxd ([/color][color=Black]#a[/color][color=Blue])
endif[/color]
[color=Navy]#endmacro[/color]

[color=Purple]w1[/color][color=DarkCyan]=[/color][color=Navy]21000[/color][color=Black]: [/color][color=Purple]w2[/color][color=DarkCyan]=[/color][color=Navy]3 [/color][color=Black]: [/color][color=Blue]gosub [/color][color=Black]Div[/color]
[color=Purple]w1[/color][color=DarkCyan]=-[/color][color=Navy]21000[/color][color=Black]: [/color][color=Purple]w2[/color][color=DarkCyan]=[/color][color=Navy]3 [/color][color=Black]: [/color][color=Blue]gosub [/color][color=Black]Div[/color]
[color=Purple]w1[/color][color=DarkCyan]=[/color][color=Navy]21000[/color][color=Black]: [/color][color=Purple]w2[/color][color=DarkCyan]=-[/color][color=Navy]3 [/color][color=Black]: [/color][color=Blue]gosub [/color][color=Black]Div[/color]
[color=Purple]w1[/color][color=DarkCyan]=-[/color][color=Navy]21000[/color][color=Black]: [/color][color=Purple]w2[/color][color=DarkCyan]=-[/color][color=Navy]3 [/color][color=Black]: [/color][color=Blue]gosub [/color][color=Black]Div[/color]
[color=Blue]end[/color]

[color=Black]Div:
SerNbr[/color][color=Blue]([/color][color=Purple]w1[/color][color=Blue])
sertxd([/color][color=Red]"/"[/color][color=Blue])
[/color][color=Black]SerNbr[/color][color=Blue]([/color][color=Purple]w2[/color][color=Blue])
sertxd([/color][color=Red]"="[/color][color=Blue])
[/color][color=Black]Divide[/color][color=Blue]([/color][color=Purple]w1[/color][color=Black],[/color][color=Purple]w2[/color][color=Blue])
[/color][color=Black]SerNbr[/color][color=Blue]([/color][color=Purple]w1[/color][color=Blue])
sertxd([/color][color=Navy]13[/color][color=Black],[/color][color=Navy]10[/color][color=Blue])
return[/color]``````
21000/3=7000
-21000/3=-7000
21000/-3=-7000
-21000/-3=7000

#### AllyCat

##### Senior Member
Hi,

Thanks, but personally I prefer to keep my programs compatible with PE5 and AxePad. I do sometimes symbolise \$7FFF as MAXPOS, or similar, but it does rather seem to be replacing one line of code with two. Also, I still need to get my head around whether w1 = Negate( Negate(w1) * 2 ) might add to the clarity of the code.

Perhaps because my PC hardware is rather ancient, I can double-click PE5, load the sunrise program, edit the date (and even the location if required), run the simulator and read the sunrise and sunset times from the Terminal before PE6 would have even finished booting up.

Yes, -3 * 10000 does indeed give a "correct" answer, but -4 * 10000 doesn't. Note that the subroutines in #1 actually use "Fractional Binary" notation, so the Signed Multiplication subroutine needs no bounds (error) checking, because 1 * 1 can't exceed 1 (or to be precise, the maximum is {+/-} 0.99997 * 0.99997 = 0.99994 ).

The Division routine does need (and includes) an error check for a result greater than unity, but that's required anyway because (the subsequent) ARCSIN of a value > 1 can't exist. That's how the program catches the incidences of "endless day" or "endless night" at high Latitudes (positive or negative).

It was actually the SMULT subroutine that caused me most concern, because it uses an (implied) GOTO (if w1 < \$8000 then {goto} negmult) which we are often being told is never necessary. But I had to settle for that because I couldn't see how to nest two IF .. THEN .. ELSE .. ENDIFs to achieve the same functionality.

Cheers, Alan.

#### newplumber

##### Senior Member
Perhaps because my PC hardware is rather ancient, I can double-click PE5, load the sunrise program, edit the date (and even the location if required), run the simulator and read the sunrise and sunset times from the Terminal before PE6 would have even finished booting up.
thats the same problem I have so I just leave PE6 up running forever ..its funny tho I thought i was the only one
then when downloading the code to the picaxe ...it gives me time to order/delivered/eat a pizza joking but
I have to say PE6 does have alot cooler features and i like it alot.

#### hippy

##### Technical Support
Staff member
Though PE6 is slower to start than PE5 it starts reasonably speedily and is perfectly okay once running for me and I am using fairly old hardware at home. The big advantage is having #MACRO and #DEFINE functions and the gains from those help any disadvantages slip away.

I personally wouldn't want to do any complex development which involved signed numbers or pseudo floating point without #MACRO and #DEFINE to help me. The gains from ease of coding from those far outweighs any start-up time.

I will admit it took me a while to get used to using PE6, some time to want to use PE6 in preference to PE5, and acceptance probably does come in using it, appreciating what it does offer, finding and realising the advantages.

And, tight-fisted as I am myself; there does come a time when one may need to spend money on a better computer system.

#### AllyCat

##### Senior Member
Hi,

Since we've rather gone Off-Topic......... Well, I doubt if I'll ever follow the Apple Mac route, which leaves me with the choice of Win10 or Linux; "Between the Devil and tbe Deep Blue Sea" really seems rather appropriate,

I do have a Win10 tablet (and larger laptops), admittedly very much at the budget end of the spectrum, which does indeed load PE6 faster than my humble netbook (and even my desktop PC). Also, the latest "update" appears to have "broken" PE5. But the required screen "real estate" is another issue; with powerful reading glasses and a Bluetooth mouse it's just about usable, but certainly not a patch on PE5 with my trusty Samsung NC10 for "mobile" use.

I don't believe that I'm quite alone (yet), which is why I still prefer to post at least my "snippetts" routines in an AxePad/PE5 compatible format...... However, I do use PE6 occasionally.

Cheers, Alan.

#### newplumber

##### Senior Member
@circuit that is awesome ...wonder if it comes with a cd of free AOL dial up internet

Hi Allycat Win10?? I'm still figuring out winXP

#### Jack Burns

##### New Member
AllyCat,

I am really impressed with your work, especially when you consider the complexity of solar calculations and the limitations of PICAXE maths.

Having tested the code in the simulators (PE 6.1.0.0 and PE 5.5.6), I am getting around a 3 minute error on sunrise and around a 2 minute error on sunset. Is this due to a limitation of PICAXE maths?

I adjusted the USER SITE DATA as shown below (I used 5148 for SITELAT so it matched the sunrise-sunset.org website). All other code is the same as your first 2 posts.
Code:
``````; USER SITE DATA
symbol TODAY = 3            ; Required Day of month (1 - 31)
symbol MONTH = 8            ; Required Month (January = 1)
symbol YEAR =  21            ; Year - 2000 (for leapyear correction)
symbol TZONE = TZE * 1          ; E or W Timezone + DST * 4. for datelined hemisphere (-48 to +60)
symbol SITELAT = 5148          ; Signed Latitude in hundredths of a degree, North positive (0 to 900)
symbol SITELONG = 0        ; Hundredths of degrees East (-18000 to +36000)
symbol SUNELEV = ZERO - 83    ; NORTH POS! Atmospheric refraction (0.6 deg) & sun radius (0.25 deg)``````

Both PE 6.1.0.0 and PE 5.5.6 produced this same output:
Code:
``````Date:DMY:3/8/21 TZone:1.00
Site Lat:51.48 Long:0.00 Sun Elev:-0.83 Day=580
LS=24063 GS=38056 Lambda=23897 EoT= 65256  Delta=3160
Num=57426 Div=19490 Num/Div=51968
Daylight= 15:16.0
SunRise= 5:28.4
Noon= 13:6.3
SunSet= 20:44.2``````

and here is the information from the internet:

When I run the same code on a PICAXE 08m2 chip (Firmware v4.A), it tells me the sun doesn't rise!!
Code:
``````Date:DMY:3/8/21 TZone:1.00
Site Lat:51.48 Long:0.00 Sun Elev:-0.83 Day=580
LS=24063 GS=38056 Lambda=24063 EoT= 0  Delta=1020
Num=0 Div=0
No Sun today``````
I'm trying to figure out where it's going wrong, but have got there yet.

Regards
Jack

#### AllyCat

##### Senior Member
Hi Jack,

I originally wrote the program over 8 years ago, with a few refinements 4 years ago, so I don't recall all the details of the calculation and potential errors. However, I did notice that the vast majority of Sunrise/Sunset calculators don't even quote figures to a resolution of better than one minute. At this time of the year, the times are changing by almost two minutes a day, so next year will be about 1/2 minute different for the same date, because of the quarter revolution of the Earth (i.e. 365.24 revolutions/year). I believe that I handled the hour/day/longitude boundaries correctly, but did say in post #1 that the tenths of a minute resolution (NOT accuracy) is only intended for analysis to avoid the "uncertainty" of 2 minutes, if figures are truncated to the nearest minute.

However, it appears that the PICaxe code is predicting a "late" sunrise and "early" sunset (with noon almost correct) which suggests that the estimation of the (part of the ) sun crossing the horizon is slightly in error. The calculation uses a figure of -0.83 degrees to take in account the radius of the sun (about 0.25 degree) and optical refraction of the Earth's atmosphere (very much a "guesstimate"). I was surprised that the refraction needed to be as much as 0.6 degree, but increasing it to 0.8 degree seems to almost "fix" the calculation, at least at this time (of the year?)

Another factor is the calculation of the Trigonometric values of SINE and COSINE which uses a lookup table (with linear interpolation) of only 16 line segments (per 90 degrees). I don't know how much (maximum) error that produces (it will vary with the exact angles) but when I revised my Code Snippet interpolation algorithm (in post #2 in 2018) I paid more attention to the accuracy and increased the lookup to 32 and 64 line segments. I don't know how much difference that might make, because it would need to be tested over a range of samples, not just a "one shot" test for an arbitrary angle / date.

I can't explain the "No Sun" result when run on a "real" PICaxe, it looks as if a "division by zero error" might not have been trapped, or even zero / zero which of course has an even more "indeterminate" value. The "Lambda" looks somewhat different and the zero for the Equation of Time (EoT) is certainly incorrect for this time of the year. The first thing I would check is that the EEPROM / DATA (used for the lookup tables) has been stored correctly (the `#NO_DATA` must be commented-out with a ";" the first time the program is run). Also swap the use of the S_Wx registers for "normal" Word registers (if sufficient are available) because I do recall seeing some issues with byte-word conversions with these S_W registers.

Cheers, Alan.

Last edited:

#### Jack Burns

##### New Member
Hi Alan,

The "No Sun" result when run on a "real" PICaxe was traced to missing Sine Table data. I made the mistake in thinking that because I had only used one location then the following line didn't apply.
#no_data ; Comment out this line if EPROM (table) data has changed

However, your following suggestion fixed the problem.
The first thing I would check is that the EEPROM / DATA (used for the lookup tables) has been stored correctly (the `#NO_DATA` must be commented-out with a ";" the first time the program is run).

Regards
Jack

#### John West

##### Senior Member
Thanks AlleyCat. Excellent work. This may be part of the answer to my PV panel sun tracking needs. In order to use minimal energy and have minimal wear and tear on the positioner, I want a tracker that knows where the sun is supposed to be. I don't want the tracker wandering around chasing clouds using a brightness detector, so a tracker using this info to help calculate sun position seems optimal.

#### AllyCat

##### Senior Member
Hi,

Yes, a "reverse" calculation (i.e. sun look-angle from the time/date) is something that's been on my "to-do" list since I wrote the above Program, particularly for PV tracking, etc.. I'm sure that a "direction prediction" should be much more reliable than attempting to "maximise" the light level, from direct measurements. However, the calculation has not proved as easy as I'd hoped; I "found" several (different) sample calculations on the web, which eventually proved to be "wrong", one from a University Lecturer's teaching notes!

Of course the first requirement is to "know" the time and date, which should be quite easy, from an I2C RTC (plenty of code available for a PICaxe). Also, I've already written code for all the required trigonometric functions, and their inverses, in the Code Snippets section. But IIRC one of the "problems" I encountered was that the inverse functions don't produce a "unique" solution (for the Azimuth), for example the ARC SIN of 0.5 can be an angle of 30 degrees or 150 degrees.

As a matter of principle, I would want any code that I write, to be accurate to (much) better than 1 degree (i.e. about the diameter of the sun). But a much lower accuracy should be sufficient for PV sun-tracking, at least of conventional flat panels. The received energy falls as the COSINE of the angle, so a 5 degree error reduces it by only about 0.5%, and 10 degrees by about 1.5%. So I think a simple "trial and error" (repeatedly "guess" the sun's elevation and put that in as the "dusk" threshold) using my original code, could estimate the sun's elevation well enough. However, I believe it was the Azimuth that caused me more issues, but it shouldn't be too difficult to estimate that to within 5 or 10 degrees?

Cheers, Alan.

#### Jack Burns

##### New Member
John,
I recently came across a 2 axis GPS solar tracker which is based on a Picaxe-20X2 and calculates the position of the sun. It certainly makes interesting reading, and all the code is available to download. I haven't tried the code myself as I was only looking to calculate sunrise and sunset, and AllyCat's code gave me all that I needed.

The solar tracker can be seen at MTM Scientific.

Regards
Jack

#### papaof2

##### Senior Member
But IIRC one of the "problems" I encountered was that the inverse functions don't produce a "unique" solution (for the Azimuth), for example the ARC SIN of 0.5 can be an angle of 30 degrees or 150 degrees.
Would it be possible to add a lookup table by month for the valid range? With 3 bytes of data (month number, minimum, maximun) for the location? Possibly best done as data entered by the user for their location?

#### AllyCat

##### Senior Member
Hi,

IIRC the "double-values" occurred each day, possibly because the maths was trying to "work backwards" from the sun elevation, which of course has the same value twice each day (except at solar noon). I believe I've seen the "MTM Scientific" PICaxe code before, but its 4000 Program Bytes and 50+ registers didn't look promising, when I'm only interested in writing code that can run on an (08)M2 .

However, their original program from David Williams is quite interesting, so I'm looking to see if it can be adapted to my approach (of using Twos Complement, Fractional Binary numbers, with interpolated tables). I prefer to use "real" maths and not "tricks" such as separating off Sign flags, or even PICaxe's "Multiply small numbers by 100" method (their Trig. function accuracy is really rather poor). So the first requirement is a good Square Root algorithm, which I did look at some months ago, but haven't documented as a Code Snippet yet. Probably another interpolated Lookup Table, but remembering that an 08M2 has only 256 bytes (or 128 words) of EEPROM Data storage.

Cheers, Alan.

#### John West

##### Senior Member
John,
I recently came across a 2 axis GPS solar tracker which is based on a Picaxe-20X2 and calculates the position of the sun. It certainly makes interesting reading, and all the code is available to download. I haven't tried the code myself as I was only looking to calculate sunrise and sunset, and AllyCat's code gave me all that I needed.

The solar tracker can be seen at MTM Scientific.

Regards
Jack
Thank you, Jack. I'll dig into it.

#### lbenson

##### Senior Member
I recently came across a 2 axis GPS solar tracker which is based on a Picaxe-20X2 and calculates the position of the sun. It certainly makes interesting reading, and all the code is available to download. . . .The solar tracker can be seen at MTM Scientific.
Great project, thanks for linking. Code is worth digging into, and around 50% comments.