Solar Tracker with GPS Receiver and Sun Position Calculator: PICAXE-20X2

M

mmruzek

New Member
This post is a quick summary of a Picaxe Solar Tracker Project that I have been working on for a long time. The reason you may be interested in this is because I am using a PICAXE-20X2 to both receive the GPS signal and calculate the Sun's position (Azimuth & Elevation) using the acquired data.

The GPS module is the EM-411, which outputs standard NMEA data sentences. This data needs to be parsed to extract the position, time of day and date. I found alot of good information on the WWW about how to do that, including on this forum. I wrote a detailed summary of the electrical connections and provide the short Picaxe BASIC code here:

http://www.mtmscientific.com/gps.html

Using the GPS data to calculate the Sun's position was alot harder! I found a program to do the calculation in standard BASIC and then had to modify it to work with Picaxe BASIC. The main challenge was doing all the math with unsigned 2 byte integers. Also, the trigonometric calculations of Sine, Cosine and Tangent were challenging in the sense that the accuracy was important. Calculation of the Equation of Time (EOT) was also challenging. Here is a link to the page that describes this part of the project and includes the Picaxe BASIC code:

http://www.mtmscientific.com/gps_2.html

I am working on a 3rd webpage that describes how to use a MEMS Acceleration sensor to move to the correct aiming direction with a Tilt-Tilt style tracker platform. I am hoping the forum finds this information interesting and useful.

Michael
 
premelec

premelec

Senior Member
Hi Michael... welcome to this forum... though I haven't got immediate use for the specifics of your project I found it of interest; I appreciate seeing what you have accomplished! You might want to post in the "finished project" section too. Best wishes...
 
B

binary1248

Senior Member
WOW, I did a sun position calculation for azimuth and elevation years ago in Basic on a desktop, and it was difficult enough even with FP arithmetic. This was back when I was practicing celestial navigation. Although I had a hand held calculator pre-programed for the sun position based on time and position (and the stars) I did it just to learn the complexities of doing the calculations. Also I was interested in moving a mirror to keep the sun's reflection on a particular spot, maybe for a solar heater of some sort.
Any way, I am very interested in your project. Thanks for posting.
Paul
.
Also calculates alignment of a heliostat mirror.
Click to expand...
Excellant
 
Last edited:
B

binary1248

Senior Member
Hmmm, I tried your code for the picaxe 20x2 but it gives me parsing errors on the symbol def.
'Calculate position of Sun in sky (Azimuth and Elevation)
'Picaxe Version M.T.Mruzek, Original by D. Williams
'Written for a PICAXE-20X2, MTM Scientific, Inc 2017
Symbol LAT=B55 'Latitude (0-90)
Symbol N_S=BIT0 'North=0 & South=1
Symbol LONG=B54 'Longitude (0-180)
Symbol E_W=BIT1 'East=0 & West=1
Symbol MONTH=B53 'Month (1-12)

what's the problem with defining Symbol LAT = B55 ' byte 55 ? I get errors on all the Symbol defs
With (Editor V6.0.9.2)
 
Last edited:
B

binary1248

Senior Member
If I do only these short list symbols it checks ok ???
.
Symbol LAT=B55 'Latitude (0-90)
Symbol N_S=BIT0 'North=0 & South=1
Symbol LONG=B54 'Longitude (0-180)
Symbol E_W=BIT1 'East=0 & West=1
Symbol MONTH=B53 'Month (1-12)
 
M

mmruzek

New Member
binary1248 said:
If I do only these short list symbols it checks ok ???
.
Symbol LAT=B55 'Latitude (0-90)
Symbol N_S=BIT0 'North=0 & South=1
Symbol LONG=B54 'Longitude (0-180)
Symbol E_W=BIT1 'East=0 & West=1
Symbol MONTH=B53 'Month (1-12)
Click to expand...
Hi, All of this code was developed using Version 5.5.5 of the Picaxe Programming Editor. The syntax of all my code checks out OK in that version. However, when I check the syntax using Version 6.0.8 I see the parsing errors you have mentioned. Does anyone know if there was a change in the allowed formats for these variable assignments? Perhaps someone knows the fix... Thanks!

Michael
 
lbenson

lbenson

Senior Member
Are you using the "#picaxe 20X2" directive?

If you are, and you still have errors, post your entire program, or at least enough to reveal the problem if someone copies it into PE6. There's nothing evidently wrong with your code.
 
B

binary1248

