pro BHCOAT,wavel,rcore,rmant,RRFRL1,RRFRL2,QQEXT,QQSCA,QBACK pi = 4.d0*atan(1.) xx = 2.d0*pi*rcore/wavel yy = 2.d0*pi*rmant/wavel ; pro BHCOAT,XX,YY,RRFRL1,RRFRL2,QQEXT,QQSCA,QBACK ;*********************************************************************** ; ; Subroutine BHCOAT calculates Q_ext, Q_sca, Q_back for coated sphere. ; All bessel functions computed by upward recurrence. ; Input: ; X = 2*PI*RCORE*REFMED/WAVEL ; Y = 2*PI*RMANT*REFMED/WAVEL ; RFREL1 = REFCOR/REFMED ; RFREL2 = REFMAN/REFMED ; where REFCOR = complex refr.index of core) ; REFMAN = complex refr.index of mantle) ; REFMED = real refr.index of medium) ; RCORE = radius of core ; RMANT = radius of mantle ; WAVEL = wavelength of light in ambient medium ; ; Routine BHCOAT is taken from Bohren & Huffman (1983) ; Obtained from C.L.Joseph ; ; History: ; 92/11/24 (BTD) Explicit declaration of all variables ; Some time ago: Original fortran converted to IDL by Mark Hervig ;*********************************************************************** ; ; Arguments: ; REAL XX,YY,QQEXT,QQSCA,QBACK ; COMPLEX RRFRL1,RRFRL2 ; Local variables: ; INTEGER IFLAG,N,NSTOP ; DOUBLE PRECISION ; & CHI0Y,CHI1Y,CHIY,DEL,PSI0Y,PSI1Y,PSIY,QEXT,RN,QSCA,X,Y,YSTOP CHI0Y=0.d0 & CHI1Y=0.d0 & CHIY=0.d0 & DEL=0.d0 & PSI0Y=0.d0 & PSI1Y=0.d0 PSIY=0.d0 & QEXT=0.d0 & RN=0.d0 & QSCA=0.d0 & X=0.d0 & Y=0.d0 & YSTOP=0.d0 ; DOUBLE COMPLEX ; & AMESS1,AMESS2,AMESS3,AMESS4,AN,ANCAP, ; & BN,BNCAP,BRACK, ; & CHI0X2,CHI0Y2,CHI1X2,CHI1Y2,CHIX2,CHIPX2,CHIPY2,CHIY2,CRACK, ; & D0X1,D0X2,D0Y2,D1X1,D1X2,D1Y2,DNBAR,GNBAR,II, ; & REFREL,RFREL1,RFREL2, ; & XBACK,XI0Y,XI1Y,XIY, ; & X1,X2,Y2 AMESS1=dcomplex(0.,0.) AMESS2=dcomplex(0.,0.) AMESS3=dcomplex(0.,0.) AMESS4=dcomplex(0.,0.) AN=dcomplex(0.,0.) ANCAP=dcomplex(0.,0.) BN=dcomplex(0.,0.) BNCAP=dcomplex(0.,0.) BRACK=dcomplex(0.,0.) CHI0X2=dcomplex(0.,0.) CHI0Y2=dcomplex(0.,0.) CHI1X2=dcomplex(0.,0.) CHI1Y2=dcomplex(0.,0.) CHIX2=dcomplex(0.,0.) CHIPX2=dcomplex(0.,0.) CHIPY2=dcomplex(0.,0.) CHIY2=dcomplex(0.,0.) CRACK=dcomplex(0.,0.) D0X1=dcomplex(0.,0.) D0X2=dcomplex(0.,0.) D0Y2=dcomplex(0.,0.) D1X1=dcomplex(0.,0.) D1X2=dcomplex(0.,0.) D1Y2=dcomplex(0.,0.) DNBAR=dcomplex(0.,0.) GNBAR=dcomplex(0.,0.) II=dcomplex(0.,0.) REFREL=dcomplex(0.,0.) ; RFREL1=dcomplex(0.,0.) ; RFREL2=dcomplex(0.,0.) XBACK=dcomplex(0.,0.) & XI0Y=dcomplex(0.,0.) & XI1Y=dcomplex(0.,0.) & XIY=dcomplex(0.,0.) X1=dcomplex(0.,0.) & X2=dcomplex(0.,0.) & Y2=dcomplex(0.,0.) ;- begin routine DEL = 1.D-10 ; was 1d-8 II = dcomplex(0.D0,1.D0) X=XX Y=YY RFREL1=RRFRL1 RFREL2=RRFRL2 ; ----------------------------------------------------------- ; del is the inner sphere convergence criterion ; ----------------------------------------------------------- x1 = rfrel1*x x2 = rfrel2*x y2 = rfrel2*y ystop = y + 4.*y^0.3333 + 2.0 refrel = rfrel2/rfrel1 nstop = ystop ; ----------------------------------------------------------- ; series terminated after nstop terms ; ----------------------------------------------------------- d0x1 = COS(x1)/SIN(x1) d0x2 = COS(x2)/SIN(x2) d0y2 = COS(y2)/SIN(y2) psi0y = cos(y) psi1y = sin(y) chi0y = -sin(y) chi1y = cos(y) xi0y = psi0y-II*chi0y xi1y = psi1y-II*chi1y chi0y2 = -SIN(y2) chi1y2 = COS(y2) chi0x2 = -SIN(x2) chi1x2 = COS(x2) qsca = 0.0 qext = 0.0 xback = dcomplex(0.0,0.0) n = 1 iflag = 0 jump200: rn = n psiy = (2.0*rn-1.)*psi1y/y - psi0y chiy = (2.0*rn-1.)*chi1y/y - chi0y xiy = psiy-II*chiy d1y2 = 1.0/(rn/y2-d0y2) - rn/y2 if (iflag eq 1) then goto, jump999 d1x1 = 1.0/(rn/x1-d0x1) - rn/x1 d1x2 = 1.0/(rn/x2-d0x2) - rn/x2 chix2 = (2.0*rn - 1.0)*chi1x2/x2 - chi0x2 chiy2 = (2.0*rn - 1.0)*chi1y2/y2 - chi0y2 chipx2 = chi1x2 - rn*chix2/x2 chipy2 = chi1y2 - rn*chiy2/y2 ancap = refrel*d1x1 - d1x2 ancap = ancap/(refrel*d1x1*chix2 - chipx2) ancap = ancap/(chix2*d1x2 - chipx2) brack = ancap*(chiy2*d1y2 - chipy2) bncap = refrel*d1x2 - d1x1 bncap = bncap/(refrel*chipx2 - d1x1*chix2) bncap = bncap/(chix2*d1x2 - chipx2) crack = bncap*(chiy2*d1y2 - chipy2) amess1 = brack*chipy2 amess2 = brack*chiy2 amess3 = crack*chipy2 amess4 = crack*chiy2 if (ABS(amess1) gt del*ABS(d1y2)) then goto, jump999 if (ABS(amess2) gt del) then goto, jump999 if (ABS(amess3) gt del*ABS(d1y2)) then goto, jump999 if (ABS(amess4) gt del) then goto, jump999 brack = dcomplex(0.,0.) crack = dcomplex(0.,0.) iflag = 1 jump999: dnbar = d1y2 - brack*chipy2 dnbar = dnbar/(1.0-brack*chiy2) gnbar = d1y2 - crack*chipy2 gnbar = gnbar/(1.0-crack*chiy2) an = (dnbar/rfrel2 + rn/y)*psiy - psi1y an = an/((dnbar/rfrel2+rn/y)*xiy-xi1y) bn = (rfrel2*gnbar + rn/y)*psiy - psi1y bn = bn/((rfrel2*gnbar+rn/y)*xiy-xi1y) qsca = qsca + (2.0*rn+1.0)*(ABS(an)*ABS(an)+ABS(bn)*ABS(bn)) xback = xback + (2.0*rn+1.0)*(-1.)^n*(an-bn) qext = qext + (2.0*rn + 1.0)*(double(an)+double(bn)) psi0y = psi1y psi1y = psiy chi0y = chi1y chi1y = chiy xi1y = psi1y-II*chi1y chi0x2 = chi1x2 chi1x2 = chix2 chi0y2 = chi1y2 chi1y2 = chiy2 d0x1 = d1x1 d0x2 = d1x2 d0y2 = d1y2 n = n + 1 ; if (n-1-nstop) 200,300,300 if (n-1-nstop lt 0) then goto,jump200 if (n-1-nstop ge 0) then goto,jump300 jump300: QQSCA = (2.0/(y*y))*qsca QQEXT = (2.0/(y*y))*qext qback = (ABS(XBACK))^2 qback = (1.0/(y*y))*qback return end