; PHIG(rsec,bsec,lam)	- scattered light, IDL ver. 5.6
;
; Calculation of the Gaussian contribution to the scattered light.
;
; INPUTS: rsec - distance from disk center in (")
;	  bsec - Gaussian parameter in (")
;	  lam  - working wavelength in A
; CALLING: result = phig(rsec,bsec,lam)
; ROUTINES:  integ_g, clv.pro, beseli.pro (IDL), qromb.pro (IDL)
;
; According to Martinez Pillet 1992, Solar Phys. 140, 207
; 5 Feb 2004, Michal
; -------------------------------------------------------
;
FUNCTION integ_g,rr
;
; Function for PHIG. 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

t0=2.*r*rr/b^2		; argument of Bessel funct. Io
if t0 lt 50. then begin
  t11=(r^2+rr^2)/b^2
  if t11 lt 80. then t1=exp(-1*t11) else t1=0.
  t2=BESELI(t0,0)*rr
endif else begin	; asymptotic approx. for Io
  t11=((r-rr)/b)^2
  if t11 lt 80. then t1=exp(-1*t11) else t1=0.
  t2=(rr*b)/sqrt(4*!pi*r*rr)
endelse

RETURN, cc*t1*t2

END

; -------------------------------------------------------
;
FUNCTION phig,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

; Gaussian normalization coefficient
yy=1./(siRs*b)
a=2*yy^2

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

RETURN,2*!pi*siRs^2*a*resint

END