%%HP: T(3)A(R)F(.); @ SUNCLK, by James Elliott @ Directory: SunClk @ Checksum: # 40750d Bytes: 8805.5 DIR Now @ Put the current date and time into TM and DT, convert to Julian @ date/time, and call SunClk, the driver program. \<< DATE TIME 'TM' STO 'DT' STO GTime 1 'STAT' STO SunClk \>> Then @ Create a temporary menu that allows specification of an arbitrary @ date and time, and display the current settings. Trap errors in case @ an invalid date or time is entered, so garbage isn't left on the @ stack. \<< { { ">DATE" @ Menu key for setting the date to view. \<< \-> i \<< DT \-> d \<< i IF DUP DUP IP == @ If only the month was entered, default the rest from the previous @ setting. THEN d FP + ELSE IF DUP FP 100 * DUP IP == @ Similarly, if the year is missing, provide the default. THEN d 100 * FP 100 / + END END 'DT' STO @ Update the date to reflect their changes, fixing things if there is @ an error. IFERR ShwTm THEN ERRM 3 DISP 2 WAIT DROP DROP d 'DT' STO ShwTm END \>> \>> \>> } { ">TIME" @ Menu key for setting the time to view. Very similar to above. \<< \-> i \<< TM \-> t \<< i 'TM' STO IFERR ShwTm THEN ERRM 3 DISP 2 WAIT DROP DROP t 'TM' STO ShwTm END \>> \>> \>> } { "A/PM" @ Toggle AM vs. PM \<< TM 12 - DUP IF 0 < THEN 24 + END 'TM' STO ShwTm \>> } { } { } { "Go" @ Use the specified time; convert it to Julian time and call SunClk, @ the driver program, restoring the normal menu. \<< 0 MENU GTime 0 'STAT' STO SunClk \>> } } TMENU @ Okay, the temporary menu is built; now show what the default time @ and date are (these are the ones used for the last display.) ShwTm \>> @ Provide brief instructions. Blank the menu while waiting for @ keypresses between screens. Help \<< CLLCD { } TMENU "Press NOW for current" 1 DISP "daylight map, or THEN" 2 DISP "to pick a day & time." 3 DISP "Set TZ to your time" 5 DISP "zone (eg. CDT = -5)." 6 DISP -1 WAIT DROP CLLCD "To interrupt, press" 1 DISP "any key other than" 2 DISP "ON/ATTN, and it will" 3 DISP "stop in a few seconds" 4 DISP "and clean up." 5 DISP 0 WAIT DROP CLLCD "'NOW' updates the map" 1 DISP "continually until you" 2 DISP "interrupt it, 'THEN'" 3 DISP "draws a map and ends." 4 DISP 0 MENU 3 FREEZE \>> @ Provide brief credits. Blank the menu while waiting for @ keypresses between screens. About \<< CLLCD { } TMENU "HP 48SX implementation" 1 DISP "by James J. Elliott" 2 DISP "" 3 DISP "Public domain; share" 5 DISP "and enjoy!" 6 DISP -1 WAIT DROP CLLCD "Based on a SunView" 1 DISP "program by:" 2 DISP "John Walker" 4 DISP "" 5 DISP 0 WAIT DROP CLLCD "Algorithms found in" 1 DISP "'Astronomical Formulae" 2 DISP "for Calculators' by" 3 DISP "Jean Meeus, Third" 4 DISP "Edition, 1985." 5 DISP 3 FREEZE 0 MENU \>> @ This variable contains your current time zone. My default is Central @ Daylight savings Time, or 5 hours earlier than UTC. TZ -5 @ This variable determines how detailed (and therefore slow) the @ computation of the daylight bands' width is. In projecting the @ sunlight region over the daylight hemisphere, this many different @ bands are computed, and linear interpolation is used between them up @ to the resolution of the graphics display. Note that it does not @ make sense to set this higher than 62, since only 62 values are @ stored! The lower the number, the chunkier the display, but the @ faster arbitrary dates can be displayed. It doesn't really affect @ the speed of ongoing updates, since this computation is only done @ when the sun's declination drifts beyond a threshold from its @ current value. Res 40 @ This is the threshold referred to above. Thres .5 @ This is just a band of black used in drawing the sunlight (with @ GXOR). It's much faster than using LINE. STRIP GROB 131 1 FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF0 @ Here's the world map. World GROB 126 64 FFFFFFF30000008FFFFFFEFFFFFFFFF3FFFFFF1000CFF18F10F00CF30FFFFFF3FFFF7000008FF7EF18FF38F30CF1EFF3FFFF70000E3FF1EFFFFF002840018FF3E30E1000007CF07EF0CF010FF5000892C8F00000004C78FF340000CF9FFF910208FF744000CC038F1008F30FFFFFFF72389FF0C000E0E3C7800EFFFFFFFFF132F000E18971E1FFF0000FFFFFFFF70083F3C3C72324CFFFF8004FFFFF9FF1F0F370FF8F7E0D1FFF700CF3FFCF1EF1E8F30EFF0F74270FFF70FFF9FFDF3FF7ECF3FFFF3E17020EFFF9FF7D7EFFFFF36EF3FFFF7F73014EFFFB0F03270FFFF10FF3FFFF7F7F00EFFF7006122FCFFF74CFF3FFFF7FFF88FFFE730024CFFFF703FFF3FFFF7EFFFEFFF8700184EFFFF701FFF3FFFFF8FF3EFFFF7883BFFFFFF720FFF3FFFFF3F78FFFFF4F3083FFFFF76EFFF3FFFFF3E09FFFFF1FFF93CFFFFFEFFFF3FFFFF3CE3FFFFF9FFF170CFFF7EFFFF3F7FFFF8E0FFFFFCFFF0F81F870EFFFF3F7EFFFB000FFFFCFFF3E93727EEFFFF3FFFFFF3040EFFFDFFF7EC73766EFFFF3FFFFFFF1CFFFFFCF3F70E7A346CFFFF3FFFFFFEFC0CFFF9F3FF1FF8B06CFFBD3FFFFFFEF008FFFBBFFF7FF8B078FFBD3FFFFFFFF3F3EFF30CFB3FF93019FFFF3FFFFFFFF7FFCFFFFCF89FFF7001FFFF3FFFFFFF33FF0FFFFCE8CFFF7820AFF73FFFFFFFB3FF4CFFFDE4EFFFF00002F73FFFFFFFFBFFF9FFF9F6FFFFF10E90EF3FFFBFFFF3FFFBF9FB76FFFFF708388F3FFFBFFFF7EFF9FFFBF64FFFFF7448FF3FFFFFFFFFAFFDFFFBFC0FFFFFF15EFE3FFFDFFFFF8FFDFFC3300FFFFFFC1CFC3FFFDFFFFFBFFCFFF3BB1DFFFF3EF97C3FFFFFFFFFBF3EFFF3FB9DFFFF9FF3FE3FFFFFFFFFBF9FFFF7F99FFFFF9FF7FF3FFFFFFFFF9F9FFFF7ECFFFFFF1E17FF3FFFFFFFFFDFCFFFFF6EFFFFFF3107FF3FFFFFFFFFD3EFFFFF0FFFFFFF3C93F33FFFFFFFFFD3FFFFFFFFFFFFFFFF18F32FFFFFFFFFD8FFFFFFFFFFFFFFFFFCF32FFFFFFFFFCEFFFFFFFFFFFFFFFFFCF03FFFFFFFFF4EFFFFFFFFFFFFFFFFFFF83FFFFFFFFF4FFFFFFFFFFFFFFFFFFFFE3FFFFFFFFF09FFFFFFFFFFFFFFFFFFFF3FFFFFFFFF0EF9FFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFF9FFFFFFFFFFFFFFFFFFFF3FFFFFFFFF7CFFFFFFFFF0E300000EFF3FFFFFFFF70CFFFFF7F00E18FFFFF08F3FFFFFF7040DFFF7000EFFFFFFFFFF3C3FFF1000C08CFF30FFFFFFFFFFFFFF3C3F70C9FFFF70E38FFFFFFFFFFFFFFF3E3F90FFFFFFFF1CFFFFFFFFFFFFFFFF1E3F7EFFFFFFFFFFFFFFFFFFFFFFFFFF1028F9FFFFFFFFFFFFFFFFFFFFFFFFFF7C3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3 @ (End of amendments) @ This reverses the pixels of a region of the line being built up, @ starting at the pixel in level 2, and ending at the pixel in level @ 1. It's called by XSpan, below. XLine \<< \-> s f \<< STRIP f 1 + s - R\->B # 0d 2 \->LIST STRIP GXOR s 2 + R\->B # 0d 2 \->LIST SWAP GXOR \>> \>> @ Reverse a region of the line, centered on the pixel in level 1, and @ twice the length specified in level 2. Handles wrapping off the left @ edge if needed. XSpan \<< \-> l n \<< l 126 MOD 'l' STO IF l n + 126 > THEN l 125 XLine 0 l n + 127 - XLine ELSE l l n + 1 - XLine END \>> \>> @ Draw the time at the bottom of the map. XTime \<< \-> s \<< s 1 \->GROB DUP SIZE DROP \-> o w \<< PICT 131 w - 2 / # 59d 2 \->LIST o GXOR \>> \>> \>> @ Update the sunlit region display if necessary, and the time display @ if appropriate. Takes the pixel location of noon in level 1. Draw \<< \-> n \<< 0 255 \-> w o \<< @ Is this an update, or the first display? If an update, set flag 6, @ so that the old lines will be taken into account. IF 'OTAB' VTYPE 0 \>= THEN 6 SF @ Has the sun moved, as far as the resolution of the map is concerned? @ If so, set flag 7. IF n ONOON \=/ THEN 7 SF ELSE 7 CF END ELSE 6 CF 7 CF END 1 64 @ Loop over the whole map, one horizontal band at a time. FOR i IF 6 FS? @ If it's an update, determine what the width was last time. THEN OTAB i DUP SUB NUM 'o' STO END @ Figure out the width of sunlight at this latitude, by looking it up @ in the width table built by PrjIll. WTAB i DUP SUB NUM 'w' STO @ If the sun has moved or the width has changed, we need to update @ this band. IF o w \=/ 7 FS? OR THEN @ Create a blank line in which the regions which need to change will @ be marked. # 131d # 1d BLANK @ If there was an old sunlit band, we'll erase it. IF o 255 < THEN ONOON o - 126 + 126 MOD o 2 * XSpan END @ If there is a new sunlit band, draw it. IF w 255 < THEN n w - 126 + 126 MOD w 2 * XSpan END @ Take the net result of the above two operations and apply it to the @ actual map in one fast GXOR, so there's no flicker. PICT SWAP # 0d i 1 - R\->B 2 \->LIST SWAP GXOR END @ If the user wants to stop, record that fact and abort the loop. IF KEY THEN DROP -1 'STAT' STO 63 'i' STO END NEXT \>> \>> @ If the clock display is enabled, draw the time at the bottom of the @ map (erasing the old time first, since GXOR is being used.) IF -40 FS? THEN OSTR XTime DT TM TSTR DUP 'OSTR' STO XTime END \>> @ Main driver program: Compute the relevant astronomical variables and @ then draw the map. SunClk @ First make sure the right mathematical modes are in effect. \<< RCLF -19 CF -15 CF -16 CF -17 SF -18 CF CLLCD @ Display an explanatory screen during the initial computation. " SunClock" 1 DISP " Figuring widths of" 3 DISP " daylight bands." 4 DISP "[Preliminaries]" 7 DISP PICT PURGE @ Put the completely dark map into PICT. PICT { # 2d # 0d } World REPL "" @ Note that no time string has yet been displayed. 'OSTR' STO @ This is the animation loop; if "NOW" was pressed, it will be @ executed multiple times. DO JTime 0 0 0 0 \-> jt gt ra dec xl \<< jt DUP GMST 'gt' STO 0 @ Figure out the location of the sun at the specified instant. SunPos DROP DROP 'dec' STO 'ra' STO @ We only care about the angular information. IF dec @ If the declination has changed more than the threshold value, or if @ it has changed sign, recompute the table of daylight widths. OldDec - ABS Thres \>= dec SIGN OldDec SIGN \=/ OR THEN dec PrjIll END @ Unless the user has already hit a key (requesting termination), draw @ the resulting map. IF STAT 0 \>= THEN { # 0d # 0d } PVIEW @ Figure out the longitude of noon. ra gt 15 * - 180 + FixAng 126 * 360 / @ Draw the current map. FLOOR DUP Draw @ Keep track of the last location of noon. 'ONOON' STO END \>> @ If the user has hit a key, end the animation loop if it was active. IF KEY THEN DROP -1 'STAT' STO END @ If we're animating, update the date and time and record the old @ widths table so they can be erased properly on the next update. IF STAT 0 > THEN DATE TIME 'TM' STO 'DT' STO GTime WTAB 'OTAB' STO END @ If we're animating, loop back and recompute/redraw. UNTIL STAT 0 \<= @ Get rid of extraneous variables. END STAT { OTAB ONOON OSTR PPAR STAT } PURGE @ If the user hasn't hit a key but we're ending because "THEN" was the @ calling program (one-shot date display), freeze the screen so the @ user can enjoy their map. IF 0 \>= THEN 7 FREEZE END STOF \>> @ Figure out what widths of sunlight fall on the latitude bands of the @ map. Given the sun's declination in level 1. PrjIll \<< \-> dec \<< 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 \-> i ilon ilat lilon lilat xt m x y z th lon lat s c @ Initialize the width table to "darkness". This was originally an @ array, but 64 real numbers took too much space, so I changed it to a @ string. 255 means no light, other numbers are 1/2 the width to @ illuminate. \<< "" 1 64 FOR c 255 CHR + NEXT @ Flag 6 is set the first time through, so interpolation isn't attempted. 'WTAB' STO 6 SF dec @ Build the transformation for the declination. D\->R NEG SIN 's' STO dec D\->R NEG COS 'c' STO \pi \->NUM 2 / NEG 'th' STO @ Increment over a semicircle of illumination. DO s NEG @ Transform the point through the declination rotation. th SIN * 'x' STO th COS 'y' STO c th SIN * 'z' STO y x @ Transform the resulting coordinate through the map projection to @ obtain screen coordiantes. ATAN2 R\->D 'lon' STO z ASIN R\->D 'lat' STO 62 lat 90 + .344444444444 * - IP 'ilat' STO lon .35 * IP 'ilon' STO IF 6 FS?C THEN @ First time through; just record the previous coordinates. ilon 'lilon' STO ilat 'lilat' STO ELSE @ Interpolate between the current and previous points if we've jumped @ more than one map line down. IF lilat ilat == THEN WTAB 62 ilat - ilon 1 MAX CHR REPL 'WTAB' STO ELSE ilon lilon - ilat lilat - / 'm' STO lilat 'i' STO WHILE i ilat \=/ REPEAT lilon i lilat - m * .5 + FLOOR + 'xt' STO WTAB 62 i - xt 1 MAX CHR REPL 'WTAB' STO ilat lilat - SIGN 'i' STO+ END END ilon 'lilon' STO ilat 'lilat' STO END \pi @ Back at the outer loop; give the user feedback about how the @ computation is progressing. \->NUM Res / 'th' STO+ "[" th \pi \->NUM 2 / + \pi \->NUM / 100 * 100 MIN FLOOR \->STR + "%]" + 7 DISP @ If the user wishes to terminate this computation, record the fact @ and abort the loop. IF KEY THEN DROP 3 'th' STO -1 'STAT' STO 1000 'OldDec' STO END @ Keep computing until the full semicircle is done. UNTIL th \pi \->NUM 2 / .001 + > END @ One edge of the map is fully illuminated, but may have been missed @ in the above loop, so tweak the widths at that end. IF dec 0 < @ South pole in daylight THEN -1 64 ELSE 1 1 @ North pole in daylight. END 'ilat' STO 'lilat' STO @ Don't bother if the user's already aborting. IF STAT 0 \>= @ Loop from the edge of the map until a non-dark line is found, @ illuminating as you go. THEN ilat 31 FOR j IF WTAB j DUP SUB NUM 255 \=/ THEN DO WTAB j 63 CHR REPL 'WTAB' STO IF j ilat == THEN lilat 'j' STO END lilat NEG 'j' STO+ UNTIL j 0 == END 31 'j' STO END lilat STEP @ All done; keep track of the declination value used in this @ computation, so it can be avoided next time if the current @ declination is close enough. dec 'OldDec' STO END \>> \>> \>> @ Calculate the position of the sun. Level 2 contains the Julian time @ of the instant for which the position is desired, and level 1 should @ be nonzero if the apparent position (corrected for nutation and @ aberration) is desired. @ The Sun's longitudinal position (in degrees) is returned in level 4; @ divide by 15 to get hours. The declination is returned in level 3. @ The radius vector (in astronomical units) is returned in level 2, @ and the Sun's longitude (true or apparent, as desired) is returned @ as degrees in level 1. SunPos \<< \-> jd ap \<< 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 \-> t t2 t3 l m e ea v th om eps ra dec rv slong \<< jd 2415020 - 36525 / @ Convert time to Julian centuries of 36525 ephemeris days, measured @ from the epoch 1900 January 0.5 ET. DUP 't' STO DUP DUP * DUP 't2' STO * 't3' STO 279.69668 36000.76892 t * + .0003025 t2 * + FixAng 'l' STO @ 'l' gets the geometric mean longitude of the Sun, referred to the @ mean equinox of the date. Now compute the sun's mean anomaly. 358.47583 35999.04975 t * + .00015 t2 * - .0000033 t3 * - FixAng 'm' STO @ Now take the eccentricity of the Earth's orbit into account. .01675104 .0000418 t * - .000000126 t2 * - 'e' STO m e @ Compute the eccentric anomaly. Kepler 'ea' STO 1 e @ Compute the true anomaly. + 1 e - / \v/ ea 2 / TAN * ATAN R\->D 2 * FixAng 'v' STO l v @ Compute the sun's true longitude. + m - 'th' STO @ Obliquity of the ecliptic: 23.452294 .0130125 t * - .00000164 t2 * - .000000503 t3 * + 'eps' STO IF ap @ Corrections for sun's apparent longitude, if desired. THEN 259.18 1934.142 t * - FixAng 'om' STO .00569 .00479 om D\->R SIN * - NEG 'th' STO+ .00256 om D\->R COS * 'eps' STO+ END th @ Record sun's longitude and radius vector, then compute coordinates. 'slong' STO 1.0000002 1 e e * - * 1 e v D\->R COS * + / 'rv' STO eps D\->R COS th D\->R SIN * th D\->R COS ATAN2 R\->D FixAng 'ra' STO eps D\->R SIN th D\->R SIN * ASIN R\->D 'dec' STO ra dec rv slong \>> \>> \>> @ Solve the equation of Kepler. Kepler \<< \-> m ecc \<< 0 0 \-> e d \<< m D\->R DUP 'm' STO 'e' STO DO e e SIN ecc * - m - 'd' STO d 1 e COS ecc * - / NEG 'e' STO+ UNTIL d ABS .000001 \<= END e \>> \>> \>> @ Compute arctan(y/x), with y in level 2 and x in level 1. Returns @ appropriate quadrant treating y and x as rectangular coordinates @ being converted to polar coordinates. ATAN2 \<< SWAP \->V2 -16 SF V\-> -16 CF SWAP DROP \>> @ Translate a big angle back in to the range 0-360 degrees. FixAng \<< DUP 360 / FLOOR 360 * - \>> @ Show the currently chosen date and time and prompt the user for @ their next action (used by the "THEN" submenu). ShwTm \<< CLLCD DT TM TSTR 4 DISP "Choose time, push GO" 3 DISP 2 FREEZE \>> GMST \<< \-> jd \<< jd .5 + FLOOR .5 - 2415020 - 36525 / 0 \-> t th0 \<< 6.6460656 2400.051262 t * + .00002581 t t * * + 'th0' STO jd .5 + DUP FLOOR - 24 * 1.002737908 * th0 + DUP 24 / FLOOR 24 * - \>> \>> \>> @ Convert the chosen date and time to Greenwich time. GTime \<< RCLF -42 CF DT 'GD' STO TM TZ - DUP IF 24 \>= THEN 24 - GD 1 DATE+ 'GD' STO END 'GT' STO STOF \>> @ Compute the chosen Julian time, as a date and day fraction. JTime \<< JDate .5 - GT HMS\-> 24 / + \>> @ Compute the Julian date for the chosen date. JDate \<< GD DUP IP SWAP FP 100 * DUP IP SWAP FP 10000 * 0 \-> m d y c \<< IF 'm>2' THEN -3 'm' STO+ ELSE 9 'm' STO+ -1 'y' STO+ END y 100 / IP DUP 'c' STO 100 * NEG 'y' STO+ d c 146097 * 4 / IP + y 1461 * 4 / IP + m 153 * 2 + 5 / IP + 1721119 + \>> \>> @ Holds the currently chosen time: TM 22.3138814208 @ Holds the currently chosen date: DT 11.151991 @ Greenwich time and date for above values: GT 3.3138814208 GD 11.161991 @ Stores the declination for the last computed daylight width table. OldDec -18.5917760031 @ The most recently computed width table, as a string for space @ reasons (about 70 bytes instead of 500, and there are two of these @ around while animation is active.) WTAB C$ 64 \255\255\255\255\255\255 !!!"""##$$$%&'')*+,/5????????? END