Text 1087, 409 rader
Skriven 2007-06-12 09:58:38 av DAVID WILLIAMS (1:250/514)
     Kommentar till en text av Miles Maxted
Ärende: Re: sun position program
================================
-> You put up some programmes for calculating gnonomic angles last  
-> October,  following some discussion about a disabled sundial here  
  
In case anyone here is interested, and can still run QBasic, here's 
another program to do with the sun's position. 
  
                      dow 
  
------------------------------------------------- 
  
' ETIMSDEC.BAS 
' David Williams, 2003 
  
' P.O. Box 48512 
' Toronto, Ontario 
' Canada. M8W 4Y6 
  
' This version dated 2007 May 25 
  
DECLARE FUNCTION ET.Dec (D, F%) 
DECLARE FUNCTION ROff$ (X) 
DECLARE SUB Waitkey (X%) 
  
CL0 = 7: CL1 = 14: CL2 = 11   ' colours used 
  
SCREEN 12 
  
COLOR CL2 
CLS 
  
PRINT "ETIMSDEC.BAS" 
PRINT 
PRINT "Shows Equation of Time and Solar Declination calculations" 
PRINT "(performed in function ET.Dec) and compares their results" 
PRINT "graphically with published values, showing close agreement." 
PRINT 
PRINT "These calculations can be used in the programs of computer-" 
PRINT "controlled solar-energy equipment, such as sun trackers and" 
PRINT "heliostats." 
PRINT 
PRINT "Equation of Time is difference between Solar Time and Mean" 
PRINT "Time. Sundials show solar time. Clocks show mean time. The" 
PRINT "Equation of Time is related to solar longitude, relative to" 
PRINT "its mean value (the value it would have if the sun's apparent" 
PRINT "motion were at uniform speed). One degree westward of solar" 
PRINT "longitude corresponds to four minutes in positive direction in" 
PRINT "the Equation of Time." 
PRINT 
PRINT "Solar Declination is latitude of sun in celestial coordinates." 
PRINT 
  
' initialize array used in graphing routine 
DIM ML(1 TO 13) 
FOR X = 1 TO 13 
  READ ML(X) 
NEXT 
DATA 33,29,33,31,33,31,33,33,31,33,31,33,0 
' Month lengths adjusted for numbers of pixels on screen 
  
Menu: 
  
COLOR CL0 
  
PRINT "1. Find Equation of Time and Solar Declination on given date" 
PRINT "2. Draw Equation of Time Graph" 
PRINT "3. Draw Solar Declination Graph" 
PRINT "4. Draw Solar Analemma" 
PRINT "5. Quit program" 
PRINT 
PRINT "Which (1 - 5)? "; 
WHILE INKEY$ <> "": WEND 
DO 
  K$ = INKEY$ 
LOOP UNTIL K$ >= "1" AND K$ <= "5" 
PRINT K$ 
  
SELECT CASE VAL(K$) 
 CASE 1: GOTO Calc 
 CASE 2, 3: GOTO Graph 
 CASE 4: GOTO Analemma 
 CASE 5: GOTO WayOut 
END SELECT 
  
Calc: 
  
PRINT 
COLOR CL1 
  
DO 
  INPUT "Date (M#,D#)"; Mth%, Day% 
LOOP UNTIL Mth% > 0 AND Mth% < 13 AND Day% > 0 AND Day% < 32 
  
D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365 
' day of 365-day year. Jan 1 = 0. Dec 31 = 364 
  
PRINT 
PRINT "Day number"; D + 1; "of year" ' conventional range 1 - 365 
  
ET = ET.Dec(D, 1) ' Equation of Time 
DC = ET.Dec(D, 0) ' Declination 
  
PRINT 
PRINT "Equation of Time: "; ROff$(ET); " minutes" 
PRINT "Solar Declination: "; ROff$(DC); " degrees" 
PRINT 
  
GOTO Menu 
  
Graph: 
  
F% = (K$ = "2") ' true if Equation of Time selected 
  
' Set up graph 
CLS 
  
IF F% THEN 
  PRINT TAB(30); "EQUATION OF TIME" 
  LOCATE 3, 15 
  PRINT "Graph shows difference in minutes between clock" 
  PRINT TAB(15); "and sundial time. Positive difference means" 
  PRINT TAB(15); "sundial is ahead of clock, and vice versa." 
  LOCATE 24, 15 
  COLOR CL1 
  PRINT "Graph shows results of calculation." 
  COLOR CL2 
  PRINT TAB(15); "Circles show published values of Equation of Time." 
  RESTORE ETData 
ELSE 
  PRINT TAB(25); "Declination of Sun, in Degrees" 
  LOCATE 3, 10 
  COLOR CL1 
  PRINT "Graph shows calculated function. "; 
  COLOR CL2 
  PRINT "Circles show published values."; 
  RESTORE DeclnData 