Senior Member
mmruzek: I did copy and run your QB version and I am impressed, especially with the heliostat data. It ran without a single hiccup. Thanks.
Paul
 
lbenson

lbenson

Senior Member
The entire code from that link (450 lines), with "#picaxe 20x2" in front of it, compiled with no problems for me in PE6--"1303 out of 4096 bytes".

Exactly what errors do you get?
 
W

westaust55

Moderator
Only the PICAXE X2 parts have byte variable b0 to b55.
M2 and X1 parts only have b0 to b27 and other parts only up to b13.

My guess is that some are trying to check the syntax without setting the PICAXE type correctly first.
Thus the line:

Symbol LAT=B55 'Latitude (0-90)

with a byte variable b55 exceeds the pre-defined variables availble and an error message will be seen:
Error: Unknown symbol - B55
 
A

AllyCat

Senior Member
Hi Michael,

Yes welcome to the forum and thanks for posting that information, which also compiles and simulates fine for me, in 5.5.6, AxePad and 6.0.8.4 (not updated since I rarely use PE6).

It reminded me that I wrote a similar program some time ago for an 08M2, to calculate sunrise/sunset times and I have been intending to extend it to sun elevation/azimuth predictions for some time. So I put in my own latitude (51 degrees) and sunrise/set times (from the Cumulus Weather station software) into your code for today's date and both came up as zero elevation. :cool:

However, I did note that the elevation at sunset remained at zero for a span of almost 20 minutes, which seems rather long even at my higher latutude (the sun gets up to 30 degrees and down again in around 10 hours at this time of year). Perhaps just due to (unsigned) integer truncation rather than rounding?

I was surprised to discover that I wrote my original program more than four years ago, but when I plugged in today's date (with no other data updates of any type) it still has sunset within one minute and sunrise within 2 minutes of the "official" times (and today must be close to the maximum EOT offset). ;) For my calculations I used two's complement maths, and with an M2 had no choice but to calculate my own trig functions:

Actually, that may have been an advantage because I could use my own angular format and chose a full word (or byte) to represent one complete revolution (i.e. 360 degrees or 2*PI radians). This has the great benefit that there is no need to test for angular numerical overflows, they just wrap around. Combined with the 2's complement format, the whole program (~1000 bytes) needed only 10 IF..THENs , and 7 of those are in the subroutines for signed multiplication/division and of course the SIN/COS/ATAN functions.

Cheers, Alan.
 
B

binary1248

Senior Member
Found my problem. If I download it from his referenced file location (sunweb.bas), it works great. I was trying to copy and past from his listing, don't know why, but now I can load it and play.
You can copy and paste the text below into the Programming Editor. If you prefer, here is the file for download: ( sunweb.bas )
Click to expand...
 
W

westaust55

Moderator
binary1248 said:
Found my problem. If I download it from his referenced file location (sunweb.bas), it works great. I was trying to copy and past from his listing, don't know why, but now I can load it and play.
Click to expand...
Maybe some embedded formatting characters that cause the PE grief.
 
M

mmruzek

New Member
AllyCat said:
Actually, that may have been an advantage because I could use my own angular format and chose a full word (or byte) to represent one complete revolution (i.e. 360 degrees or 2*PI radians). This has the great benefit that there is no need to test for angular numerical overflows, they just wrap around. Combined with the 2's complement format, the whole program (~1000 bytes) needed only 10 IF..THENs , and 7 of those are in the subroutines for signed multiplication/division and of course the SIN/COS/ATAN functions.
Click to expand...
Hello Alan, The code and methods you describe for doing the Sunrise/Sunset time calculation sound really interesting. Would you be able to post the code? Thanks! Michael
 
A

AllyCat

Senior Member
Hi,

Here's the code that I pulled out of my archive; I'm afraid I can only offer it "as is" (E&OE, etc.). To be honest, I didn't fully understand all the variable symbol names when I converted the code, and now four years later I'm not sure what some of the more cryptic comments meant. ;)

I do remember that the source "NOAA Solrad" Excel spreadsheet (on which many other similar programs appear to be based) was packed with macros to convert from the native angular radians to (user) degrees, and to support mathematical under- and over-flows. That's why I adopted the 360 degrees := 65535 / 65536 "fractional binary" numerical conversion. The MSbit of most words or bytes is the sign (2's complement), but to keep the print formatting simple, most data is shown unsigned (e.g. 255 = -1) .

