function soldis,julian

;----------------------------------------------------------------------------
;
;  Routine calculates the distance from the Earth to the Sun accounting for 
;  the mean anomaly of Earth's orbit at a particular date in the year.
;
; Input: 
;   julian....julian date as yyddd or yyyyddd,  it can handle either format 
;
; Output:
;   soldis....Earth - Sun distance in km

; Source: Mark Hervig
;
;----------------------------------------------------------------------------


;- some constants

   pi     = 3.1415926535898D0
   eccent = 0.01672592d0        ; Eccentricity of Earth orbit
   axsmaj = 1.4957D8            ; Semi-major axis of Earth orbit (km)
   solyr  = 365.2563855d0       ; Number of days in a solar year
   year1  = 1984.

;- get the date squared away,  we need date as yy.x,  where x = ddd/365.

   yy     = fix(julian/1000)
   dd     = fix(julian - yy*1000)

   if (yy le  99 and yy gt 50) then yy = 1900 + yy
   if (yy le  50)              then yy = 2000 + yy

   date_frac = float(yy) + float( long(dd + 0.5) ) / 365.

;- Calculate the Mean Anomaly of Earth's Orbit

   ANOMEAN = (2.D0 * PI / SOLYR) *365.* (date_frac - (year1 + 3./365.) )

   ANOMEAN = ANOMEAN - 2*PI*long(ANOMEAN/(2*PI))

   E0      = ANOMEAN + ECCENT* SIN(ANOMEAN) + 0.5D0* ECCENT^2 *SIN(2.*ANOMEAN)

   jump1:

   DELTAE = (ANOMEAN - E0 + ECCENT* SIN(E0)) / (1.D0 - ECCENT* COS(E0) )

   if ( ABS(DELTAE/E0) GT 1D-10 AND DELTAE NE 0.0D0) then begin

     e0 = e0 + deltae
     goto,jump1

   endif

;- now the earth-sun distance 

   SOLDIS = AXSMAJ* (1.D0 - ECCENT* COS(E0) )   ; km

;- done

   return,soldis

   end