Text 1198, 178 rader
Skriven 2007-07-07 20:56:06 av DAVID WILLIAMS (1:250/514)
Ärende: sun position and heliostat
==================================
Having nothing much better to do today, I wrote a QBasic program that 
calculates the position of the sun in the sky, as seen from anywhere on 
earth at ant time. It also calculates the alignment of a mirror if it 
is to reflect sunlight in a desired direction, as is done in a 
heliostat. 
  
In case anyone is interested, here it is. 
  
------------------------------------------------- 
  
' 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 
  
' Version dated 2007 Jul 07 
  
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) 
DECLARE FUNCTION Ang2 (X, Y) 
  
CONST PI = 3.1415926536# 
CONST DR = 180 / PI ' 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 - 12 - TZN) + (MIN + ET) / 4 - LG) / DR ' long diff 
CALL P2D(LD, DC / DR, sX, sY, sZ) 
RR = SQR(sY * sY + sZ * sZ) 
CL = (90 - LT) / DR ' colatitude in radians 
AA = Ang(sY, sZ) + CL 
sY = RR * COS(AA) 
sZ = RR * SIN(AA) 
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 * DR); " 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 - 180) / DR, tEL / DR, 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 * DR); " 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) * PI / 2 
  ELSE 
    A = ATN(Y / X) 
    IF X < 0 THEN A = A + PI 
  END IF 
  Ang = A 
  
END FUNCTION 
  
FUNCTION Ang2 (X, Y) 
  ' does Ang in degrees, between zero and 360 
  
  A = Ang(X, Y) * DR + 180 
  Ang2 = A - 360 * INT(A / 360) 
  
END FUNCTION 
  
SUB D2P (X, Y, Z, AZ, EL) 
  ' convert from X,Y,Z to AZ, EL 
  
  EL = Ang(SQR(X * X + Y * Y), Z) 
  AZ = Ang2(Y, X) 
  
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 PI and DR (180/PI) are already initialized 
  
IF W = 0 THEN ' first call, initialize constants 
  
 W = 2 * PI / 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)) / PI 
  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 
  
  Z = SIN(EL) 
  C = COS(EL) 
  X = C * SIN(AZ) 
  Y = C * COS(AZ) 
  
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 
  
------------------------------------------------------ 
  
P.S. I tried to send this a few minutes ago, and the system locked up. 
You may get an incomplete version, along with this one, which is 
complete. 
--- Platinum Xpress/Win/WINServer v3.0pr5
 * Origin: The Bayman BBS,Toronto, (416)698-6573 - 1:250/514 (1:250/514)