function bhmie_matt_optimize, x,refrel,nang,s1,s2,qext,qsca,qback,gsca ;*********************************************************************** ; Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine ; to calculate scattering and absorption by a homogenous isotropic ; sphere. ; Given: ; X = 2*pi*a/lambda ; REFREL = (complex refr. index of sphere)/(real index of medium) ; NANG = number of angles between 0 and 90 degrees ; (will calculate 2*NANG-1 directions from 0 to 180 deg.) ; if called with NANG<2, will set NANG=2 and will compute ; scattering for theta=0,90,180. ; Returns: ; S1(1 - 2*NANG-1) = -i*f_22 (incid. E perp. to scatt. plane, ; scatt. E perp. to scatt. plane) ; S2(1 - 2*NANG-1) = -i*f_11 (incid. E parr. to scatt. plane, ; scatt. E parr. to scatt. plane) ; QEXT = C_ext/pi*a**2 = efficiency factor for extinction ; QSCA = C_sca/pi*a**2 = efficiency factor for scattering ; QBACK = (dC_sca/domega)/pi*a**2 ; = backscattering efficiency [NB: this is (1/4*pi) smaller ; than the "radar backscattering efficiency"; see Bohren & ; Huffman 1983 pp. 120-123] ; GSCA = for scattering ; ; Original program taken from Bohren and Huffman (1983), Appendix A ; Modified by B.T.Draine, Princeton Univ. Obs., 90/10/26 ; in order to compute ; Converted to IDL by P. J. Flatau, Scripps Inst. Oceanography. UCSD 96/11/14 ; 91/05/07 (BTD): Modified to allow NANG=1 ; 91/08/15 (BTD): Corrected error (failure to initialize P) ; 91/08/15 (BTD): Modified to enhance vectorizability. ; 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1 ; 91/08/15 (BTD): Changed definition of QBACK. ; 92/01/08 (BTD): Converted to full double precision and double complex ; eliminated 2 unneed lines of code ; eliminated redundant variables (e.g. APSI,APSI0) ; renamed RN -> EN = double precision N ; Note that DOUBLE COMPLEX and DCMPLX are not part ; of f77 standard, so this version may not be fully ; portable. In event that portable version is ; needed, use src/bhmie_f77.f ; 93/06/01 (BTD): Changed AMAX1 to generic function MAX ; 96/11/14 (PJF): Converted to double precision (again) ; added declaration standarization (strong typing) ; some code "polishing" using Toolpack ; Converted to IDL ; ; 08/11/21 (MMS) Massive Overhaul by Matthew Slivkoff ; More advanced conversion to IDL ; Got rid of about 3 totally pointless loops (inherited from fortran) ; Got rid of about 2 loops within the main for loop ; by utilizing idl's array operations ; Got rid of mxnang. Why does it exist if the user specifies nang??. nmxx=150000 d=dcomplexarr(nmxx) ;*** Obtain double pi: pii = 4.D0*atan(1.D0) dx = x ;flatau I have converted everything to double so drefrl is really not needed drefrl = refrel y = x*drefrl ymod = abs(y) ;*** Series expansion terminated after NSTOP terms ; Logarithmic derivatives calculated from NMX on down xstop = x + 4.D0*x^0.3333D0 + 2.D0 nmx = max(xstop,ymod) + 15 ;flatau IDL fix (convert to integer) nmx=fix(nmx) print, ' nmx ', nmx nstop = xstop ; IF (nmx GT nmxx) THEN begin print, 'error: nmx > nmxx=', nmxx, ' for |m|x=', ymod STOP ENDIF dang=.5D0*pii/double(nang-1) amu=cos(dindgen(nang)*dang) pi0=dblarr(nang) pi1=replicate(1.D0, nang) pi=dblarr(nang) tau=dblarr(nang) s1=dcomplexarr(2*nang-1) s2=dcomplexarr(2*nang-1) ;*** Logarithmic derivative D(J) calculated by downward recurrence ; beginning with initial value (0.,0.) at J=NMX ; ; d(nmx) = dcomplex(0.D0,0.D0) nn = nmx - 1 for n=1, nn do begin ;DO 40 n = 1, nn en = nmx - n + 1 d(nmx-n) = (en/y) - (1.D0/ (d(nmx-n+1)+en/y)) endfor ;40 CONTINUE ; ;*** Riccati-Bessel functions with real argument X ; calculated by upward recurrence ; psi0 = cos(dx) psi1 = sin(dx) chi0 = -sin(dx) chi1 = cos(dx) xi1 = dcomplex(psi1,-chi1) qsca = 0.D0 gsca = 0.D0 p = -1.D0 for n=1, nstop do begin; DO 80 n = 1, nstop en = n fn = (2.D0*en+1.D0)/ (en* (en+1.D0)) ; for given N, PSI = psi_n CHI = chi_n ; PSI1 = psi_{n-1} CHI1 = chi_{n-1} ; PSI0 = psi_{n-2} CHI0 = chi_{n-2} ; Calculate psi_n and chi_n psi = (2.D0*en-1.D0)*psi1/dx - psi0 chi = (2.D0*en-1.D0)*chi1/dx - chi0 xi = dcomplex(psi,-chi) ; ;*** Store previous values of AN and BN for use ; in computation of g= IF (n GT 1) THEN begin an1 = an bn1 = bn ENDIF ; ;*** Compute AN and BN: an = (d(n)/drefrl+en/dx)*psi - psi1 an = an/ ((d(n)/drefrl+en/dx)*xi-xi1) bn = (drefrl*d(n)+en/dx)*psi - psi1 bn = bn/ ((drefrl*d(n)+en/dx)*xi-xi1) ; ;*** Augment sums for Qsca and g= animag=double(an*complex(0.d0,-1.d0)) bnimag=double(bn*complex(0.d0,-1.d0)) qsca = qsca + (2.D0*en+1.D0)* (abs(an)^2+abs(bn)^2) gsca = gsca + ((2.D0*en+1.D0)/ (en* (en+1.D0)))* $ (double(an)*double(bn)+animag*bnimag) IF (n GT 1) THEN begin an1imag=double(an1*complex(0d0,-1.d0)) bn1imag=double(bn1*complex(0d0,-1.d0)) animag=double(an*complex(0d0,-1.d0)) bnimag=double(bn*complex(0d0,-1.d0)) gsca = gsca + ((en-1.D0)* (en+1.D0)/en)*$ (double(an1)*double(an)+an1imag*animag+$ double(bn1)*double(bn)+bn1imag*bnimag) ENDIF ;stop pi=pi1 tau=en*amu*pi1-(en+1.D0)*pi0 s1(0:nang-1)=s1(0:nang-1)+fn*(an*pi+bn*tau) s2(0:nang-1)=s2(0:nang-1)+fn*(an*tau+bn*pi) ;*** Now do angles greater than 90 using PI and TAU from ; angles less than 90. ; P=1 for N=1,3,...; P=-1 for N=2,4,... p = -p temp1=reverse(fn*p*(an*pi-bn*tau)) temp2=reverse(fn*p*(bn*pi-an*tau)) s1(nang:2*nang-2)=s1(nang:2*nang-2)+temp1(1:n_elements(temp1)-1) s2(nang:2*nang-2)=s2(nang:2*nang-2)+temp2(1:n_elements(temp1)-1) ; stop psi0 = psi1 psi1 = psi chi0 = chi1 chi1 = chi xi1 = dcomplex(psi1,-chi1) ; ;*** Compute pi_n for next value of n ; For each angle J, compute pi_n+1 ; from PI = pi_n , PI0 = pi_n-1 pi1=((2.D0*en+1.D0)*amu*pi-(en+1.D0)*pi0)/en pi0 = pi endfor ; 80 CONTINUE ; ;*** Have summed sufficient terms. ; Now compute QSCA,QEXT,QBACK,and GSCA gsca = 2.D0*gsca/qsca qsca = (2.D0/ (dx*dx))*qsca qext = (4.D0/ (dx*dx))*double(s1(0)) qback = (abs(s1(2*nang-2))/dx)^2/pii output={s1:s1,s2:s2,gsca:gsca,qsca:qsca,qext:qext,qback:qback} RETURN,output END