;-----------------------------------------------------------------------------
;
;       procedure:  solcor_sbsp
;
;       purpose:  do solar ephemeris
;
;       author:  unknown
;
;-----------------------------------------------------------------------------
        pro solcor_sbsp, djd, ras, decs, gsdt, bzro, p, solong, rsun

        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

            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
        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

        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