Forum - MCS Electronics

 

FAQFAQ SearchSearch RegisterRegister Log inLog in

Sun position and Heliostat offset.

 
Post new topic   Reply to topic    www.mcselec.com Forum Index -> Share your working BASCOM-AVR code here
View previous topic :: View next topic  
Author Message
Dave

Bascom Member



Joined: 05 Feb 2005
Posts: 264
Location: McMinnville, OR

blank.gif
PostPosted: Sun Jul 14, 2019 6:02 pm    Post subject: Sun position and Heliostat offset. Reply with quote

Here is a nice little program for calculating the suns position and also heliostat position.
I translated it from a QBasic program that can be found here.
https://www.instructables.com/community/heliostat-and-sun-tracker-basic-program/
Thanks to David Williams for original QBasic code.



Code:
'Thanks to David Williams for sharing the QBasic version of this program.
'The original QBasic code can be found here:
'https://www.instructables.com/community/heliostat-and-sun-tracker-basic-program/

' SunAlign.BAS (Version for QBasic and similar dialects)

' Calculates position of sun in sky, as azimuth (compass bearing
' measured clockwise from True North) and angle of elevation, as
' seen from any place on earth, on any date and any time.
' Also calculates alignment of a heliostat mirror.

' David Williams

' Initially dated 2007 Jul 07
' This version 2008 Jan 13

' All angles in radians except in i/o routines DegIn and DegOut

'Translated to BASCOM by:
'Dave Carr
'July 14, 2019

$regfile = "m328pdef.dat"
$crystal = 8000000
$baud = 38400
$hwstack = 100
$swstack = 100
$framesize = 100

DECLARE SUB C2P(byref X as single, byref Y as single, byref Z as single, byref AZ as single, byref EL as single)
DECLARE SUB P2C(byval AZ as single, byval EL as single, byref X as single ,byref Y as single ,byref Z as single)
DECLARE FUNCTION Ang(byval X as single ,byval Y as single)as single
DECLARE SUB DegIn(byval P as string,byref X as single)
DECLARE SUB DegOut(byval P as string,byref X as single)

CONST c_PY = 3.1415926536 ' "PI" not assignable in some BASICs
CONST c_DR = 180 / c_PY ' degree / radian factor
const c_W = 2 * c_PY / 365 ' earth's mean orbital angular speed in radians/day
const c_WR = c_PY / 12' earth's speed of rotation relative to sun (radians/hour)
const c_C = -23.45 / c_DR ' reverse angle of earth's axial tilt in radians
const c_ST = SIN(c_C) ' sine of reverse tilt
const c_CT = COS(c_C) ' cosine of reverse tilt
const c_E2 = 2 * .0167 ' twice earth's orbital eccentricity
const c_SN = 10 * c_W ' 10 days from December solstice to New Year (Jan 1)
const c_SP = 12 * c_W ' 12 days from December solstice to perihelion

dim s_C as single
dim st_S as string * 1
dim s_Lt as single
dim s_LG as single
dim i_TZN as integer
dim b_HR, b_MIN as byte
dim b_Mth, b_Day as byte
dim s_CL as single
dim w_D as word
dim s_A, s_B, s_L, s_DC, s_LD, s_SL as single
dim s_tAZ, s_tEL, s_tX, s_tY, s_tZ, s_mAZ, s_mEL as single
dim s_sX, s_sY, s_sZ, s_sAZ, s_sEL as single
dim s_temp1, s_temp2, s_temp3 as single
CLS

Menu:

   PRINT "1. Calculate sun's position"
   PRINT "2. Calculate mirror orientation"
   PRINT "3. Calculate both"
   PRINT "4. Quit program"
   PRINT
   PRINT "Which? (1 - 4)";

   DO
      st_S = INKEY()
   LOOP UNTIL st_S >= "1" AND st_S <= "4"
   PRINT st_S

   IF st_S ="4" THEN END

