# GPS Solar Tracker (2)

## MTM Scientific, Inc.

This series of webpages describes a GPS enabled solar tracker based upon the STMAX tilt-tilt platform design. The STMAX solar tracker controller is described in much more detail in our book "Build a Solar Tracker" and at http://www.mtmscientific.com/stmax.html.

The first webpage of this series described how to acquire time and position data from a GPS module. This webpage describes how to use the acquired GPS information to calculate the position of the sun in the sky. (Our final webpage on this topic describes how to use the MEMS ADC counts to aim the tracker: How to Aim. )The calculation method used here is an adaptation of a software program written by David Williams. The original program was written for generic BASIC, and calculates the Elevation and Azimuth of the sun based upon the observer's position and time. We find it quite remarkable that such a relatively simple program is able to perform this calculation so well. Here is the program souce code: sunalign.bas Figure 1. Definition of Elevation and Azimuth

You can run the original program using many different versions of BASIC, for example Microsoft's QuickBASIC Version 4.50. The original program is interactive, and will ask for the time, date and space coordinates it needs. We suggest that you familiarize yourself with the original program to gain a basic understanding of how it works. Figure 2. Screenshot showing original program being used in QuickBASIC

There are several substantial programming challenges involved in adapting the original program to Picaxe BASIC, used by the Picaxe-20X2 microcontroller. The main challenges are: 1) Picaxe BASIC performs all math calculations using only positive single or double byte integers, 2) The Picaxe BASIC trigonometric functions use units of degrees (instead of radians), and 3) The devices using Picaxe BASIC have limited memory space for variables and data. Given these limitations you might be surprised to know this calculation is even possible using a Picaxe-20X2!

Since Picaxe BASIC only performs positive integer math, it is necessary to write code with each variable assigned a bit flag to indicate the sign. Subsequent math operations use the flags in logical branches for calculations. In essence we are creating a method to handle negative and positive integer numbers. The fact that we must perform math using integers requires keeping the arithmetic range between 0 and 255 for a byte variable, or 0 and 65535 for  a word variable.

The limited memory space of devices used with Picaxe BASIC requires the careful assignment and use of variables. The single most important factor in that regard is to use (and reuse) scratchpad variables as often as possible during the calculations. We emphasize this is in the code by often naming variables as temporary.

Because negative integers are not allowed in Picaxe BASIC we use a method of addition and substraction in several places in the code called 'midpoint math'. Essentially this is a simple technique of shifting the zero origin to the midpoint of a word variable integer (32768).

We also found an opportunity to improve the accuracy of the Picaxe BASIC sine function near a zero crossing. Picaxe BASIC uses a unique coding system to indicate positive and negative trigonometric values, refer to the manuals for more details.  Improving the accuracy of the trigonometric functions in general might be an opportunity for future work.

The calculation of sun position also requires calculation of an entity called the 'Equation of Time' (EOT). Essentially this is the offset in minutes relating the sun's actual position relative to an average theoretical position. We have modified the original code to use a different method of calculating the EOT. See the comments in the code for a web reference. Figure 3. Graph of Equation of Time (EOT)

