PRO lpsun,iy,im,id,sday,ha,dec ; ;returns hour angle and declination of the sun at La Palma SSO for a set of ;times ; ; time variable for Newcomb's folmulae: ; fraction of Julian centuries elapsed since 1900.05 (noon 1st of January) ; = J.D.2415020.0) jd=julian(iy,im,id) h0=(jd-2415020.0)/36525.0 h=(jd+double(sday)/86400.D0-2415020.0)/36525.0 hh=h*h ; Newcomb's formulae. (page 98 Explanatory suppl. to the ephemeris) ; mean obliquity of the ecliptic ehel=0.4093198-2.2703E-4*h-2.86E-8*hh ;type,'mean obliquity =',ehel ; eccentricity of earth's orbit eks=0.01675104-0.0000418*h ;type,'eccent. =',eks ; mean longitude of sun sml=279.6967+36000.769*h ; sml= ( (sml/360.0) - (fix(sml/360.0)) )*360.0 sml = MODULO(sml,360.) ;type,'mean longitude of sun =',sml sml=sml*!pi/180.0 ; mean anomaly anm=358.4758+35999.0498*h-0.00015*hh anm=anm*!pi/180.0 ; true longitude of sun (sl) cc=(1.91946-0.00479*h)*sin(anm)+0.020*sin(2*anm) cc=cc*!pi/180 sl=sml+cc ; declination of apparent sun dec=asin(sin(ehel)*sin(sl)) ; right ascension of apparent sun ra=atan((cos(ehel)*sin(sl)),cos(sl)) ; convert ra in radians to hours ra=ra*24.0/(2.*!pi) ra=ra+24.*(ra lt 0) ; convert dec in radians to degrees dec=dec*360.0/(2.*!pi) ;sidereal time ; sidereal time in Greenwich at 0 UT (Newcomb) sidg0=6.6461D0+2400.0513D0*h0 sidg0 = MODULO(sidg0,24.) ;type,'sidereal time in Greenwich at 0 UT',gmt(sidg0*3600.),sidg0 ; sidereal time in Greeenwich at any instant sidg=(double(sday)*1.002737908D0/3600.0D0)+sidg0 ;type,'computed sidereal time at Greenwich',dms(sidg%24),sidg,sidg%24 ; longitude of observatory in degrees, negative if west lo=-17.880 ; convert longitude of obs (degrees) to time measure and find local sidereal ; time, longitude is positive when west sidl=sidg+lo*(24.0/360.0) ha=sidl-ra ;restrict to the -12 to +12 range ha = MODULO(ha+12.,24.) - 12. ha = MODULO(ha-12.,24.) + 12. ;print,'hour angle at Greenwich ',dms(sidg-ra) ;print,'hour angle at La Palma ',dms(ha) end