( SolarCalc)

-1 DPL !
2 FLD !
: deg->rad ( fp degrees...fp radians)
  Pi F* 180 FINT F/ ;

: rad->deg 180 FINT F* Pi F/ ;

: .deg rad->deg F. ;

50.7 FLOAT deg->rad VARIABLE Latitude
-3.4 FLOAT deg->rad VARIABLE Longitude
0 VARIABLE BST
0 VARIABLE clock_hr    0 VARIABLE clock_min
12 ARRAY sum_month

: load_sum_month
   0 31 59 90 120 151 181 212 243 273 304 334
    0 11 DO I sum_month ! -1 +LOOP ;
load_sum_month

( keyboard entry stuff repeated here - temp)
0 VARIABLE month  0 VARIABLE day 
0 VARIABLE hour_angle  0 VARIABLE declination
0 VARIABLE sol_elevation  0 VARIABLE azimuth

: input QUERY INTERPRET ;

: input_date&time
  CR ."  month (1-12):" input month !
  CR ."    day (1-31):" input day !
  CR ."   hour (6-18):" input clock_hr !
  CR ."     min (<60):" input clock_min !
;

( new words recommence here:)
( clock time ... local time wrt solar noon)
: localT  ( longitude,hr,min...local)
  FINT 60 FINT F/ SWAP FINT F+  ( longitude ,fp hour...)
  BST @ IF 1 FINT F- ENDIF   ( =GMT hour)
  12 FINT F-    ( clock time hrs wrt GMT noon)
  SWAP rad->deg 15 FINT F/
  F+ ;

: hour_ang Longitude @ clock_hr @ clock_min @ localT -15 FINT F* 
   deg->rad hour_angle  ! ;

: calc_N   ( nmonth,nday....N) ( N is #days after Jan 1st)
  SWAP 1- 0 11 confine sum_month @  + ;

: calc_declination ( N...dec)  ( dec in rads)
  10 + 2* FINT Pi  F* [ 365 FLOAT ] LITERAL F/ FCOS
  [ -24.44 FLOAT ] LITERAL F* deg->rad declination ! ;

( solar elevation - all f.p.rads)
: sol_ht  ( Latitude,declination,hour-angle....solar elevation)
  FCOS    ( Lat,dec,COShr)
  >R                          ( Lat,dec...)                  ( R: Coshr)
  DUP FSIN SWAP FCOS          ( Lat,SINdec,Cosdec...)        ( R: Coshr)
  ROT DUP FSIN SWAP FCOS      ( SINdec,COSdec,SINLat,COSlat)
  ROT                         ( SINdec,SINlat,COSlat,COSdec)
  R> F* F* >R                 ( SINdec,SINlat...)
  F* R> F+ 
  FASN sol_elevation ! ;

: calc_azimuth   ( declination,hr_ang,ht...azimuth)
  FCOS >R
  FSIN SWAP FCOS F*
  R> F/ FASN FMNF azimuth ! ;
 
: calc_solposn   ( ...)
  hour_ang
  month @ day @ calc_N calc_declination
  Latitude @  declination @ hour_angle @ sol_ht
  declination @ hour_angle @ sol_elevation @ calc_azimuth
  ;

: show_result
 CR ." hour_angle=" hour_angle @ .deg SPACE
 ." Declination=" declination @ .deg
 CR ." Sol Elevation=" sol_elevation @ .deg SPACE
 ." Azimuth=" azimuth @ .deg CR ;

: sun_point
 input_date&time
 calc_solposn 
;

: scan_hrs  ( hr1,hr2....)
  CR ."  month (1-12):" input month !
  CR ."    day (1-31):" input day !
  1+ SWAP DO I DUP CR . clock_hr ! 0 clock_min ! 
             calc_solposn show_result
             60 10 DO I clock_min ! ." +" I . ." min:"
                      calc_solposn show_result
                10 +LOOP
              LOOP ;