END IF 
  
READ ML%, VT%, VB%, LL%, UL% 
  
COLOR CL0 
LOCATE 16, 67 
PRINT "-="; 
  
FOR T = LL% TO UL% STEP 5 
  LINE (137, 247 - 6.4 * T)-(530, 247 - 6.4 * T) 
  LOCATE 16 - T / 2.5, 14 
  PRINT RIGHT$(" " + STR$(T), 3); 
  IF T = 0 THEN PRINT " ="; 
NEXT 
LOCATE ML%, 20 
PRINT "Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec"; 
  
' Do graphing 
T = 0 
D = 0 
MC = 0 
FOR M = 0 TO 384 
  M1 = 147 + M 
  IF M = T THEN 
    ' Draw vertical line 
    COLOR CL0 
    LINE (M1, VT%)-(M1, VB%) 
    MC = MC + 1 
    T = T + ML(MC) 
    MS = 0 
  ELSE 
    ' Plot point(s) 
    COLOR CL1 
    G = ET.Dec(D, F%) ' Function used 
    PSET (M1, 247 - 6.4 * G) 
    IF MS = 0 OR MS = 10 OR MS = 20 THEN 
      ' Plot circle representing published value 
      COLOR CL2 
      READ E    ' published values 
      E = FIX(E) + (E - FIX(E)) * 5 / 3 
      CIRCLE (M1, 247 - 6.4 * E), 2 
    END IF 
    MS = MS + 1 
    D = D + 1 
  END IF 
  IF ML(MC) = 33 AND M = T - 18 THEN M = M + 1 
NEXT 
  
Waitkey 25 
  
GOTO Menu 
  
ETData: 
  
DATA 7, 97, 343, -15, 20: ' plotting parameters for E. of T. graph 
  
DATA -3.12,-7.38,-11.08,-13.33,-14.19,-13.49,-12.34,-10.18,-7.28 
DATA -4.08,-1.16,1,2.51,3.4,3.34,2.25,.39,-1.28 
DATA -3.33,-5.16,-6.15,-6.16,-5.14,-3.16,-.12,3.08,6.4 
DATA 10.05,13.02,15.12,16.2,16,14.16,11.11,7.02,2.13 
' Above data show published values of Equation of Time on 
' 1st, 11th and 21st days of months. Decimal point separates 
' minutes, to left, from seconds, to right, e.g. 3.16 means 
' 3 minutes and 16 seconds. Seconds are multiplied by 5/3 to 
' get fractional minutes after READ instruction. 
  
DeclnData: 
  
DATA 5, 65, 407, -25, 25: ' plotting parameters for Decl'n graph 
  
DATA -23.04,-21.56,-20.05,-17.2,-14.18,-10.52,-7.49,-3.57,0 
DATA 4.18,8.07,11.39,14.54,17.43,20.04,21.58,23.02,23.26 
DATA 23.09,22.11,20.36,18.1,15.27,12.19,8.3,4.48,.57 
DATA -2.57,-6.48,-10.29,-14.14,-17.15,-19.47,-21.43,-22.57,-23.26 
' Above data show published values (in degrees and minutes) of solar 
' declination on 1st, 11th and 21st days of month. See note after 
' Equation of Time data. 
  
' Published values are from the book: 
' Sundials, Their Theory and Construction 
' by: Albert E. Waugh 
' Dover Publications, New York, 1973 
' ISBN 0-486-22947-5 
  
Analemma: 
  
CLS 
COLOR CL0 
  
RESTORE Anadata 
  
PRINT "SOLAR ANALEMMA" 
PRINT 
PRINT 
PRINT "Shows apparent path of sun," 
PRINT "through the year, relative" 
PRINT "to its mean position" 
PRINT "(shown by cross)." 
  
LINE (475, 52)-(475, 90)  ' compass lines 
LINE (456, 71)-(494, 71) 
  
LOCATE 3, 60 
PRINT "N"; 
LOCATE 7, 60 
PRINT "S"; 
LOCATE 5, 56 
PRINT "E"; 
LOCATE 5, 64 
PRINT "W"; 
  
LOCATE 9, 1 
PRINT TAB(53); "Looking at sky," 
PRINT TAB(53); "if North is upward," 
PRINT TAB(53); "East is to left." 
  
LOCATE 24, 1 
PRINT "For dimensions, see graphs" 
PRINT "of Solar Declination and" 
PRINT "Equation of Time." 
PRINT "+4 minutes in Equation of Time" 
PRINT "equal 1 degree Westward." 
  
LOCATE 15, 1 
PRINT TAB(50); "Celestial" 
PRINT TAB(50); "Equator" 
  