' Note: For brevity, no error checks on user inputs

   PRINT
   PRINT "Use negative numbers for directions opposite to those shown."
   PRINT
   DegIn "Observer's latitude (degrees North)", s_LT
   DegIn "Observer's longitude (degrees East)", s_LG
   INPUT "Time Zone (+/- hours from GMT/UT)", i_TZN
   INPUT "Hours (24-hr format)", b_HR
   input "Minutes ", b_MIN
   INPUT "Month ", b_Mth
   input "Day of month ", b_Day

   PRINT

   's_CL = c_PY / 2 - s_LT ' co-latitude
   s_CL = c_PY / 2
   s_CL = s_CL - s_LT
   'print "co-lat ";s_CL

   'w_D = INT(30.6 * ((b_Mth + 9) MOD 12) + 58.5 + b_Day) MOD 365
   ' day of year (D = 0 on Jan 1)
   s_C = b_Mth + 9
   s_C = s_C mod 12
   s_C = 30.6 * s_C
   s_C = s_C + 58.5
   s_C = s_C + b_Day
   s_C = int(s_C)
   w_D = s_C mod 365
   'print "Day of year = ";w_D

   's_A = c_W * w_D + c_SN ' orbit angle since solstice at mean speed
   s_A = c_W * w_D
   s_A = s_A + c_SN
   'print "Orbit angle ";s_A

   's_B = s_A + c_E2 * SIN(s_A - c_SP) ' angle with correction for eccentricity
   s_C = s_A - c_SP
   s_B = sin(s_C)
   s_B = s_B * c_E2
   s_B = s_B + s_A
   'print "Eccentricity ";s_B

   's_C = (s_A - ATN(TAN(s_B) / c_CT)) / c_PY
   s_C = tan(s_B)
   s_C = s_C/c_CT
   s_C = atn(s_C)
   s_C = s_A - s_C
   s_C = s_C/c_PY
   's_SL = c_PY * (s_C - INT(s_C + .5))' solar longitude relative to mean position
   s_SL = s_C + .5
   s_SL = int(s_SL)
   s_SL = s_C - s_SL
   s_SL = c_PY * s_SL
   'print "solar longitude ";s_SL

   's_C = c_ST * COS(s_B)
   s_C = cos(s_B)
   s_C = c_ST * s_C
   's_DC = ATN(s_C / SQR(1 - s_C * s_C)) ' solar declination (latitude)
   ' arcsine of C. ASN not directly available in QBasic
   s_DC = s_C * s_C
   s_DC = 1 - s_DC
   s_DC = sqr(s_DC)
   s_DC = s_C / s_DC
   s_DC = atn(s_DC)
   'print "solar declination ";s_DC

   's_LD = (b_HR - i_TZN + b_MIN / 60) * c_WR + s_SL + s_LG ' longitude difference
   s_LD = b_MIN / 60
   s_C = b_HR - i_TZN
   s_LD = s_C + s_LD
   s_LD = s_LD * c_WR
   s_LD = s_LD + s_SL
   s_LD = s_LD + s_LG
   'print "longitude difference ";s_LD

   CALL P2C(s_LD, s_DC, s_sX, s_sY, s_sZ) ' polar axis (perpend'r to azimuth plane)
   CALL C2P(s_sY, s_sZ, s_sX, s_sAZ, s_sEL) ' horizontal axis
   s_temp1=s_sAZ - s_CL
   CALL P2C(s_temp1, s_sEL, s_sY, s_sZ, s_sX) ' rotate by co-latitude

   IF s_sZ < 0 THEN
      'BEEP
      PRINT "Sun Below Horizon"
      PRINT
      GOTO NewCalc
   END IF

   IF st_S <> "2" THEN ' calculate and display sun's position

      CALL C2P(s_sX, s_sY, s_sZ, s_sAZ, s_sEL) ' vertical axis

      DegOut "Sun's azimuth: ", s_sAZ
      DegOut "Sun's elevation: ", s_sEL

      PRINT

   END IF

   IF st_S > "1" THEN ' calculate and display mirror orientation

      PRINT "For target direction of light reflected from mirror:"
      DegIn "Azimuth of target direction (degrees)", s_tAZ
      DegIn "Elevation of target direction (degrees)", s_tEL

      PRINT

      CALL P2C(s_tAZ, s_tEL, s_tX, s_tY, s_tZ) ' target vector X,Y,Z
      s_temp1=s_sX + s_tX
      s_temp2=s_sY + s_tY
      s_temp3=s_sZ + s_tZ
      CALL C2P(s_temp1,s_temp2 ,s_temp3 , s_mAZ, s_mEL)