As I live almost exactly on the Greenwich Meridian, I didn't worry too much about the Longitude or Time Zone. IIRC the Longitude can be used relative to the local time zone (e.g. 60 degrees West for a TZ = GMT - 4 hours). Ideally, the Longitude should be changed to a word, to support +/- 180 degrees and a resolution of at least 0.1 degree (which still represents almost half a minute in time). Note that the (dawn/dusk) sun elevation is set to -1 degree (255) because the upper edge (not centre) of the sun is normally considered to determine daylight. So, strictly the Elevation and Latitude also need to be word rather than byte variables.

The subroutines should be (more) understandable; The Trig functions use a 16-line (interpolated) lookup table in EEPROM. A "Universal" interpolation subroutine to handle SIN/COS, ATAN2, Log/Exponential and NTC (Thermistor) curves, etc. has been on my "todo" list for some time, to post in the Snippetts section. :rolleyes:

Code: 
; Sunrise, Solar Noon and Sunset calculation for 08M2, etc.
; AllyCat  15th Feb 2013
; Minor updates 17/2/17
;** = Formulae (adapted) from NOAA Solrad Excel spreadsheet

#picaxe 08m2					; And most others
;#no_data
;#define debug					; About 100 bytes
#define report					; About 130 bytes

symbol EDat  = 0			; Start of data in EEPROM
symbol DayP  = EDat		; Pointer to Day (Binary 1-31) in EPROM 
symbol MonP  = EDat + 1 		; Pointer to Month (Binary 1-12) in EPROM
symbol YearP = EDat + 2		; Pointer to Year (Binary referenced to year 2000) in EPROM
symbol LatP  = EDat + 3		; Latitude of site (degrees, South neg) in EPROM
symbol LongP = EDat + 4		; Longitude of site relative to timezone (degrees,E+) in EPROM 
symbol AltP  = EDat + 5		; Required elevation of Sun (Typ ~ -1 degree for dawn/dusk)
symbol Day = w5			; Daynumber within 4 year cycle
symbol LG1 = w6
symbol G = w7
symbol L = w8 
symbol Lambda = w9
symbol E = w10
symbol GHA = w10			; Greenwich Hour Angle
symbol delta = w8
symbol cost = w7
symbol cost1 = w6
symbol Altit = w9
symbol tempw = w2	
symbol t = w7
symbol tempw1 = w11	; s_w1		;)
symbol tempw2 = w12	; s_w2		;} Used in Trig routines
symbol tempw3 = w13	; s_w3		;)
symbol temp = b3
symbol temp2 = b2
symbol noon = w6
symbol pass  = temp2			; Loop counter
symbol MDat = 8				; Month data

