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


Definition of Elevation and Azimuth
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.

Screenshot: Original program in action
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.

Graphical depiction Equation of Time (EOT)
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