This code can be easily tested using the simulation mode of the Picaxe Programming Editor. (You don't even need a Picaxe-20X2 chip!) We've done some basic testing of the code for different times, dates and locations. You may want to do more testing for your location. The setup parameters for the input variables are assigned at the start of the code, and can easily be changed as desired.

Here is the Picaxe BASIC code for calculating the position of the sun in the sky based on the observer's location and time. You can copy and paste the text below into the Programming Editor. If you prefer, here is the file for download: ( sunweb.bas )

'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)
Symbol DAY=B52      'Day of Month (1-31)
Symbol HOUR=B51   'Hour of Day UTC (0-23)
Symbol MINUTE=B50 'Minute of Day UTC (0-59)
'Above variables will all eventually be parsed from GPS
Symbol DAYNO=W18   'Day of Year (1-365)
Symbol BORG=W17      'Symbol "B" original program
Symbol ELEV=W2         'Elevation, Positive from Horizon
'Above variables important enough to name
Symbol TEMP1=W24    'Scratchpad (0-65535)
Symbol TEMP2=W23    'Scratchpad (0-65535)
Symbol TEMP3=W22    'Scratchpad (0-65535)
Symbol TEMP4=W21    'Scratchpad (0-65535)
Symbol TEMP5=W20    'Scratchpad (0-65535)
Symbol TEMP6=W19    'Scratchpad (0-65535)
Symbol TEMP7=W16    'Scratchpad (0-65535)
Symbol TEMP8=W15    'Scratchpad (0-65535)
Symbol TEMP9=W14    'Scratchpad (0-65535)
Symbol TEMP10=W13  'Scratchpad (0-65535)
Symbol TEMP11=W12  'Scratchpad (0-65525)
Symbol TEMP12=W11  'Scratchpad (0-65525)
Symbol TEMP13=W10  'Scratchpad (0-65525)
Symbol TEMP14=W9    'Scratchpad (0-65525)
Symbol TEMP15=W8    'Scratchpad (0-65525)
Symbol TEMP16=W7    'Scratchpad (0-65525)
Symbol TEMP17=W6    'Scratchpad (0-65525)
Symbol TEMP18=W5    'Scratchpad (0-65525)
'Above variables are scratchpad variables, bit flags are here
Symbol SSIGN=BIT2      'SINE SIGN, Negative=1, Positive=0
Symbol CSIGN=BIT3      'COSINE SIGN, Negative=1, Positive=0
Symbol ESIGN=BIT4      'EOT SIGN, Negative=1, Positive=0
Symbol S2SIGN=BIT5    'SINE 2 SIGN, Negative=1, Positive=0
Symbol T4SIGN=BIT6    'TEMP4 Sign, Negative=1, Positive=0
Symbol T3SIGN=BIT7    'TEMP3 Sign, Negative=1, Positive=0
Symbol T5SIGN=BIT8    'TEMP5 Sign, Negative=1, Positive=0
Symbol T7SIGN=BIT9    'TEMP7 Sign, Negative=1, Positive=0
Symbol T8SIGN=BIT10    'TEMP8 Sign, Negative=1, Postive=0
Symbol T12SIGN=BIT11  'TEMP12 Sign, Negative=1, Positive=0
Symbol T13SIGN=BIT12  'TEMP13 Sign, Negative=1, Positive=0
Symbol T14SIGN=BIT13  'TEMP14 Sign, Negative=1, Positive=0
Symbol T15SIGN=BIT14  'TEMP15 Sign, Negative=1, Positive=0
Symbol T19SIGN=BIT15  'TEMP19 Sign, Negative=1, Positive=0
Symbol T20SIGN=BIT16  'TEMP20 Sign, Negative=1, Positive=0
Symbol T6SIGN=BIT18    'TEMP6 Sign, Negative=1, Positive=0
Symbol T9SIGN=BIT19    'TEMP9 Sign, Negative=1, Positive=0
Symbol T10SIGN=BIT20   'TEMP10 Sign, Negative=1, Positive=0
Symbol UP=BIT21             'UP=1=TRUE=SUN_UP, UP=0=FALSE
Symbol T1SIGN=BIT22    'TEMP1 Sign, Negative=1,Positive=0
Symbol T2SIGN=BIT23    'TEMP2 Sign, Negative=1,Positive=0
'Note: Must preserve N/S (TEMP4) and E/W (TEMP7)
'****************************************************************
'Use this section to enter data for any Sun position calculation
'Example: 42N, 84W, 10/12/16, 16:00 UTC is 153.8 AZ, 37.4 ELEV
'Serial Out is 9600, N, 8, 1 (Picaxe default settings)
LAT=42
N_S=0
LONG=84
E_W=1
MONTH=10
DAY=12
HOUR=16
MINUTE=0
'****************************************************************
'Begin Sun sky position calculation
'Calculate day number based on date
TEMP1=MONTH + 9
TEMP1=TEMP1 % 12    'Modulus operator returns remainder
TEMP1=TEMP1 * 306    'Do math 10X scale for accuracy
TEMP2=DAY * 10          'Do math 10X scale for accuracy
TEMP1=TEMP1 + 585 + TEMP2
TEMP1=TEMP1 / 10       'Remove the 10X scale
TEMP1=TEMP1 % 365   'Day number of year
DAYNO=TEMP1
SERTXD("DAYNO=", #DAYNO,32,13,10)

'A and B must be calculated
TEMP1=TEMP1 + 10      'Solstice is 10 days before January 1.
TEMP2=TEMP1 / 73       'Correction factor (360 deg in 365 days)
TEMP1=TEMP1-TEMP2 'Degrees around orbit since Solstice
TEMP3=TEMP1              'Equal to 'A' original program
IF TEMP1<12 THEN      'Subtraction precaution
TEMP1=TEMP1+360      'Subtraction precaution
ENDIF        'Finish
TEMP2=TEMP1-12       'Correction for solstice to perihelion
TEMP1=SIN TEMP2     'Sine Function with 100X scale
IF TEMP1>127 THEN  'See Picaxe SIN function
SSIGN=1        'Store sign information
TEMP1=TEMP1-128    'Remove coding system
ENDIF        'Finish
TEMP1=TEMP1 / 30    'Identical to multiplying by 0.0334
TEMP1=TEMP1 / 100  'Remove 100X scale
IF SSIGN=0 THEN
TEMP2=TEMP3+TEMP1
ELSE
TEMP2=TEMP3-TEMP1
ENDIF        'Equal to 'B' original program in degrees
BORG=TEMP2
'End of A and B calculation
'**************************************************************
'Calculate Equation of Time correction from date
'See http://mb-soft.com/public3/equatime.html (near bottom)
IF DAYNO<81 THEN
DAYNO=DAYNO+365
ENDIF
TEMP2=DAYNO-81
TEMP2=TEMP2*72    'Same as 360/365
TEMP2=TEMP2/73    'TEMP2 equal "B"
TEMP3=SIN TEMP2
IF TEMP3>100 THEN
SSIGN=1
TEMP3=TEMP3-128
ENDIF
TEMP3=3*TEMP3    'Same as 1.5(*3 & /2)
TEMP3=TEMP3/2     'Last term of equation
TEMP4=COS TEMP2
IF TEMP4>100 THEN
CSIGN=1
TEMP4=TEMP4-128
ENDIF
TEMP4=15*TEMP4    'Close to 7.53(*15 & /2)
TEMP4=TEMP4/2    'Second to last term of equation
TEMP5=2*TEMP2
IF TEMP5>360 THEN
TEMP5=TEMP5-360
ENDIF
TEMP5=SIN TEMP5
IF TEMP5>100 THEN
S2SIGN=1
TEMP5=TEMP5-128
ENDIF
TEMP5=TEMP5*69    'Close to 9.87
TEMP5=TEMP5/7

TEMP6=32768    'Midpoint value for 16 bit word math
IF S2SIGN=1 THEN
TEMP6=TEMP6-TEMP5
ELSE
TEMP6=TEMP6+TEMP5
ENDIF
IF CSIGN=1 THEN    'Negative sign in equation!
TEMP6=TEMP6+TEMP4
ELSE
TEMP6=TEMP6-TEMP5
ENDIF
IF SSIGN=1 THEN    'Negative sign in equation!
TEMP6=TEMP6+TEMP3
ELSE
TEMP6=TEMP6-TEMP3
ENDIF
IF TEMP6>32767 THEN    'Remove midpoint math
ESIGN=0
TEMP6=TEMP6-32768
ELSE
ESIGN=1
TEMP6=32768-TEMP6
ENDIF
TEMP6=TEMP6/100    'Remove the 100X scale
SERTXD("EOT=", #TEMP6,32,13,10)
SERTXD("EOT_SGN=", #ESIGN,32,13,10)
'*************************************************************
'Back to method in the original D. Williams program...
TEMP1=COS BORG
IF TEMP1>100 THEN
SSIGN=1
TEMP1=TEMP1-128
ELSE
SSIGN=0
ENDIF
TEMP1=TEMP1*2    'Close to -.3939
TEMP1=TEMP1/5    '(*2 & /5)
IF SSIGN=1 THEN
SSIGN=0
ELSE
SSIGN=1
ENDIF
TEMP2=TEMP1*TEMP1    'Squared: Always positive
TEMP2=10000-TEMP2       '10000=100*100 (Scale Factor)
TEMP2=SQR TEMP2         'Square root with result 100X scaled
TEMP1=TEMP1*100
TEMP2=TEMP1/TEMP2
T2SIGN=SSIGN
IF TEMP2<101 THEN
TEMP2=ATAN TEMP2
ELSE
TEMP2=10000/TEMP2
TEMP2=ATAN TEMP2
TEMP2=90-TEMP2
ENDIF
'Take SSIGN as the sign for this angular result in degrees
'Now calculate the term AZ in original
TEMP3=32768    '      'Midpoint math
IF ESIGN=0 THEN    'EOT is positive
TEMP4=MINUTE+TEMP6
ELSE            'EOT is negative
IF MINUTE>TEMP6 THEN
TEMP4=MINUTE-TEMP6
T4SIGN=0
ELSE
TEMP4=TEMP6-MINUTE
T4SIGN=1
ENDIF
ENDIF
TEMP4=TEMP4/4    'Has same sign as T4SIGN
IF T4SIGN=0 THEN
TEMP3=TEMP3+TEMP4
ELSE
TEMP3=TEMP3-TEMP4
ENDIF
IF E_W=0 THEN    'Sum the longitude
TEMP3=TEMP3+LONG
ELSE
TEMP3=TEMP3-LONG
ENDIF
TEMP5=15*HOUR
TEMP3=TEMP3+TEMP5    'Sum the hour
IF TEMP3>32767 THEN
TEMP3=TEMP3-32768
T3SIGN=0
ELSE
TEMP3=32768-TEMP3
T3SIGN=1
ENDIF
'**********************************************************
'**********************************************************
'Calculate GOSUB 800 from original
'EL is TEMP2 with sign SSIGN & AZ is TEMP3 with sign T3SIGN
'Perform 2 calculations with EL first
IF SSIGN=1 THEN
TEMP2=360-TEMP2
ENDIF
TEMP4=SIN TEMP2
IF TEMP4>100 THEN
TEMP4=TEMP4-128
T4SIGN=1
ELSE
T4SIGN=0
ENDIF
'TEMP4 equals 'Z' with factor of 100 and sign T4SIGN
'Calculate COS term in GOSUB 800
'Already did the 360 correction check above
TEMP5=COS TEMP2
IF TEMP5>100 THEN
TEMP5=TEMP5-128
T5SIGN=1
ELSE
T5SIGN=0
ENDIF
IF T5SIGN=1 THEN
T5SIGN=0
ELSE
T5SIGN=1
ENDIF
'Perform 2 calculation with AZ next
IF T3SIGN=1 THEN
TEMP3=360-TEMP3
ENDIF
TEMP7=SIN TEMP3
IF TEMP7>100 THEN
TEMP7=TEMP7-128
T7SIGN=1
ELSE
T7SIGN=0
ENDIF
'Improves Sine function near zero crossing
IF TEMP3>180 AND TEMP3<200 THEN
TEMP7=TEMP3-180
TEMP7=18 * TEMP7
TEMP7=TEMP7/10
T7SIGN=1
ENDIF
IF TEMP3>160 AND TEMP3<180 THEN
TEMP7=180-TEMP3
TEMP7=18*TEMP7
TEMP7=TEMP7/10
T7SIGN=0
ENDIF
TEMP7=TEMP5*TEMP7
TEMP7=TEMP7/100    'Still has a factor of 100X
'Compare T7SIGN and T5SIGN
IF T7SIGN=T5SIGN THEN
T7SIGN=0
ELSE
T7SIGN=1
ENDIF
'**************************************************
'Do the final calculation with AZ
TEMP8=COS TEMP3
IF TEMP8>100 THEN
TEMP8=TEMP8-128
T8SIGN=1
ELSE
T8SIGN=0
ENDIF
TEMP8=TEMP5*TEMP8
TEMP8=TEMP8/100    'Still has a factor of 100X
'Compare T8SIGN and T5SIGN
IF T8SIGN=T5SIGN THEN
T8SIGN=0
ELSE
T8SIGN=1
ENDIF
'This is end of the GOSUB 800 calculations
'*****************************************************
'*****************************************************
TEMP9=TEMP8*TEMP8    'Equals 'Y' squared, Always positive
TEMP10=TEMP4*TEMP4'Equals 'Z' squared, Always positive
TEMP11=TEMP9+TEMP10
TEMP11= SQR TEMP11'Square root with factor 100X
GOSUB SEVENTEN '************LINE 480 in original
'Calculate 'A' original program
'Use LAT with N_S sign with TEMP13 and T13SGN
'Use AA=TEMP12 with sign T12SIGN
'LAT can not be greater than 90
IF N_S=0 THEN
TEMP13=90-LAT
T13SIGN=0
ELSE
TEMP13=90+LAT
T13SIGN=0
ENDIF
'Calculate 'A'
IF T12SIGN=0 THEN
TEMP13=TEMP12+TEMP13
T13SIGN=0
ELSE
IF TEMP12<TEMP13 THEN
TEMP13=TEMP13-TEMP12
T13SIGN=0
ELSE
TEMP13=TEMP12-TEMP13
T13SIGN=1
ENDIF
ENDIF
'Calculate 'Y' and 'Z' with R=TEMP11 positive & 100X
TEMP14=COS TEMP13
IF TEMP14>100 THEN
TEMP14=TEMP14-128
T14SIGN=1
ELSE
T14SIGN=0
ENDIF
TEMP14=TEMP14*TEMP11
TEMP14=TEMP14/100    'Still has 100X factor
IF T13SIGN=1 THEN    'Negative case before dawn...
TEMP13=360-TEMP13    'SIN function repeats 360 degrees
ENDIF            'See Picaxe manual for more info...
TEMP15=SIN TEMP13
IF TEMP15>100 THEN
TEMP15=TEMP15-128
T15SIGN=1
ELSE
T15SIGN=0
ENDIF
'Improves Sine function near zero crossing at 180
IF TEMP13>180 AND TEMP13<200 THEN
TEMP3=TEMP3-180
TEMP15=17 * TEMP13
TEMP15=TEMP15/10
T15SIGN=1
ENDIF
IF TEMP13>160 AND TEMP13<180 THEN
TEMP13=180-TEMP13
TEMP15=17*TEMP13
TEMP15=TEMP15/10
T15SIGN=0
ENDIF
TEMP15=TEMP15*TEMP11
TEMP15=TEMP15/100    'Still has 100X factor
'****************************************************
'****************************************************
'Calculate GOSUB 740 in original program
'AX=TEMP8 with T8SIGN, AY=TEMP4 with T4SIGN
TEMP16=TEMP14*TEMP14
TEMP16=TEMP16    'Has 10000X factor, reduced by SQR
TEMP17=TEMP7*TEMP7
TEMP17=TEMP17    'Has 10000X factor, reduced by SQR
TEMP17=TEMP16+TEMP17
TEMP8=SQR TEMP17
T8SIGN=0
TEMP4=TEMP15
T4SIGN=T15SIGN
GOSUB SEVENTEN    'Called by the GOSUB 740
'Note: SEVENTEN calculated elevation and sign
SERTXD("SUN ELEVATION=", #TEMP12,32,13,10)
SERTXD("ELEV SIGN=", #T12SIGN,32,13,10)
IF T12SIGN=1 THEN
SERTXD("SUN BELOW HORIZON", 32,13,10)
ENDIF
ELEV=TEMP12        'Store for aiming calculation
'Calculate AX and AZ for second call to GOSUB 710
'Assign TEMP8 and TEMP4 before calling
TEMP8=TEMP14
T8SIGN=T14SIGN
TEMP4=TEMP7
T4SIGN=T7SIGN
GOSUB SEVENTEN 'This is inside GOSUB 740
'***************************************************
'Manage before or after solar noon
IF T12SIGN=0 THEN
TEMP12=180+TEMP12
ELSE
TEMP12=180-TEMP12
ENDIF
SERTXD("SUN AZIMUTH=", #TEMP12,32,13,10)
'Sun position calculation is completely finished!
'This is the program ending.
END_P:
SERTXD("PROGRAM ENDING",32,13,10)
END
'***************************************************
'SUBROUTINES
'***************************************************
SEVENTEN:
'Start calculation of GOSUB 710
'AX=Y=TEMP8 with T8SIGN and AY=Z=TEMP4 with T4SIGN
'Calculates 'AA' in original, Returns angle in degrees
IF TEMP8=0 THEN
TEMP12=9000
T12SIGN=T8SIGN
ELSE
TEMP12=TEMP4*100    'Maximum is 100*100
TEMP12=TEMP12/TEMP8
IF T4SIGN != T8SIGN THEN
T12SIGN = 1
ELSE
T12SIGN=0
ENDIF
'*********************************************
IF TEMP12<101 THEN    'Allowed range of ATAN
TEMP12=ATAN TEMP12
ELSE
TEMP18=10000/TEMP12
TEMP18=ATAN TEMP18
TEMP18=90-TEMP18
TEMP12=TEMP18
ENDIF
'**********************************************
IF T8SIGN=1 THEN
TEMP12=TEMP12+180
ENDIF
ENDIF
RETURN