data EDat,(18,2,17,51,0,255)	; Initial Data=DD,MM,YY,Lat,Long(rel TZ),Sun Alt (2's comp deg)
Report:
	read Edat,b4,b5,b6,b7,b8,b9
#ifdef report	
	sertxd("Date:DMY=",#b4,"/",#b5,"/",#b6,lf)
	sertxd("Site Lat=",#b7," Long rel TZ =",#b8,lf)
	sertxd("Sun Elev=",#b9,lf)
#endif
data MDat,(31,28,31,30,31,30,31,31,30,31,30,31)	; Days each month
	Day = 0					; Or from previous years
	bptr = MDAt
	pass = 1					; Month counter
	do while b5 > pass				; Month number	 
		read bptr,temp
		Day = Day + temp
		inc bptr 
		inc pass
	loop
	Day = b6 // 4 * 365 + Day + b4; Include  year and day in current month
	If Day < 61 and b5 < 3 then
		Day = Day - 1		; 1st March in leap year is Day 60 (ref 1 Jan = 0)  
	endif	
#ifdef report
sertxd("Day=",#Day,cr,lf)
#endif
; Convert to part of a revolution	; Multiply by 179.428 (65536^2 per 4yr cycle)
;** LG1 = MOD(360 * Day / 365.25)	; 35.72895277
	LG1 = Day ** 28035			; Fractional part of each day (0.428/65536)
	LG1 = Day * 179 + LG1 			; Subtract 1 days (ref zero is noon on Date 1) 
;sertxd(#LG1,"=LG1 ")	
;**	L = MOD(280.56 + LG1,360)	 
	L = 51074 + LG1				; 280.56 / 360
;**	G = MOD(357.41 + LG1,360)
	G = 65064 + LG1				; 357.41 / 360
	w0 = G					; L + 2 * DegSin(G)
#ifdef debug
sertxd("L=",#L," G=",#G)
#endif
	call sine	
;**	Lambda = MOD(L + 2 * DegSin(G),360)	
	Tempw = 349				; 1.915 scaled to 180 degrees
	call smult		
	Lambda = w0				; sinG coefficient
	w0 = G + G
	call sine
	Tempw = 4 				; 0.02 scaled (unnecessary?)
	call smult
	Lambda = Lambda + w0
	E = -Lambda				; Used in next formula
	Lambda = Lambda + L	
#ifdef debug
sertxd(" Lambda=",#Lambda)		
#endif	
;**	E1 = -1.9 * DegSin(G)			; Done within Lambda	
;**	E2 = 2.5 * DegSin(2 * lambda)	
;**	E3 = 0.053 * DegSin(4 * lambda)	
;**	E = E1 + E2 - E3			
	w0 = Lambda + Lambda			; "2.5 * DegSin(2 * lambda)"
	call sine
	tempw = 449				; 2.5 / fullscale (180)   
	call smult
	E = E + w0
	w0 = Lambda * 4
	call sine
	Tempw = 9				; 0.053 fullscale 180
	call smult 
	GHA = E - w0
#ifdef debug					; Final value
sertxd(" GHA= ",#GHA," ")
#endif
;**	delt1	= DegSin(obliq)
;**	delt2	= delt1 * DegSin(lambda)
;**	delta	= DegArcsin(delt2)	
;**	obliq	= 23.43			
	w0 = 4267				; Obliqiuity scaled (65536 * 23.43 / 360)
	call sine
	Delta = w0				; Start of calculation
	w0 = Lambda
	call sine
	tempw = Delta
	call smult	
	call arcsin
	Delta = w0
#ifdef debug			
sertxd(" Delta=",#w0," ")
#endif
;**	cost1	= DegCos(Lat) * DegCos(delta)
;**	cost2	= DegSin(Lat) * DegSin(delta)
	temp = LatP				; Get site Latitude
	call readdeg				; Start cos(Lat) * cos(delta)	
	call cosine
	cost = w0				; Store to use in numerator
	w0 = delta
	call cosine
	tempw = cost
	call smult					; +1 > denominator > -1   (+ or - sign)
	cost1 = w0
'	temp = LatP				; Still stored ?
	call readdeg				; Lat pointer still stored
	call sine					; Now start sin(Lat) * sin(delta)
	cost = w0				; Might be tempw (if not used in trig functions)
	w0 = delta
	call sine
	tempw = cost
	call smult					; +1 > numerator > -1 in w0
	cost = w0	
;**	cost3	= DegSin(Alt) - cost2			
;**	cost	= cost3 / cost1	
	temp = AltP				; 90 degs -> $4000
	call readdeg				; Get Sun Elevation
	call sine
	w0 = w0 - cost				; Numerator
	cost = w0
#ifdef debug
sertxd(" Num=",#w0)
#endif
	w0 = cost1
	tempw = $ffff			; Now calc reciprocal of denominator (avoids signed div)
	if w0 > $7fff then
		w0 = -w0
		call sdiv
		w0 = -w0
	else
		call	 sdiv
	endif
#ifdef debug	
sertxd(" Div=",#w0)				; MAYBE GETTING RATHER LARGE FOR HIGH LATS
#endif
	tempw = cost + cost			; Numerator 
	call smult
#ifdef debug
sertxd(" Num/Div=",#w0)
#endif
	if w0 > $3fff AND w0 < $c000 then
		sertxd("No sunrise/set",cr,lf)	
	endif
	w0 = w0 + w0				; Signed, OK provided mod < $4000
	call	arccos
	t = w0
	temp = LongP				; Get LSite Longitude
	call readdeg 
	Noon = $8000 - GHA - w0
	w0 = Noon - t
sertxd(cr,lf,"SunRise=")
call printime
#ifdef report
	w0 = noon
	sertxd("Noon=")
	call printime
	w0 = Noon + t
	sertxd("SunSet=")
	call printime
#endif
end	

; SUBROUTINES:
;=============
My text has taken this post over the forum limit, so I have removed the subroutines and will put them in the next post.

Cheers, Alan.
 
Last edited:
A

AllyCat

Senior Member
Hi,

And here are the subroutines. I now notice that the program doesn't actually need ATAN and signed division routines so they are not included here.

Code: 
; SUBROUTINES:
;=============
readdeg:		; Convert degrees (signed) to fractional binary revolution
	tempw = 91 * 256		; 90 degs -> $4000 (change to 90.00 degs?)
	read temp,b1		; Read into high byte of w0
	b0 = 0
	call smult		
return

sdiv:				; Scaled division
	w0 = cost1 / 64		; Scale so typical value (0.5) just fits into lower byte
	w0 = tempw / w0
	w0 = w0 * 64
return
	
smult:				; Signed (fractional word) multiplication: w0 = w0 ** tempw
	if tempw > $7fff then
		tempw = -tempw
		if w0 < $8000 then 
			goto negmult
		else
			w0 = -w0
posmult:
			w0 = w0 ** tempw * 2	; Half-scales are +/-1 and +/- 180 degrees	
			return
		endif
	else
		if w0 < $8000 then
			goto posmult
		else
			w0 = - w0
		endif
	endif		
negmult:
	w0 = w0 ** tempw * 2
	w0 = -w0
return	
		
symbol width = $4000 / 16			; Piecewise Linear horizontal segment length
symbol sintab = 20				; 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,6,128) ; Sine table ($7FFF = 1)

cosine:				; Enter with angle in w0 (2^16 = 360 degrees)
	w0 = w0 + $4000		; Now fall through into sine 
sine:				; Enter with angle in w0 (2^16 = 360 degrees)

	tempw1 = w0
	if tempw1 > $7fff then			; Shift down if 180 - 360 degrees
		tempw1 = tempw1 - $8000  
	endif
	if tempw1 > $3fff then			; Reverse 2nd and 4th quadrants
		tempw1 = $8000 - tempw1
	endif
	
	b0 = tempw1 / width * 2 + sintab		; Segment number + table offset
	read b0,word tempw2,word tempw3		; Read LHS value

	b0 = tempw1 / width 			; Horizontal start of segment
	tempw3 = tempw3 - tempw2 * 8		; Height of segment (scaled)
	tempw1 = -b0 * width + tempw1 * 8		; Horizontal element (scaled)
	tempw1 = tempw1 ** tempw3		; High word
	tempw1 = tempw1 + tempw2 max $7fff	; Positive sine result
	if w0 > $7fff then
		tempw1 = -tempw1		
	endif
	w0 = tempw1	 
return

arccos:		; Input -1 to +1 in w0, output 0 - 360 degrees
	tempw1 = w0
	b1 = 2				; Cosine flag
	goto arcin
arcsin:
	tempw1 = w0
	b1 = 0					; Sine flag
arcin:
	if tempw1 > $7fff then
		tempw1 = -tempw1
		inc b1			; Negative input flag
	endif	
	b0 = sintab
aslp:					; Piecewise linear (interpolated) lookup table
	b0 = b0 + 2
	read b0,b26,b27			; RHS
	if tempw3 < tempw1 then goto aslp
	b0 = b0 - 2
	read b0,b24,b25			; LHS
	b0 = b0 - sintab / 2		; Number of skipped elements 
	tempw3 = tempw3 - tempw2 / 16		; Ramp height (max 32k/10) scaled
	tempw1 = tempw1 - tempw2 * 16		; Element height
	tempw1 = tempw1 / tempw3 * 4		; Element width *256 then scaled *4
	tempw2 = b0 * width + tempw1		; Horizontal angle
	on b1 goto possin,negsin,poscos		; Fall through for last entry
negcos:
	w0 = $4000 + tempw2 			; Range 90 to 180 degrees (-cos)
	return
possin:	
	w0 = tempw2				; Range 0 to 90 degrees (+sin)
	return
negsin:
	w0 = -tempw2				; Range 270 - 360 degrees (-sin)
	return
poscos:						; Range 0 to 90 degrees (+cos)
	w0 = $4000 - tempw2
return

Printime:					; Enter with angle (65535 = 24 hours) in w0
	w12 = w0 ** 14400 + 2		; Time of day in tenths of a minute (rounded)
	b26 = w12 / 600			; Hours
	b27 = w12 // 600 / 10		; Minutes
	w12 = w12 // 10			; Tenths of minute
	sertxd(" ",#b26,":",#b27,".",#w12," ")
return
Finally, as others have said, please post your code directly somewhere on this forum; it's so frustrating when useful code is "lost" by a broken link. :(

Cheers, Alan.
 
Last edited:
You must log in or register to reply here.
Top