Text 1201, 282 rader
Skriven 2007-07-28 13:55:18 av DAVID WILLIAMS (1:250/514)
Ärende: heliostat code
======================
I don't seem to have posted this here, although I've done so in some 
Usenet newsgroups. 
  
I recently wrote some BASIC code that calculates the position of the 
sun in the sky, as azimuth and angle of elevation, as seen from 
anywhere on earth, at any time on any date. It also calculates the 
direction in which a mirror should be aimed if it is to reflect 
sunlight in any desired direction. A heliostat, of course, has a mirror 
that does that, and is turned so as to continue reflecting in the right 
direction as the sun moves across the sky. 
  
I'll append two versions of the program below. The first is in QBasic, 
but will also run in some other implementations. The second is in a 
very simple BASIC, and should run on almost any version of the 
language. However, it is necessarily messier than the first one. 
  
Enjoy. Let me know if you have any questions. 
  
                      dow 
  
------------------------------------------------ 
  
' SunAlign.BAS 
' Calculates position of sun in sky, as azimuth (true 
' compass bearing) 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 
' P.O. Box 48512 
' Toronto, Canada 
' M8W 4Y6 
  
' Initially dated 2007 Jul 07 
' This version 2007 Jul 22 
  
DECLARE FUNCTION ET.Dec (D, F%) 
DECLARE FUNCTION ROff$ (X) 
DECLARE SUB D2P (X, Y, Z, AZ, EL) 
DECLARE SUB P2D (AZ, EL, X, Y, Z) 
DECLARE FUNCTION Ang (X, Y) 
  
CONST PY = 3.1415926536# ' "PI" not assignable in some BASICs 
CONST DR = 180 / PY ' degree / radian factor 
  
CLS 
PRINT "Use negative numbers for directions opposite to those shown." 
INPUT "Observer's latitude (degrees North)"; LT 
INPUT "Observer's longitude (degrees West)"; LG 
INPUT "Date (M#,D#)"; Mth%, Day% 
INPUT "Time (HH,MM) (24-hr format)"; HR, MIN 
INPUT "Time Zone (+/- hours from GMT/UT)"; TZN 
  
D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365 
ET = ET.Dec(D, 1) ' Equation of Time 
DC = ET.Dec(D, 0) ' Declination 
  
LD = 15 * (HR - TZN) + (MIN + ET) / 4 - LG  'longitude difference 
CALL P2D(LD, DC, sX, sY, sZ) 
CALL D2P(sY, sZ, sX, qAZ, qEL) 
CALL P2D(qAZ + LT - 90, qEL, sY, sZ, sX) 
CALL D2P(sX, sY, sZ, sAZ, sEL) 
  
PRINT 
IF sEL < 0 THEN PRINT "Sun Below Horizon": END 
PRINT "Sun's azimuth: "; ROff$(sAZ); " degrees" 
PRINT "Sun's elevation: "; ROff$(sEL); " degrees" 
  
PRINT 
PRINT "Calculate heliostat mirror alignment? (y/n) "; 
DO 
 K$ = UCASE$(INKEY$) 
LOOP UNTIL K$ = "Y" OR K$ = "N" 
PRINT K$ 
  
IF K$ = "Y" THEN 
   PRINT 
   INPUT "Azimuth of target direction (degrees)"; tAZ 
   INPUT "Elevation of target direction (degrees)"; tEL 
   CALL P2D(tAZ, tEL, tX, tY, tZ) 
   CALL D2P(sX + tX, sY + tY, sZ + tZ, mAZ, mEL) 
  
   PRINT 
   PRINT "Mirror aim direction (perpendicular to surface):" 
   PRINT "Azimuth: "; ROff$(mAZ); " degrees" 
   PRINT "Elevation: "; ROff$(mEL); " degrees" 
END IF 
  
END 
  
FUNCTION Ang (X, Y) 
 ' calculates angle between positive X axis and vector to (X,Y) 
  
  IF X = 0 THEN 
    A = SGN(Y) * PY / 2 
  ELSE 
    A = ATN(Y / X) 
    IF X < 0 THEN A = A + PY 
  END IF 
  Ang = A 
  
