NB. For a given latitude, find the angle A between NB. the meridian and the plane of the Milky Way, NB. i.e., the +/- azimuth of the Galactic plane, NB. when the Milky Way passes through the zenith. NB. Also, find the times of the zenith crossings. NB. (N Galactic pole: RA=192.859 deg, Dec=27.128 deg) NB. => latitude must be < 62.872 deg (i.e., 90-27.128) NB. Usage: Milky_Way'' Milky_Way=: monad define 'day mth yr hr'=. daytime'' now=. htr hrx NB. "now" includes the DST offset 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=. (jul - 2415020)%36525.0 theta=. 0.276919398+(100.0021359+0.000001075*dt)*dt x=. theta - <. theta STr=. 2p1*x + 1.002737908*ut%24.0 'longit latit'=. place'' print 'Longitude & Latitude: ';(longit%dtrf);(latit%dtrf) LSTr=: STr - longit LST=: rth LSTr print 'Local Sidereal Time (LST) = '; LST inc=. dtrf* 62.872 NB. inclination of Galactic plane A=. rtd Ar=. asin (cos inc)% cos latit print 'Angle A (Azimuth of Milky Way)(+ or -) = '; A side_a=: asin (tan latit)% tan inc SWNE=. rth zSWNE=. (dtrf* 192.859)+ 0.5p1+side_a NWSE=. rth zNWSE=. (dtrf* 192.859)- 0.5p1+side_a capt1=. 'RA of Milky Way crossing zenith towards ' print (capt1,'NE'); SWNE print (capt1,'NW'); NWSE sidtosol=. >: %365.2422 NB. convert sidereal to solar HANE=. LSTr- zSWNE NB. H=LST-RA: HA of MW on zenith HANW=. LSTr- zNWSE 'T1 ap1'=. AMPM rth now - sidtosol* HANE 'T2 ap2'=. AMPM rth now - sidtosol* HANW capt2=. 'Time Milky Way crosses zenith toward North ' print (capt2,'East'); (}: T1);ap1 (capt2,'West'); (}: T2);ap2 ) 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 'Date: ';day;mth;yr;' Time:';hr;min;sec hrx=: hr+ 60%~(min+ sec% 60) else. 'day mth yr'=. dmy print 'Enter: hour minute (EST/EDT)' 'hr min'=. enter'' hrx=: hr+ 60%~ min end. key=. prompt_char'Is this Daylight Savings Time? (y/n):' hr=. hrx - (key='y') day,mth,yr,hr ) place=: monad define print'Where in the world are you?' print'1) College Park, MD 2) Salem, Ohio' print'3) Buffalo, NY 4) Newark, Ohio' print'5) Chillicothe, Ohio 6) Marietta, Ohio' print'7) Somewhere else...' key=. enter'' select. key case. 1 do. lonlat=. 76.954166 39.001666 case. 2 do. lonlat=. 80.8640 40.9100 case. 3 do. lonlat=. 78.8667 42.8833 case. 4 do. lonlat=. 82.4472 40.0528 case. 5 do. lonlat=. 82.9667 39.3167 case. 6 do. lonlat=. 81.4500 39.4000 case. 7 do. print'Enter longitude & latitude (degrees; West +, East -)' lonlat=. enter'' end. dtrf*lonlat ) NB. Convert 24h time to 12h AM/PM. AMPM=: monad define 'hr min sec'=. y ampm=. 'AM' if. hr >: 12 do. ampm=. 'PM' hr=. hr - 12 end. (hr,min,sec);ampm ) dtrf=: 1p1%180 NB. convert decimal degrees to radians