10 ' SUNCALC 20 'Program to compute sun position, sunrise, solar power collection, etc. 30 ' 40 'general comments 50 ' 60 'To use, adjust the constants in the next section to reflect your earth 70 ' position, the day and time on which you want to find the sun, and 80 ' the direction you want to face a solar collector. Default values are 90 ' noon on the solistice in DeKalb, IL, USA with sunlight falling on 100 ' a horizontal stretch of ground. 110 'Program runs in UBASIC, available in Simtel archives www.simtel.net 120 ' On a PC, type 'ubibm', then 'load "suncalc'; type 'system' when 130 ' done. Other forms of BASIC will probably work with some fiddling. 135 ' You'll need to know: #pi=3.1416; atan(x)=arctangent function. 140 'Note: 'energy collection' is a multiplier. It's the fraction of the 150 ' energy received by an equal-sized collector in constant full-frontal 160 ' sun. Assuming no air or shading, this is 1.40 kilowatt/m^2, spread 170 ' among all wavelengths. I don't know what it reduces to with air 180 ' absorption. If your collector is wavelength-sensitive, there is a 190 ' different number which depends on the absorption spectrum of air. 200 'Times in this program are relative to local solar max (here, "noon") 210 ' This moment is actually (12+xx) o'clock GMT, where xx=Longitude/15. 220 ' (Of course this is not usually an integer, e.g. might be 17:45 GMT) 230 ' Convert to local time by subtracting the appropriate number of whole 240 ' hours (ROUGHLY the integer nearest to xx) -- don't forget to adjust 250 ' for daylight savings time as well. I decided not to get involved 260 ' with the conversion; you'll see Longitude is not really used below. 270 'Notes on the model used (see accompanying text file): We ignore 280 ' the finitude of sun's distance to earth (23491 earth radii) 290 ' eccentricity of earth's orbit (dist=3.40% greater July4 than Jan4) 300 ' non-sphericity of earth (radius=0.34% larger at equator than poles) 310 'Written by dave rusin, rusin@math.niu.edu, on 4/10/94 320 ' Released to the public domain. No warranties, of course. 330 ' This program is not intended to be beautiful, fast, robust, or 340 ' amazingly accurate, but is mathematically sound. For improvements, 350 ' or for details of calculations, please contact the author. 360 ' 370 'user input section 380 ' 390 point 4:' specifies extra accuracy in UBASIC computations. 400 Slant=23.45:' angle of inclination of axis to earth orbit 410 Latitude=+42:' data for DeKalb, Illinois, USA 420 Longitude=+89:' consult an almanac or atlas to get your coordinates 430 Days=91:' days since March equinox 440 Timex=12.0:' hours since midnight at which you want to find the sun 450 Nx=0:Ex=0:Ux=1.0:'outward normal on solar collector (North,East,Up) 460 ' "up"=(0,0,1);"northwest"=(1,-1,0), etc. 470 ' 480 'now solve the problems at hand 490 ' 500 print 510 print "Earth computations, latitude";Latitude;","; 520 print Days;"days after March equinox." 530 print using(2,3),"Solar collector facing ";Nx;", ";Ex;", ";Ux 540 print "Times relative to local noon (=solar max)" 550 gosub 630:' internal calculations 560 gosub 1400:'calculations for a moment 570 gosub 1640:'calculations for a day 580 gosub 1960:'calculations for a year 590 end 600 ' 610 'subroutines follow 620 ' 630 'internal data conversions 640 ' 642 'need to normalize so nx^2+ex^2+ux^2=1. 645 Length=sqrt(nx*nx+ex*ex+ux*ux) 646 Nx=Nx/Length:Ex=Ex/Length:Ux=Ux/Length 650 Alpha=Slant*#pi/180:'convert to radians 660 Theta=Latitude*#pi/180:'note: we don't actually use Longitude 670 Psi=Days*2*#pi/365.25 680 Noon=atan(tan(Psi)*cos(Alpha)) 690 if cos(Theta)*cos(Psi)*cos(Noon)<0 then Noon=Noon+#pi 700 'glossing over cases where phi, psi, or theta = +- pi/2 710 Phix=Noon+(Timex-12)*2*#pi/24 720 return 730 ' 740 'compute apparent solar position: 750 ' (N,E,U)=cartesian coords, (Ang1,Ang2)=spherical. 760 ' Also compute energy flux here. 770 ' 780 N1=-sin(Theta)*cos(Psi)*cos(Phi) 790 N2=-cos(Alpha)*sin(Theta)*sin(Psi)*sin(Phi) 800 N3=+sin(Alpha)*cos(Theta)*sin(Psi) 810 N=N1+N2+N3 820 E=-cos(Psi)*sin(Phi)+cos(Alpha)*sin(Psi)*cos(Phi) 830 U1=cos(Theta)*cos(Psi)*cos(Phi) 840 U2=cos(Alpha)*cos(Theta)*sin(Psi)*sin(Phi) 850 U3=sin(Alpha)*sin(Theta)*sin(Psi) 860 U=U1+U2+U3 870 if E=0 then Ang1=#pi/2 880 if E<>0 then Ang3=atan(abs(N/E)) 890 if E>0 then Ang1=atan(N/E) 900 if E<0 then Ang1=atan(N/E)+#pi 910 if U=1 then Ang2=#pi/2 920 if U=-1 then Ang2=-#pi/2 930 if abs(U)<1 then Ang2=atan(U/sqrt(1-U*U)) 940 'now u=sin(ang2),e=cos(ang2)*cos(ang1),n=cos(ang2)*sin(ang1) 950 if U<0 then Energy=0:goto 980:'you get no energy with sun below horizon 960 Energy=(N*Nx+E*Ex+U*Ux)*100:if Energy<0 then Energy=0 970 ' (assuming collector is one-sided: no energy if sun is behind it). 980 return 990 ' 1000 'Technical problem: how to solve A*cos(phi)+B*sin(phi)+C=0 for phi. 1010 ' 1020 F=sqrt(A*A+B*B) 1030 if F=0 then Nroots=0:goto 1150 1040 'actually the equation with a=b=0 DOES have solutions if c=0 too 1050 'but having inf many solutions is in our setting like having none 1060 A=A/F:B=B/F:C=-C/F:' now a^2+b^2=1 1070 if B=0 then X=A*#pi/2:goto 1100 1080 X=atan(A/B):' then a=+-sin(x),b=+-cos(x) 1090 if B*cos(X)<0 then X=X+#pi:'now what we need to solve is sin(x+phi)=c 1100 if abs(C)>1 then Nroots=0:goto 1150 1110 if abs(C)=1 then Phi1=C*#pi/2:Root1=Phi1-X:Nroots=1:goto 1150 1120 Phi1=atan(C/sqrt(1-C*C)):' =arcsin(c) 1130 Root1=Phi1-X:Root2=#pi-Phi1-X:Nroots=2:goto 1150 1140 'difference is 2*arccos(-C/sqrt(A^2+B^2)). This gives day length. 1150 return 1160 ' 1170 'calculate day's energy accumulation 1180 ' 1190 'probably possible to get a closed form integral of sun*normal,but 1200 ' seems easier just to do numerical integration: 1210 Total=0 1220 Numpoints=100 1230 for I=1 to Numpoints 1240 Phi=2*#pi*I/Numpoints 1250 gosub 740 1260 Total=Total+Energy 1270 next I 1280 Total=Total/Numpoints 1290 return 1300 ' 1310 'ho-hum routine to print times nicely 1320 ' 1330 print using(3,0),int(X);":"; 1340 X=X-int(X) 1350 X=int(0.5+60*X) 1360 if X<10 then print " 0"; print using(1,0),X;:goto 1380 1370 print using(3,0),X; 1380 return 1390 ' 1400 'show data for the moment 1410 ' 1420 print 1430 print "Observations about this moment" 1440 print "******************************" 1450 print "At ";:X=Timex:gosub 1310:print ","; 1460 print " sun appears to be at coordinates: "; 1470 Phi=Phix 1480 gosub 740 1490 print using(2,3),"N=";N;" E=";E;" U=";U 1500 if U=1 then print "straight overhead":goto 1600 1510 if U=-1 then print "straight underneath":goto 1600 1520 print "Sighting is "; 1530 if E=0 then print "due ";:goto 1550 1540 print using(4,2),abs(Ang3)*180/#pi;" degrees "; 1550 if N>0 then print "north, and ";:goto 1570 1560 if N<0 then print "south, and "; 1570 print using(4,2),abs(Ang2)*180/#pi;" degrees "; 1580 if U>=0 then print "up.":goto 1600 1590 if U<0 then print "down ... but you can't see it." 1600 print "Relative energy collection rate is ";using(4,2),Energy; 1610 print " percent maximal" 1620 return 1630 ' 1640 'show data for the day 1650 ' 1660 print 1670 print "Observations about this day" 1680 print "***************************" 1690 A=cos(Theta)*cos(Psi) 1700 B=cos(Alpha)*cos(Theta)*sin(Psi) 1710 C=sin(Alpha)*sin(Theta)*sin(Psi) 1720 gosub 1000 1730 'note:similarly one can find e.g. when sun is due west, etc. 1740 if Nroots>0 then 1780 1750 print "No dawn or dusk that day -- sun is all day "; 1760 if sin(Psi)*cos(Alpha-Theta)>0 then print "shining" else print "hidden" 1770 goto 1920 1780 Dawn=(Root1-Noon)*12/#pi+12:Dawn=Dawn-24*int(Dawn/24) 1790 Dusk=(Root2-Noon)*12/#pi+12:Dusk=Dusk-24*int(Dusk/24) 1800 if Nroots=2 then 1830 1810 print "sun touches horizon just at";:X=Dawn:gosub 1310:'nroots=1 1820 print ".":goto 1880 1830 if Dawn>Dusk then X=Dawn:Dawn=Dusk:Dusk=X:Root1=Root2 1840 print "This day has";using(3,2),Dusk-Dawn;" hours of daylight: "; 1850 print "Dawn at ";:X=Dawn:gosub 1310 1860 print ", dusk at ";:X=Dusk:gosub 1310 1870 print "." 1880 Phi=Root1:gosub 740 1890 print "Sunrise at bearing ";using(3,2),Ang1*180/#pi;" degrees." 1900 Phi=Noon:gosub 740 1910 print "High sun at angle ";using(3,2),Ang2*180/#pi;" degrees" 1920 gosub 1170 1930 print "Energy collected is ";using(4,2),Total;" percent maximal" 1940 return 1950 ' 1960 'show data for the year 1970 ' 1980 print 1990 print "Observations about this place - energy collection per day" 2000 print "*********************************************************" 2010 print "Solistices:"; 2020 Days=0:gosub 630:gosub 1170:print using(4,2),Total;" % maximal" 2030 print "June equinox:"; 2040 Days=91:gosub 630:gosub 1170:print using(4,2),Total;" % maximal" 2050 print "December equinox:"; 2060 Days=273:gosub 630:gosub 1170:print using(4,2),Total;" % maximal" 2070 return