END FUNCTION 
  
SUB D2P (X, Y, Z, AZ, EL) 
  ' convert from X,Y,Z to AZ, EL 
  ' Outputs in degrees 
  
  EL = Ang(SQR(X * X + Y * Y), Z) * DR 
  A = Ang(Y, X) * DR 
  IF A < 180 THEN AZ = A + 180 ELSE AZ = A - 180 
  
END SUB 
  
FUNCTION ET.Dec (D, F%) STATIC 
  ' Calculates equation of time, in minutes, or solar declination, 
  ' in degrees, on day number D of year. (D = 0 on January 1.) 
  ' F% selects function: True (non-zero) for Equation of Time, 
  ' False (zero) for Declination. 
  ' STATIC means variables are preserved between calls of function 
  ' This version assumes PY (PI) and DR (180/PI) are initialized 
  
IF W = 0 THEN ' first call, initialize constants 
  
 W = 2 * PY / 365 ' earth's mean orbital angular speed in radians/day 
 C = -23.45 / DR ' reverse angle of earth's axial tilt in radians 
 ST = SIN(C) ' sine of reverse tilt 
 CT = COS(C) ' cosine of reverse tilt 
 E2 = 2 * .0167 ' twice earth's orbital eccentricity 
 SP = 12 * W ' 12 days from December solstice to perihelion 
 D1 = -1 ' holds last D. Saves time if D repeated for both functions 
  
END IF 
  
IF D <> D1 THEN ' new value of D 
  A = W * (D + 10) ' Solstice 10 days before Jan 1 
  B = A + E2 * SIN(A - SP) 
  D1 = D 
END IF 
  
IF F% THEN  ' equation of time calculation 
  C = (A - ATN(TAN(B) / CT)) / PY 
  ET.Dec = 720 * (C - INT(C + .5)) 
  ' in 720 minutes, earth rotates PI radians relative to sun 
  
ELSE ' declination calculation 
  C = ST * COS(B) 
  ET.Dec = ATN(C / SQR(1 - C * C)) * DR 
  ' arcsine of C in degrees. ASN not directly available in QBasic 
  
END IF 
  
END FUNCTION 
  
SUB P2D (AZ, EL, X, Y, Z) 
  ' convert from AZ,EL to X,Y,Z 
  ' Inputs in degrees 
  
  E = EL / DR 
  A = AZ / DR 
  Z = SIN(E) 
  C = -COS(E) 
  X = C * SIN(A) 
  Y = C * COS(A) 
  
END SUB 
  
FUNCTION ROff$ (X) 
 ' neatly rounds number to one place of decimals 
  
S$ = LTRIM$(STR$(INT(10 * ABS(X) + .5))) 
IF S$ = "3600" THEN S$ = "0" 
IF LEN(S$) = 1 THEN S$ = "0" + S$ 
IF X < 0 THEN IF VAL(S$) THEN S$ = "-" + S$ 
ROff$ = LEFT$(S$, LEN(S$) - 1) + "." + RIGHT$(S$, 1) 
  
END FUNCTION 
  
---------------------------------------------------- 
  
100 REM SunAlign.BAS 
  
110 REM Calculates position of sun in sky, as azimuth (true 
120 REM compass bearing) and angle of elevation, as seen from 
130 REM any place on earth, on any date and any time. 
140 REM Also calculates alignment of a heliostat mirror. 
  
150 REM David Williams 
160 REM P.O. Box 48512 
170 REM Toronto, Canada 
180 REM M8W 4Y6 
  
190 REM Initially dated 2007 Jul 07 
200 REM This version 2007 Jul 22 
  
210 REM Note: For brevity, no error checks on user inputs 
  
220 CLS 
230 PRINT "Use negative numbers for opposite directions." 
240 INPUT "Observer's latitude (degrees North)"; LT 
250 INPUT "Observer's longitude (degrees West)"; LG 
260 INPUT "Date (M#,D#)"; Mth, Day 
270 INPUT "Time (HH,MM) (24-hr format)"; HR, MIN 
280 INPUT "Time Zone (+/- hours from GMT/UT)"; TZN 
  