' angle bisection by vector addition

      PRINT "Mirror aim direction (perpendicular to surface):"
      DegOut "Azimuth: ", s_mAZ
      DegOut "Elevation: ", s_mEL

      PRINT

   END IF

NewCalc:

   PRINT
   PRINT "New Calculation"
   PRINT

   GOTO Menu

FUNCTION Ang (byval X as single,byval Y as single)
' calculates angle from positive X axis to vector to (X,Y)
   'SELECT CASE SGN(X)
   'CASE 1: Ang = ATN(Y / X)
   'CASE -1: Ang = ATN(Y / X) + PY
   'CASE ELSE: Ang = SGN(Y) * PY / 2
   'END SELECT
   local s as single
   s = sgn(x)
   SELECT CASE s
      CASE 1
         s=y/x
         Ang = ATN(s)
      CASE -1
         s=y/x
         s=atn(s)
         Ang = s+c_PY
      CASE ELSE
         s=sgn(y)
         s=s*c_PY
         Ang = s/2
   END SELECT

END FUNCTION

SUB C2P (byref X as single,byref Y as single,byref Z as single,byref AZ as single,byref EL as single)
' Cartesian to Polar. Convert from X,Y,Z to AZ,EL
   local A as single
   local s as single
   local s2 as single
   'EL = Ang(SQR(X * X + Y * Y), Z)
   s=x*x
   s2=y*y
   s=s+s2
   s=sqr(s)
   el=ang(s,z)
   'A = Ang(Y, X)
   A = Ang(Y, X)
   'IF A < c_PY THEN AZ = A + c_PY ELSE AZ = A - c_PY
   IF A < c_PY THEN
      AZ = A + c_PY
   ELSE
      AZ = A - c_PY
   end if
END SUB

SUB DegIn (byval P as string,byref X as single)
' Input angle in degrees and convert to radians
   local N as single
   PRINT P;
   INPUT N
   X = N / c_DR
END SUB

SUB DegOut (byval P as string,byref X as single)
' converts radians to degrees, rounds to nearest 0.1, and prints
   local st as string * 20
   local s as single
   local i as integer
   'S$ = LTRIM$(STR$(INT(10 * ABS(X * c_DR) + .5)))
   s=X*c_DR
   s=abs(s)
   s=s*10
   s=s + .5
   i=int(s)
   st=str(i)
   st=ltrim(st)
   IF St = "3600" THEN St = "0"
   IF LEN(St) = 1 THEN St = "0" + St
   'IF X < 0 THEN IF VAL(St) THEN St = "-" + St
   IF X < 0 THEN
      s=val(st)
      IF s>0 THEN St = "-" + St
   end if
   'PRINT P; LEFT(St, LEN(St) - 1); "."; RIGHT$(St, 1); " degrees"
   i=len(st)-1
   print P;left(st,i);".";right(st,1);" degrees"
END SUB

SUB P2C (byval AZ as single,byval EL as single,byref X as single,byref Y as single,byref Z as single)
' Polar to Cartesian. Convert from AZ,EL to X,Y,Z
   local s as single
   Z = SIN(EL)
   s_C = COS(EL)
   s_C=s_C*-1
   X = s_C * SIN(AZ)
   Y = s_C * COS(AZ)
