pro azam_solcor, djd, ras, decs, gsdt, bzro, p, solong, rsun ;----------------------------------------------------------------------------- ; ; procedure: solcor_sbsp ; ; purpose: do solar ephemeris ; ; author: unknown ; modified 17 Jan 2012 by BWL for processing vector,array data for azam ; ;----------------------------------------------------------------------------- savdjd = 1000000.d savras = 1000000.d savdecs = 1000000.d savgsdt = 1000000.d savbzro = 1000000.d savp = 1000000.d savsolong = 1000000.d savrsun = 1000000.d pi = 3.14159265358979323846d tpi = 2.d*pi ;if( djd eq savdjd ) then begin ;ras = savras ;decs = savdecs ;gsdt = savgsdt ;bzro = savbzro ;p = savp ;solong = savsolong ;rsun = savrsun whrdj = where( (djd - savdjd) eq 0.d,countj) if countj gt 0 then begin ras(whrdj) = savras(whrdj) decs(whrdj) = savdecs(whrdj) gsdt(whrdj) = savgsdt(whrdj) bzro(whrdj) = savbzro(whrdj) p(whrdj) = savp(whrdj) solong(whrdj) = savsolong(whrdj) rsun(whrdj) = savrsun(whrdj) return endif jd = long(djd) fjd = djd-jd d=(jd-2415020)+fjd iyr=long(d/365.25) fiyr=double(iyr) g=-.026601523+.01720196977*d-1.95e-15*d*d-tpi*fiyr xlms=4.881627938+.017202791266*d+3.95e-15*d*d-tpi*fiyr obl=.409319747-6.2179e-9*d ecc=.01675104-1.1444e-9*d e=ecc*sin(g)/(1.-ecc*cos(g)) e=g+e-0.5*e^3 rsun=961.18/(1.-ecc*cos(e)) f=d-365.25*fiyr gsdd=1.739935476+(tpi*f+1.342027e-4*d)/365.25 gsdt=gsdd+tpi*(fjd-0.5) xlts=xlms+2.*ecc*sin(g)+1.25*ecc*ecc*sin(2.*g) sndc=sin(xlts)*sin(obl) decs=asin(sndc) csra=cos(xlts)/cos(decs) ras=acos(csra) ; if(sin(xlts) lt 0.) then ras=tpi-ras whrsi = where( sin(xlts) lt 0.,counti) if counti gt 0 then ras(whrsi) = tpi - ras(whrsi) omega=1.297906+6.66992e-7*d thetac=xlts-omega bzro=asin(.126199*sin(thetac)) p=-atan(cos(xlts)*tan(obl))-atan(.127216*cos(thetac)) ;xlmm=atan2(.992005*sin(thetac),cos(thetac)) xlmm=atan(.992005*sin(thetac),cos(thetac)) jdr=jd-2398220 irot=long((jdr+fjd)/25.38) frot=(jdr+fjd)/25.38-irot solong=xlmm-tpi*frot+pi-3.e-4 ;if(solong lt 0.) then solong = solong+tpi ;if(solong ge tpi) then solong = solong-tpi whrso = where (solong lt 0.,counts) if counts gt 0 then solong(whrso) = solong(whrso)+tpi whrso = where (solong ge tpi,counts) if counts gt 0 then solong(whrso) = solong(whrso)-tpi savdjd = djd savras = ras savdecs = decs savgsdt = gsdt savbzro = bzro savp = p savsolong = solong savrsun = rsun ; print, djd, ras, decs, gsdt, bzro/!dtor, p/!dtor, solong, rsun ; print,'Julian date: ',djd ; print,'right ascension disk center(deg): ',ras/!dtor ; print,'declination disk center(deg): ',decs/!dtor ; print,'right ascension from central meridian(deg): ',gsdt/!dtor ; print,'b0-angle(deg): ',bzro/!dtor ; print,'p-angle(deg): ',p/!dtor ; print,'carrington longitude disk center(deg): ',solong/!dtor ; print,'solar radius (arcsecs): ',rsun return end