290 PY = 4 * ATN(1): REM "PI" not assignable in some BASICs 
300 DR = 180 / PY: REM degree/radian factor 
310 W = 2 * PY / 365: REM earth's mean orbital speed in radians/day 
320 C = -23.45 / DR: REM reverse angle of axial tilt in radians 
330 ST = SIN(C): REM sine of reverse tilt 
340 CT = COS(C): REM cosine of reverse tilt 
350 E2 = 2 * .0167: REM twice earth's orbital eccentricity 
360 SP = 12 * W: REM 12 days from December solstice to perihelion 
  
370 D = INT(30.6 * ((Mth + 9) MOD 12) + 58.5 + Day) MOD 365 
380 A = W * (D + 10): REM Solstice 10 days before Jan 1 
390 B = A + E2 * SIN(A - SP) 
  
400 C = (A - ATN(TAN(B) / CT)) / PY 
410 ET = 720 * (C - INT(C + .5)): REM equation of time 
420 REM in 720 minutes, earth rotates PI radians relative to sun 
  
430 C = ST * COS(B) 
440 EL = ATN(C / SQR(1 - C * C)) * DR: REM solar declination 
  
450 AZ = 15 * (HR - TZN) + (MIN + ET) / 4 - LG: REM longitude diff 
460 GOSUB 800 
470 R = SQR(Y * Y + Z * Z) 
480 AX = Y: AY = Z: GOSUB 710 
490 A = AA + (90 - LT) / DR 
500 Y = R * COS(A) 
510 Z = R * SIN(A) 
520 GOSUB 740 
  
530 PRINT : REM AZ & EL are now sun's azimuth & elevation in degrees 
540 IF EL < 0 THEN PRINT "Sun Below Horizon": END 
550 R = AZ: GOSUB 870: PRINT "Sun's azimuth: "; R; " degrees" 
560 R = EL: GOSUB 870: PRINT "Sun's elevation: "; R; " degrees" 
  
570 PRINT 
580 INPUT "Calculate heliostat mirror alignment (y/n)"; K$ 
590 IF K$ = "N" OR K$ = "n" THEN END 
  
600 SX = X: SY = Y: SZ = Z 
  
610 PRINT 
620 INPUT "Azimuth of target direction (degrees)"; AZ 
630 INPUT "Elevation of target direction (degrees)"; EL 
640 GOSUB 800 
650 X = X + SX: Y = Y + SY: Z = Z + SZ: GOSUB 740 
  
660 PRINT : REM AZ & EL are now aim azimuth & elevation in degrees 
670 PRINT "Mirror aim direction (perpendicular to surface):" 
680 R = AZ: GOSUB 870: PRINT "Azimuth: "; R; " degrees" 
690 R = EL: GOSUB 870: PRINT "Elevation: "; R; " degrees" 
  
700 END 
  
710 IF AX = 0 THEN AA = SGN(AY) * PY / 2: RETURN 
720 AA = ATN(AY / AX): IF AX < 0 THEN AA = AA + PY 
730 RETURN 
  
740 AX = SQR(X * X + Y * Y): AY = Z: GOSUB 710 
750 EL = AA * DR 
760 AX = Y: AY = X: GOSUB 710 
770 AZ = AA * DR 
780 IF AZ < 180 THEN AZ = AZ + 180 ELSE AZ = AZ - 180 
790 RETURN 
  
800 E = EL / DR 
810 A = AZ / DR 
820 Z = SIN(E) 
830 C = 0 - COS(E): REM Won't work without "0" in Liberty Basic 
840 X = C * SIN(A) 
850 Y = C * COS(A) 
860 RETURN 
  
870 R = INT(10 * R + .5) / 10: REM roundoff to 1 decimal place 
880 RETURN 
  
--------------------------------------------------- 
--- Platinum Xpress/Win/WINServer v3.0pr5
 * Origin: The Bayman BBS,Toronto, (416)698-6573 - 1:250/514 (1:250/514)