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
		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
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.
; AllyCat, 2013 - 2017.  Formulae adapted from "NOAA" Excel spreadsheet
; 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
		read b1,b2
		if b2 <> @bptr then
			write b1,@bptr	; Update the EPROM data byte when necessary
		endif
		inc bptr
	next		 

	read Edat,dd,mm,yy,zz				; Read Time dats from EPROM
	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 readdeg
	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.
 

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

hopefully I'm not adding garbage to your perfect post
 

Attachments

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
your last friends friend mark
 

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:
24758
24757


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,

Many thanks for your prompt and detailed reply.

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).

Thanks again for your help.

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.
 
Top