END SUB
Back to top
View user's profile
albertsm

Administrator



Joined: 09 Apr 2004
Posts: 4747
Location: Holland

blank.gif
PostPosted: Tue Jul 16, 2019 9:51 am    Post subject: Reply with quote

hi Dave

Thanks for sharing.
I wonder if the accuracy would not be better when using double instead of single?

_________________
Mark
Back to top
View user's profile Visit poster's website
Dave

Bascom Member



Joined: 05 Feb 2005
Posts: 264
Location: McMinnville, OR

blank.gif
PostPosted: Tue Jul 16, 2019 3:07 pm    Post subject: Reply with quote

Hello Mark

Quote:
I wonder if the accuracy would not be better when using double instead of single?


I don't think so. After I first got the code working I ran it in the simulator and compared the output to online calculators. The output was off by a few tenths of a degree so I did a quick search and replace of all singles to doubles and the output was still off by about the same amount. I think David Williams took some shortcuts to keep the code small at the cost of accuracy. But for pointing a solar panel at the sun or reflecting sun light to a shaded part of my garden a few tenths a degree is fine.

I'm working to optimize the code to take advantage on some of bascoms built in features like fusing and asin. I'll share that at a latter date.
Back to top
View user's profile
Duval JP

Bascom Member



Joined: 22 Jun 2004
Posts: 840
Location: France

france.gif
PostPosted: Wed Jul 17, 2019 9:20 am    Post subject: Reply with quote

Hi Dave
I really like your program, my hobby being astronomical clock building, this kind of calculation is always interesting.
To make it usable by everyone, I suggest you to use an SD card and AVR-DOS for storing latitude and longitude coordinates.
As Mark suggests : the use of double give you more accuracy, a little error+ a little error give you a big one.
JP Wink

_________________
pleasure to learn, to teach, to create
Back to top
View user's profile Visit poster's website
albertsm

Administrator



Joined: 09 Apr 2004
Posts: 4747
Location: Holland

blank.gif
PostPosted: Wed Jul 17, 2019 9:33 am    Post subject: Reply with quote

Ok, for a solar panel it is fine Smile

Wish my roof was better pointed towards the Sun.

_________________
Mark
Back to top
View user's profile Visit poster's website
Dave

Bascom Member



Joined: 05 Feb 2005
Posts: 264
Location: McMinnville, OR

blank.gif
PostPosted: Wed Jul 17, 2019 11:32 pm    Post subject: Reply with quote

JP

Thanks for your suggestions.
Have you posted any pics of your astronomical clocks? I'd like to see them.

D.
Back to top
View user's profile
Duval JP

Bascom Member



Joined: 22 Jun 2004
Posts: 840
Location: France

france.gif
PostPosted: Thu Jul 18, 2019 9:29 am    Post subject: Reply with quote

Hi Dave,

see https://www.mcselec.com/index2.php?option=com_forum&Itemid=59&page=viewtopic&t=13096&highlight=astronomical

astro-IV will soon arrive with a simpler method to enter latitude, longitude and birthdays, without or with an Excel sheet

JP

Wink

_________________
pleasure to learn, to teach, to create
Back to top
View user's profile Visit poster's website
Dave

Bascom Member



Joined: 05 Feb 2005
Posts: 264
Location: McMinnville, OR

blank.gif
PostPosted: Thu Jul 18, 2019 4:11 pm    Post subject: Reply with quote

JP,

WOW! Great work! Keep doing what you are doing. I will be looking closely at your code.

Thanks,
Dave
Back to top
View user's profile
Display posts from previous:   
Post new topic   Reply to topic    www.mcselec.com Forum Index -> Share your working BASCOM-AVR code here All times are GMT + 1 Hour
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
You cannot attach files in this forum
You cannot download files in this forum