; PHIL(rsec,bsec,lam)	- scattered light, IDL ver. 5.6
;
; Calculation of the Lorentzian contribution to the scattered light.
;
; INPUTS: rsec - distance from disk center in (")
;	  bsec - Lorentzian parameter in (")
;	  lam  - working wavelength in A
; CALLING: result = phil(rsec,bsec,lam)
; ROUTINES:  integ_l, clv.pro, qromb.pro (IDL)
;
; According to Martinez Pillet 1992, Solar Phys. 140, 207
; 5 Feb 2004, Michal
; -------------------------------------------------------
;
FUNCTION integ_l,rr
;
; Function for PHIL. Calculates the integrand.
; rr - distance from the disk center in solar radii Rs
;      (integration variable)
;
common share, Rs,lamb,r,b

; conversion from rr to limb-distance in (") for CLV
ldis=Rs*(1.-rr)
cc=CLV(ldis,lamb)	; center-to-limb variation

r2=r^2
rr2=rr^2
t0=sqrt((b^2+r2+rr2)^2-4*r2*rr2)	; denominator

RETURN, cc*rr/t0

END

; -------------------------------------------------------
;
FUNCTION phil,rsec,bsec,lam
;
common share, Rs,lamb,r,b
;
; unit conversions
Rs = 960. 	; solar radius in arcsec
siRs = 0.004655 ; sin(Rs) = Rs in radians
r=rsec/Rs
b=bsec/Rs
lamb=lam

; Lorentzian normalization coefficient
srs2=siRs^2	; Rs^2 = (sin(Rs))^2
	; k=b^2*srs2+1	is approx. 1 + 1.e-8
	; k2=sqrt(k), approx.=1
zz=alog(4.d0/double(b^2*srs2))	;alog((k2+1)/(k2-1))
a=float(2./(srs2*zz))

; Integration (QROMB from IDL)
resint=QROMB('integ_l',0.,1.)

RETURN,2*!pi*srs2*a*resint

END