FUNCTION F_MIE,R
;
; Evaluate integrands of extinction and scattering cross-section integrals,
; and transformation matrix integrals
;
	COMMON SIZDIS,DRG,LNSG
	COMMON BH_MIE,DREFRL,K,NANGS
;
;  Compute scattering efficiency
	X = K*R
	BHMIE,X,DREFRL,NANGS,S1,S2,QEXT,QSCA,QBACK,GSCA
;
;  Compute transformation matrices
	S1R = DOUBLE(S1)
	S1I = IMAGINARY(S1)
	S2R = DOUBLE(S2)
	S2I = IMAGINARY(S2)
	NS1 = S1R*S1R + S1I*S1I
	NS2 = S2R*S2R + S2I*S2I
	F11 = 0.5D0*(NS1 + NS2)
	F21 = 0.5D0*(NS1 - NS2)
	F33 = S1R*S2R + S1I*S2I
	F43 = S1I*S2R - S1R*S2I
;
;  Compute particle size density times particle radius, less normalization
	NDISR = EXP(-0.5D0*(ALOG(R/DRG)/LNSG)^2)
;
;  Compute integrands and return
	RETURN,[[QEXT,QSCA]*(R*NDISR),[F11,F21,F33,F43]*(NDISR/R)]
;
	END
FUNCTION ERR_MIE,S,OS,EPS,J
;
; Determine if array of integrals in QSIMPS has converged
;
;  Input:
;     S   - Array of integrals in present iteration
;     OS  - Array of integrals in previous iteration
;     EPS - Array of maximum relative errors for Simpson's Rule
;           integrations
;     J   - Index of present interation
;
;  Output:
;     Status returned as function value:
;       0 - No convergence
;       1 - Successful convergence
;
	COMMON BH_MIE,DREFRL,K,NANGS
;
	IF J LT 3 THEN RETURN,0
	ERR_K = ABS(S(0:1) - OS(0:1)) GT EPS(0)*ABS(OS(0:1))
	IF MAX(ERR_K) EQ 1 THEN RETURN,0
	ERR_P11 = ABS(S(2:2*NANGS) - OS(2:2*NANGS)) GT $
	          EPS(1)*ABS(OS(2:2*NANGS))
	IF MAX(ERR_P11) EQ 1 THEN RETURN,0
	RETURN,1
;
	END
PRO	MIE_LOG_NRM,LAMDA,REFREL,NANG,RG,SG,KEXT,KSCA,P, $
	            ERR_SIMP=EPS,MIN_DIST=EPSN,MAX_RAD=RMAX
;
; Compute the extinction and scattering cross-sections,
; and phase matrices for a log normal size distribution
;
;  Input:
;     LAMDA  - Wavelength (microns)
;     REFREL - Relative index of refraction ( COMPLEX(NR,NI) )
;     NANG   - Number of angles between 0 and 90 degrees, equally spaced
;     RG     - Radius of log normal distribution (microns)
;     SG     - Spread of log normal distribution
;     EPS    - Array of maximum relative errors for Simpson's Rule
;              integrations, passed as a value of keyword ERR_SIMP.
;              The default is [1.E-5,1.E-3].
;     EPSN   - Value of size distribution relative to its maximum,
;              at the integration limits, passed as a value of keyword
;              MIN_DIST.  The default is 1.E-9.
;     RMAX   - Maximum radius for integration, passed as a value of
;              keyword MAX_RAD.  If set, this overrides the upper
;              integration limit determined by EPSN or its default.
;
;  Output:
;     KEXT   - Extinction cross-section (microns^2)
;     KSCA   - Scattering cross-section (microns^2)
;     P      - Phase function:
;                 P(0:2*(NANG-1),0) = P11
;                 P(0:2*(NANG-1),1) = P21
;                 P(0:2*(NANG-1),2) = P33
;                 P(0:2*(NANG-1),3) = P43
;
	COMMON SIZDIS,DRG,LNSG
	COMMON BH_MIE,DREFRL,K,NANGS
;
;  If no input for ERR_SIMP keyword, use default
	IF NOT KEYWORD_SET(EPS) THEN EPS = [1.E-5,1.E-3]
;
;  If no input for MIN_DIST keyword, use default
	IF NOT KEYWORD_SET(EPSN) THEN EPSN = 1.E-9
;
;  Save keyword values as double precision
	DEPS = DOUBLE(EPS)
	DEPSN = DOUBLE(EPSN)
;
;  Save input parameters for use by integration function F_MIE
	DRG = DOUBLE(RG)
	LNSG = ALOG(DOUBLE(SG))
	K = 2.D0*!DPI/DOUBLE(LAMDA)
	DREFRL = DCOMPLEX(REFREL)
	NANGS = FIX(NANG)
;
;  Calculate limits for integrals
	U = EXP(SQRT(-2.0D0*ALOG(DEPSN))*LNSG)
	V = EXP(-LNSG^2)
	R0 = DRG*V
	R1 = R0/U
	R2 = R0*U
	IF KEYWORD_SET(RMAX) THEN R2 = DOUBLE(RMAX)
;
;  Evaluate integrals
	I_MIE = QSIMPS('F_MIE',R1,R2,'ERR_MIE',DEPS)
;
;  Compute extinction and scattering cross-sections
	NORM = SQRT(2.D0*!DPI)*LNSG
	KEXT = !DPI*I_MIE(0)/NORM
	KSCA = !DPI*I_MIE(1)/NORM
;
;  Compute phase matrices
	P = DBLARR(2*NANGS-1,4)
	P(0) = I_MIE(2:8*NANGS-3)*(4.D0*!DPI/(K*K*KSCA*NORM))
;
	RETURN
	END