X0 = 320 
Y0 = 240 
  
LINE (X0, Y0 + 5)-(X0, Y0 - 5) ' cross at (0,0) 
LINE (X0 + 5, Y0)-(X0 - 5, Y0) 
  
FOR S = -1 TO 1 STEP 2 
  LINE (X0, Y0 + S * 210)-(X0, Y0 + S * 225) ' lines marking N-S 
  LINE (X0 + S * 30, Y0)-(X0 + S * 55, Y0)   ' and E-W 
NEXT 
  
LOCATE 10, 1 
PRINT "Month colours:"; 
  
M = 0 
DM = 0 
LM = 0 
  
' draw analemma 
FOR DY = 0 TO 364 
  IF DM = LM THEN 
    M = M + 1 
    C = M XOR 14 * (M AND 1) 
    IF C = 8 THEN C = 14 
    COLOR C 
    READ LM, MN$ 
    PRINT TAB(17); MN$ 
    DM = 0 
  END IF 
  DM = DM + 1 
  FOR D = DY TO DY + .99 STEP .2 
    E = ET.Dec(D, 1) 
    L = ET.Dec(D, 0) 
    PSET (X0 + 2.5 * E, Y0 - 10 * L) 
  NEXT 
NEXT 
  
COLOR CL0 
Waitkey 48 
  
GOTO Menu 
  
Anadata: 
  
DATA 31,Jan,28,Feb,31,Mar,30,Apr,31,May,30,Jun 
DATA 31,Jul,31,Aug,30,Sep,31,Oct,30,Nov,31,Dec 
' month lengths and names 
  
WayOut: 
  
COLOR CL2 
CLS 
PRINT "These calculations of Equation of Time and Solar Declination" 
PRINT "are simplified and approximate. However, they are quite good." 
PRINT "All differences between calculated and published values of" 
PRINT "the solar declination are small compared with the angular" 
PRINT "size of the sun in the sky. Similarly, differences between" 
PRINT "calculated and published values of the equation of time are" 
PRINT "small compared with the time the sun takes to traverse its" 
PRINT "own diameter as it moves across the sky. The non-zero size" 
PRINT "of the sun, rather than any inaccuracies of calculation," 
PRINT "is the limiting factor on how accurately the sun can be" 
PRINT "tracked in solar energy applications, using this routine." 
PRINT 
PRINT "The Equation of Time is treated here as the correction to be" 
PRINT "subtracted from a sundial reading to get local mean (clock)" 
PRINT "time. It is therefore positive when sundial is ahead of clock." 
PRINT "This is the usual sign convention, but the opposite usage is" 
PRINT "sometimes found. Take care when comparing values of the" 
PRINT "Equation of Time from different sources." 
PRINT 
PRINT "Note that atmospheric refraction of light affects the apparent" 
PRINT "position of the sun in the sky when it is close to the" 
PRINT "horizon. Sunlight is then too weak to be used for most" 
PRINT "solar-energy applications, so refraction is not usually" 
PRINT "included in calculations related to solar energy. This program" 
PRINT "does not take account of refraction." 
  
COLOR 15 ' default white 
  
END 
  
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 
  
IF PI = 0 THEN ' first call, initialize constants 
  
 PI = 4 * ATN(1) ' ATN = arctangent 
 W = PI / 182.5'earth's mean orbital angular speed in radians per day 
 DR = 180 / PI ' degree/radian factor 
 C = -23.45 / DR ' reverse angle of earth's axial tilt 
 ST = SIN(C) ' sine of reverse tilt 
 CT = COS(C) ' cosine of reverse tilt 
 E2 = 2 * .0167 ' twice earth's orbital eccentricity 
 ZI = 12 * W ' 12 days from 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 - ZI) 
  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 
  
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 
  
FUNCTION ROff$ (X) 
 ' neatly rounds number to one place of decimals 
  
S$ = LTRIM$(STR$(INT(10 * ABS(X) + .5))) 
IF LEN(S$) = 1 THEN S$ = "0" + S$ 
IF VAL(S$) THEN S$ = MID$("-x+", SGN(X) + 2, 1) + S$ 
ROff$ = LEFT$(S$, LEN(S$) - 1) + "." + RIGHT$(S$, 1) 
  
END FUNCTION 
  
SUB Waitkey (X%) 
 ' prints prompt at position X% on bottom row, then waits 
  
LOCATE 30, X% 
PRINT "*** Press any key to continue ***"; 
WHILE INKEY$ <> "": WEND 
WHILE INKEY$ = "": WEND 
CLS 
  
END SUB 
  
------------------------------------------------- 
--- Platinum Xpress/Win/WINServer v3.0pr5
 * Origin: The Bayman BBS,Toronto, (416)698-6573 - 1:250/514 (1:250/514)