NB. Program to get altitude & azimuth of sun, moon, and planets.
NB. Also plots lunar crescent in correct phase & orientation.
NB. Usage: Skys''   <--- program will prompt for other info
NB. Requires data files in SS_DAT. This is only accurate to ~ 1 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.
NB. This version (Skys as opposed to Sky) plots all stars in the Bright
NB. Star Catalog brighter than 5.2 magnitudes and with altitudes and
NB. azimuths between user supplied limits, and also plots the location
NB. of the sun (yellow), moon (blue) and planets (red and orange).
NB. This uses the new {{ ... }} syntax. Run with j602 (or later) version of J.
load 'local_defs.ijs'
require 'files'
cn=. *@(2&#.@(*@])) * |@] NB. convert _d m s -> _d _m _s
cno=. cn@([: 0 1 2&{ ] , 0 0"_) NB. add zeros if s or m & s missing
dtr =: ((o.%648000)&*)@(360 60 60&#.)@cno f. NB. Convert deg m s -> rad
rtd =: **360 60 60&#:@((648000%(o.1))&*)@| NB. rad -> deg m s
NB. mread=: (0&".)@('m'&fread)
RtD=: %DtR=: 1p1%180
require 'plot'
Skys=: monad define  NB. ----------- Start "Skys" -----------------
'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 '/your_path/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
xx=. (*HAr)* rth | HAr
print 'Planet  ';'Right Ascension ';'Declination ';'Rho   ';'Hour Angle'
all=. snames;(<"1 rth RAr);(<"1 rtd Dec);(<"0 Rho);(<<"1 xx)
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<Az360<2pi 
test=. -/(0 3){ Az360  NB. diff in azimuths for crescent orientation below
EW=.'West';'East'
EWx=. <(HAr<0){EW
yy=. ,<"1 rtd Alt
zz=. ,<"1 rtd Az
all=. snames;yy;zz;EWx
print' '
print 'Planet';'Altitude  ';'Azimuth  ';'(East or  West)'
i=. _1
while. (i=. >: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)
xx=. cos (3{RAr)-(0{RAr)
M_S_ang=. acos ((sin 0{Dec)*(sin 3{Dec))+((cos 0{Dec)*(cos 3{Dec)*xx)
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
xx=. ((cos zds)-(cos M_S_ang)*(cos zdm))%((sin zdm)*(sin M_S_ang))
ocr=. 1p1 - acos xx
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'
NB. --------- stars ---------------------------------------------
'r d m'=: |: mread'/your_path/SS_DAT/stars5.2.dat'
NB. r=Right Ascwntion, d=Declination, both in degrees; m=magnitude
sd=. DtR*d   NB. sd = stellar declination in radians
sr=. DtR*r   NB. sr = stellar right ascension in radians
NB. precess J2000 star coordinates to coordinates of date:
PMAT=: Get_PMAT jul- 2400000.5       NB. precession matrix
'sr sd'=. |: PMAT&precess"1 |:>sr;sd
sHAr=. LSTr -sr                  NB. stellar hour angles in radians
sHAr=. sHAr + 2p1*(sHAr<_1p1)
sHAr=. sHAr - 2p1*(sHAr>1p1)     NB. -12h < HA < 12h
szd=. acos sqz=. ((sin latit)*(sin sd))+(cos latit)*(cos sd)*(cos sHAr)
NB. szd = stellar zenith distances in radians, sAlt = stellar altitudes
sAlt=. 0.5p1- szd=. szd+ (szd<0)*1p1
top=. (sin sd)- sqz* sin latit
sAz=: acos top% (cos latit)*sin szd
NB. HA positive => Azimuth must be 180<Az<360
sAz360=. (sAz*(sHAr<:0))+(2p1-sAz)*(sHAr>0)  NB. 0<sAz360<2pi 
print'  Enter Altitude limits A1 & A2 in degrees:'
'Alt1 Alt2'=. DtR* enter''
NB. 'Alt1 Alt2 Az1 Az2'=. 0,0.4p1,0.5p1,1.25p1
id1=. (Alt1<sAlt)*(Alt2>sAlt)
al0=: -: RtD* Alt1 + Alt2    NB. central altitude in degrees
print'  Enter Azimuth limits Az1 & Az2 in degrees:'
'Az1 Az2'=. DtR* enter''
id2=. (Az1<sAz360)*(Az2>sAz360)
az0=: -: RtD* Az1 + Az2      NB. central azimuth in degrees
id3=. id1*id2                NB. indices within alt & az limits
NB. ----------------------------------------------------------------
sal=: RtD*(I. id3){sAlt    NB. stellar altitudes in degrees
saz=: RtD*(I. id3){sAz360  NB. stellar azimuths in degrees
'psal psaz'=. sal;saz
'Gpsal Gpsaz'=: (al0,az0)&Gnomo (>sal;saz)
m=. (I. id3){m 
i0=. I. m> 4.2
i1=. I. (m> 3.2)*(m<: 4.2)
i2=. I. (m> 2.2)*(m<: 3.2)
i3=. I. (m >: 1.2)*(m <: 2.2)
i4=. I. m < 1.2
id3=. (Alt1<Alt)*(Alt2>Alt)   NB. reuse id3, id4 for planets
id4=. (Az1<Az360)*(Az2>Az360)
id5=. id3*id4         NB. indices within alt & az limits
xal=. (I. id5){Alt    NB. visible planet altitudes in radians
xaz=. (I. id5){Az360  NB. visible planet azimuths in radians
NB. --------------------- Make Plot -----------------------------------
null=. enter''
pd 'reset'
pd 'type point'
if. (1<:+/id5) do.
  pd 'pensize 2.0;color orange'   NB. make planets initially orange
  'pxal pxaz'=. (RtD*xal);(RtD*xaz)
  pd pxaz;pxal
end.
  pd 'pensize 3.0;color red'   NB. make Venus & Jupiter (2 & 5) red
if. (5{id5)=1 do.              NB. Jupiter
  'pxaz pxal'=. (5{Az360);(5{Alt)
  'pxal pxaz'=. (RtD*pxal);(RtD*pxaz)
  pd (1$pxaz);(1$pxal)         NB. make position an array of one item
end.
if. (2{id5)=1 do.              NB. Venus   
  'pxaz pxal'=. (2{Az360);(2{Alt)
  'pxal pxaz'=. (RtD*pxal);(RtD*pxaz)
  pd (1$pxaz);(1$pxal)
end.
if. (3{id5)=1 do.                 NB. if Moon then do
  'pxaz pxal'=. (3{Az360);(3{Alt)
  'pxal pxaz'=. (RtD*pxal);(RtD*pxaz)
   aspect=: (Alt2-Alt1)%(Az2-Az1) NB. altitude extent divided by azimuth
   screen=: 0.7566                NB. height of display divided by width
   pd 'pensize 6;color white'     NB. overwrite moon position with white
  pd (1$pxaz);(1$pxal)           
  pd 'pensize 0.5;color blue'     NB. draw lunar crecent in blue
NB.  pd (pxaz+2*screen*sxra);(pxal+2*aspect*syra)
  pd (pxaz+screen*sxra);(pxal+aspect*syra)
NB.   'Gmoonx Gmoony'=: (al0,az0)&Gnomo (>(pxal+2*syra);(pxaz+2*sxra))
end.
pd 'pensize 6.0;color yellow'     NB. overwrite sun with yellow
if. (0{id5)=1 do.
  'pxaz pxal'=. (0{Az360);(0{Alt)
  'pxal pxaz'=. (RtD*pxal);(RtD*pxaz)
   pd (1$pxaz);(1$pxal)
end.
pd 'pensize 1.0;color black'
pd (i0{psaz);(i0{psal)
pd 'pensize 1.5;color black'
pd (i1{psaz);(i1{psal)
pd 'pensize 2.0;color black'
pd (i2{psaz);(i2{psal)
pd 'pensize 2.5;color black'
pd (i3{psaz);(i3{psal)
pd 'pensize 3.0;color black'
pd (i4{psaz);(i4{psal)
pd 'pdf'
NB. ------------------- End of Main Verb "Skys" ---------------------
)

NB. Julian date. Usage: julian d m y
julian=: 3 : 0 "1 
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
)

NB. Modified Julian date: Mod. JD = JD - 2400000.5
NB. see Green p 251 (who dropped a 0)
mjd=: monad :'_2400000.5+ julian y'

daytime=: {{ 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'0) Pikesville, MD'
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. 0 do. lonlat=. 76.711000   39.374556
case. 1 do. lonlat=. 76.937000   38.962139
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=: {{ snames=: ''
delbar=: 0$0
ldxs=: 0$0
VUW=: 0 3$0
i=. 0
while. (i=. >:i) <: 7 do.
  file=. '/your_path/jph/J6/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
xx=. < ;. 1 (>0{ y)
sname=: > 1 { xx
NB. Delta-bar and longitude index:
'del ldx'=: ". > 2 3 { xx
xx=. ". , (>1{ y)
nterm=.  0 1 2 { xx     NB. # of V,U & W terms
ncol=.  3 { xx          NB. # number of columns
cindx=. (4+i. ncol){xx  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 
) 

NB. Routine to precess star coordinates from J2000 
NB. to coordinates of date. RA & Dec in radians.
NB. Usage: (precession matrix) precess RA,Dec       
precess=: dyad define
'ra dec'=. y
v=. ((cos dec)*(cos ra),(sin ra)),sin dec 
'v0 v1 v2' =. x mm v 
a=. (o.(v0<0))+ atan v1%v0
a, asin v2 
)

Get_PMAT=: monad define  NB. get precession matrix for epoch of date
dmy1=. 'd1 m1 y1'=. 1.5 1 2000   NB. J2000.0 (epoch of star catalog)
'eF eD'=. (mjd dmy1), y
NB. Base epoch is J2000.0 = 2451545.0 <-> mjd 1.5 1 2000 = 51544.5
epoch=. 36525%~(eF-51544.5),(eD-eF)   NB. T & t in Julian centuries
Pr_Mat epoch
)

NB. Precession Matrix (Seidelmann "Explanatory Supplement" pp. 103-104)
Pr_Mat=: monad define  NB. Th = theta, Z = z_A, Z0 = Zeta_A 
'T t'=. y              NB. T = Fixed - 2000, t = Date - Fixed
NB. zeta_A, z_A and theta_A from p 104, Table 3.211.1:
Z0=. 0,(2306.2181 1.39656 _0.000139 p. T),(0.30188+_0.000344*T),0.017998
cZ0=. cos Z0 =. dtr 0,0,(Z0 p. t)
Z=. 0,(2306.2181 1.39656 _0.000139 p. T),(1.09468+0.000066*T),0.018203
cZ=. cos Z =. dtr 0,0,(Z p. t)
Th=. 0,(2004.3109 _0.85330 _0.000217 p. T),(_0.42665+_0.000217*T),_0.041833
cT=. cos Th =. dtr 0,0,(Th p. t)
'sZ0 sZ sT'=. sin Z0, Z, Th
NB. matrix P from p 103, eq (3.21-8):
f=.(cZ0*cZ*cT)-sZ0*sZ
f=.f,(-(cZ0*sZ)+sZ0*cZ*cT),(-cZ*sT)
f=.f,(sZ0*cZ)+cZ0*sZ*cT
f=.f,(cZ0*cZ)-sZ0*sZ*cT
3 3$f,(-sZ*sT),(cZ0*sT),(-sZ0*sT),cT
)

moon_plot=: monad define
'phase ang'=. y
sx=. sin xx=. }: 1p1*int01 2000
xa=. sx ,|. phase*sx
ya=. (,|.) cos xx
sth=. sin ang
cth=. cos ang
sxra=: xra=. (cth*xa)+(sth*ya)   NB. make sxra,syra global 
syra=: yra=. (-sth*xa)+(cth*ya)  NB. for use in sky plot
options=. 'title The Moon;axes 0 0;grids 0 0;labels 0;aspect 1.0'
NB. options=. options,';pensize 3;frame 1;xrange _1 1;yrange _1 1' 
options=. options,';pensize 8;frame 1;xrange _1 1;yrange _1 1' 
pd 'reset'
pd options
pd xra;yra
pd 'pdf'
)

NB. Gnomonic projection
NB. Tangent plane centered at al0,az0
NB. Altitude and Azimuth of stars are al & az
Gnomo=: {{ 'al0 az0'=. DtR* x
'al az'=. DtR* y
cosc=. (sin al0)*(sin al)+(cos al0)*(cos al)*cos az-az0
xx=. cosc%~ (cos al)*sin az-az0
yy=. cosc%~ ((cos al0)*(sin al))-((sin al0)*(cos al)*cos az-az0 )
>xx;yy }}