Forum - MCS Electronics

Author Message
Dave Joined: 05 Feb 2005
Posts: 264
Location: McMinnville, OR  Posted: Sun Jul 14, 2019 6:02 pm    Post subject: Sun position and Heliostat offset. 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  albertsm  Joined: 09 Apr 2004
Posts: 4747
Location: Holland  Posted: Tue Jul 16, 2019 9:51 am    Post subject: hi Dave Thanks for sharing. I wonder if the accuracy would not be better when using double instead of single?_________________Mark   Dave Joined: 05 Feb 2005
Posts: 264
Location: McMinnville, OR  Posted: Tue Jul 16, 2019 3:07 pm    Post subject: 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.  Duval JP  Joined: 22 Jun 2004
Posts: 840
Location: France  Posted: Wed Jul 17, 2019 9:20 am    Post subject: 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 _________________pleasure to learn, to teach, to create   albertsm  Joined: 09 Apr 2004
Posts: 4747
Location: Holland  Posted: Wed Jul 17, 2019 9:33 am    Post subject: Ok, for a solar panel it is fine Wish my roof was better pointed towards the Sun._________________Mark   Dave Joined: 05 Feb 2005
Posts: 264
Location: McMinnville, OR  Posted: Wed Jul 17, 2019 11:32 pm    Post subject: JP Thanks for your suggestions. Have you posted any pics of your astronomical clocks? I'd like to see them. D.  Duval JP  Joined: 22 Jun 2004
Posts: 840
Location: France  Posted: Thu Jul 18, 2019 9:29 am    Post subject: 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 _________________pleasure to learn, to teach, to create   Dave Joined: 05 Feb 2005
Posts: 264
Location: McMinnville, OR  Posted: Thu Jul 18, 2019 4:11 pm    Post subject: JP, WOW! Great work! Keep doing what you are doing. I will be looking closely at your code. Thanks, Dave  Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First

 Jump to: Select a forum BASCOM AVR/8051----------------BASCOM-AVRBASCOM-8051BASCOM-ARDUINOShare your working BASCOM-8051 code hereShare your working BASCOM-AVR code hereBASCOM BETA-SLA BASCOM Related----------------EASY TCP/IPAVR-DOSAR7212KokkeKat FAT-free SD card libBASCOM Project Blog Other Stuff----------------VariousPCB'sRoboticsNew WebSiteAnnouncementsAVR Archive----------------BASCOM-AVR ArchiveBASCOM-8051 ArchiveBASCOM-AVR Unsupported versionsEasy TCP/IP ArchiveBASCOM-EDB
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