NB. Program to get altitude & azimuth of sun, moon, and planets. NB. Also plots lunar crescent in correct phase & orientation. NB. Usage: Sky'' <--- program will prompt for other info NB. Requires data files in SS_DAT. This is only accurate to ~arcmin NB. The computation is based on trignometric series given in a paper by NB. Van Flandern, T.C., and Pulkkinen, K.F. 1979. Ap.J.Supp., 41, 391. load 'local_defs.ijs' require 'files' mread=: (0&".)@('m'&fread) NB. read a numeric array require 'plot' Sky=: monad define 'day mth yr hr'=. daytime'' print day,mth,yr,hr jul=. julian day,mth,yr ut=. hr + 5 jd=. jul + ut%24 if. (ut >: 24.0) do. ut=. ut - 24 jul=. jul + 1 end. dt=. 36525%~ jul-2415020 print 'Julian day = '; precise jul theta=. 0.276919398+(100.0021359+1.075e_6*dt)*dt STr=. 2p1*(24%~1.002737908*ut)+ 1|theta 'longit latit'=. place'' LSTr=: STr - longit LST=: rth LSTr print 'LST = '; LST T=: jd - 2451545.0 NB. days from 2000 Jan 1.5 tc=: 1 + T%36525 NB. centuries from 1900.0 NB. coefficients of fundamental parameters: 'c1 c2'=. |: mread 'SS_DAT/ss0.dat' A=: 2p1* 1|(c1+ T*c2) NB. L,F, and G for each planet VUW=. Get_VUW'' snames=: ;: snames L=. (<: ldxs) { A 'V U W'=. |: VUW RAr=: L + asin(W % %:(U - *:V)) Dec=: asin V % %:U Rho=. delbar * %:U HAr=: LSTr - RAr HAr=: HAr + 2p1*(HAr<_1p1) HAr=: HAr - 2p1*(HAr>1p1) NB. -12h < HA < 12h NB. ---------- Lunar Parallax Correction ------------ NB. 'del_ha del_dec'=. Parallax (3{HAr),(3{Dec),latit,(3{Rho) 'del_ha del_dec'=. Parallax (>3{ &. > HAr;Dec;Rho),latit xx=. (*del_ha)* rth |del_ha print 'Lunar parallax corrections: ';xx;(rtd del_dec) HAr=: HAr(3})~ (3{HAr)+ del_ha RAr=: RAr(3})~ (3{RAr)- del_ha Dec=: Dec(3})~ (3{Dec)+ del_dec x=. (*HAr)* rth | HAr print 'Planet ';'Right Ascension ';'Declination ';'Rho ';'Hour Angle' all=. snames;(<"1 rth RAr);(<"1 rtd Dec);(<"0 Rho);(<<"1 x) i=. _1 while. (i=. >:i)<: 6 do. print >i{ &. >all end. null=. enter'' NB. ---- Altitude and Azimuth --------- zd=. acos qz=. ((sin latit)*(sin Dec))+(cos latit)*(cos Dec)*(cos HAr) Alt=: 0.5p1- zd=. zd+ (zd<0)*1p1 'zdm zds'=. (3 0){ zd NB. moon & sun for moon phase below top=. (sin Dec)- qz* sin latit bot=. (cos latit)*sin zd Az=. acos top%bot Az360=. (Az*(HAr <: 0))+(2p1-Az)*(HAr > 0) NB. 0:i)<: 6 do. print >i{ &. >all end. null=. enter'' NB. ----------- Orientation of lunar crescent -------------------- print 'Zenith distance of moon ';(rtd zdm);' and of sun ';(rtd zds) x=. cos (3{RAr)-(0{RAr) M_S_ang=. acos ((sin 0{Dec)*(sin 3{Dec))+((cos 0{Dec)*(cos 3{Dec)*x) phase=. -: 1 + cos 1p1 - M_S_ang if. test<0 do. test=. test + 2p1 end. NB. sun azimuth measured from moon if. test < 1p1 do. waxwan=. ', waxing' else. waxwan=. ', waning' end. print 'Moon-Sun angle ';(rtd M_S_ang);' Phase ';phase; waxwan x=. ((cos zds)-(cos M_S_ang)*(cos zdm))%((sin zdm)*(sin M_S_ang)) ocr=. 1p1 - acos x if. test<0 do. test=. test + 2p1 end. NB. sun azimuth measured from moon if. test > 1p1 do. ocr=. - ocr end. NB. if relative azimuth > pi, sun to left moon_plot (1-+:phase),(0.5p1-ocr) 'Nadir-Moon-Sun angle ';(rtd ocr);'(+)-> sunny side right, (-)-> left' ) julian=: 3 : 0 " 1 NB. Usage: julian d m y m1=. m+12*i=. 3<~ m=. 1{y y1=. (-i)+(y0<0)+ y0=. 2{y i=. 0<: #. |. *y -15 10 1582 b=. i*(2+(-a)+<. 4%~a=. <.y1%100) c=. (*c)* <. |c=.(365.25*y1)-0.75*(y1<0) 1720994.5+b+c+(0{y)+<.30.6*m1+1 ) daytime=: monad define print 'Enter: day month year' print '(Or hit enter for current date & time)' dmy=. enter'' if. (#dmy)=0 do. 'day mth yr hr min sec'=. 2 1 0 3 4 5 { 6!:0'' print day; mth; yr print hr; min; sec hr=. hr + 60 %~ (min + sec % 60) else. 'day mth yr'=. dmy print 'Enter: hour minute (EST/EDT)' 'hr min'=. enter'' hr=. hr + 60 %~ min end. print 'Is this Daylight Savings Time? (y/n):' key=. ([: 1!:1 1:)'' hr=. hr - (key='y') day,mth,yr,hr ) place=: monad define print'Where in the world are you?' print'1) College Park, MD' print'2) Salem, Ohio' print'3) Buffalo, NY' print'4) Newark, Ohio' print'5) Flint Hill, VA' print'6) Somewhere else...' key=. enter'' select. key case. 1 do. lonlat=. 76.954166 39.001666 case. 2 do. lonlat=. 80.863997 40.910000 case. 3 do. lonlat=. 79.000000 43.000000 case. 4 do. lonlat=. 82.447222 40.052778 case. 5 do. lonlat=. 78.0755 38.7785 case. 6 do. print'Enter longitude & latitude (degrees; West +, East -)' lonlat=. enter'' end. (1p1%180)*lonlat ) Get_VUW=: monad define snames=: '' delbar=: 0$0 ldxs=: 0$0 VUW=: 0 3$0 i=. 0 while. (i=. >:i) <: 7 do. file=. 'SS_DAT/ss',(":i),'.dat' data=. < ;. _2 fread file VUW=: VUW, Series data snames=: snames, sname delbar=: delbar, del ldxs=: ldxs, ldx end. VUW ) Series=: monad define x=. < ;. 1 (>0{ y) sname=: > 1 { x NB. Delta-bar and longitude index: 'del ldx'=: ". > 2 3 { x x=. ". , (>1{ y) nterm=. 0 1 2 { x NB. # of V,U & W terms ncol=. 3 { x NB. # number of columns cindx=. (4+i. ncol){ x NB. column indices start=. +/\2,}:nterm Vblock=. ". >((0{start)+i.(0{nterm)){y Ublock=. ". >((1{start)+i.(1{nterm)){y Wblock=. ". >((2{start)+i.(2{nterm)){y Coef=. 0{"1 Vblock Tcoef=. tc^(1{"1 Vblock) sincos=. (2{"1 Vblock)&o. Trig=. 3}."1 Vblock AA=. (<: cindx){ A Arg=. +/"1 AA *"1 Trig V=. +/ Coef*Tcoef*sincos Arg Coef=. 0{"1 Ublock Tcoef=. tc^(1{"1 Ublock) sincos=. (2{"1 Ublock)&o. Trig=. 3}."1 Ublock Arg=. +/"1 AA *"1 Trig U=. +/ Coef*Tcoef*sincos Arg Coef=. 0{"1 Wblock Tcoef=. tc^(1{"1 Wblock) sincos=. (2{"1 Wblock)&o. Trig=. 3}."1 Wblock Arg=. +/"1 AA *"1 Trig W=. +/ Coef*Tcoef*sincos Arg V,U,W ) Parallax=: monad define 'ha dec rho lat'=. y cq=. rho %~ cos lat sq=. rho %~ sin lat cha=. cos ha cd=. cos dec sd=. sin dec td=. sd%cd del_ha=. xx=. atan(cq*sin ha)%(cd - cq*cha) xx=. (cos ha+xx)*(sd-sq)%(cd*cha)-cq del_dec=. atan(xx-td)%(1+xx*td) del_ha,del_dec ) moon_plot=: monad define 'phase ang'=. y sx=. sin x=. }: 1p1*int01 2000 xa=. sx ,|. phase*sx ya=. (,|.) cos x sth=. sin ang cth=. cos ang xra=. (cth*xa)+(sth*ya) yra=. (-sth*xa)+(cth*ya) options=. 'title The Moon;axes 0 0;grids 0 0;labels 0;aspect 1.0' options=. options,';pensize 3;frame 1;xrange _1 1;yrange _1 1' pd 'reset' pd options pd xra;yra pd 'show' )