function solim,th,wavein,distnt

;------------------------------------------------------------------------------
;  Routine returns the solar limb darkening curve calculated according
;  the the "Allen" formulation.  The limb darkening curve is a fractional
;  solar intensity (1 at sun center 0 off the sun edge) versus angular
;  displacement from sun center as viewed from the Earth.                    
;                                         
;  input:                                                                    
;    th........ANGLE IN ARCMIN FROM CENTER OF SOLAR DISK AS SEEN FROM EARTH
;    wavein....wavelength,  microns
;    distnt....earth - sun distance,  km,  if zero 1.5e8 km will be used 
;
;  output:
;    solim.....fractional solar intensity (compared to sun center)
;
;  source: Mark Hervig
;
;------------------------------------------------------------------------------

;- setup some variables

   radsun  = 6.9595d5
   pi      = !pi
   convert = pi/10800.
   denom   = 0.07658d0

   WAVNOW = 0.5
   URUN   = 0.95
   VRUN   = -0.20
   CONST  = 0.05
   ACON   = 0.25
   MODi   = 1         ;/.TRUE./

;- The polynomial coefficients

   WAVE = [.32,.35,.37,.40,.45,.5,.6,.8,1.,1.5,2.,3.,5.,10.05] * 1d0
   U    = [.86,.96,1.01,.93,.99,.95,.86,.71,.63,.56,.48,.36,.23,.14] * 1d0
   V    = [.05,-.08,-.15,-.06,-.17,-.20,-.21,-.20,-.19,-.21,-.18,-.13,-.08,-.04] * 1d0

;- The polynomial coefficients, this set has an extrapolation to shorter wavelengths (<0.32)

   WAVE = [.26,.32,.35,.37,.40,.45,.5,.6,.8,1.,1.5,2.,3.,5.,10.05] * 1d0
   U    = [.66,.86,.96,1.01,.93,.99,.95,.86,.71,.63,.56,.48,.36,.23,.14] * 1d0
   V    = [.31,.05,-.08,-.15,-.06,-.17,-.20,-.21,-.20,-.19,-.21,-.18,-.13,-.08,-.04] * 1d0
   
   iallen = n_elements(wave)

;- check the earth-sun distance

   if distnt lt 1e7 then distnt = 1.5d8

;- do the calculation

   SIGMA  = ASIN(1000.d0/DISTNT)/3.
   SIGSQR = 0.5d0/(SIGMA^2)

   MODi   = 1
   RATIO  = DISTNT/RADSUN
   EDGE   = ASIN(1./RATIO)

   if (WAVNOW NE WAVEIN) then begin

     if (WAVEIN GT WAVE(IALLEN-1) OR WAVEIN LT WAVE(0)) then begin
       print,'wavelength out of bounds in solim.pro',wavein
       solim = 0.0
       goto,jump2
     endif

     for i=1,iallen-1 do begin

       if (wave(i) ge wavein) then begin

         URUN   = (U(I)*(WAVEIN-WAVE(I-1))+U(I-1)*(WAVE(I)-WAVEIN))/  (WAVE(I)-WAVE(I-1))
         VRUN   = (V(I)*(WAVEIN-WAVE(I-1))+V(I-1)*(WAVE(I)-WAVEIN))/  (WAVE(I)-WAVE(I-1))
         CONST  = 1 - URUN
         ACON   = CONST - VRUN
         WAVNOW = WAVEIN
         MODi   = 1
         goto,jump1

       endif

     endfor

     print,'wavelength out of bounds in solim.pro',wavein
     solim = 0.0
     goto,jump2
     jump1:

   endif

   if (modi eq 1) then begin

     modi   = 0
     OUTCON = PI*URUN/(EDGE*ACON)

   endif

   THETA = ABS(TH)*CONVERT

   if (THETA LE EDGE) then begin

     SINTHSQ = (RATIO * SIN(THETA))^2
     SOLIM  = CONST + URUN*SQRT(1.-SINTHSQ) - VRUN*SINTHSQ

   endif else begin

     DELT   = THETA - EDGE
     COEF   = DELT*(DELT*SIGSQR + OUTCON)
     SOLIM = 0.0

     if (COEF LE 50.) then SOLIM = ACON * EXP(-COEF)

   endelse

;- done

   jump2:
   return, solim

   end