C ARTURO QUIRANTES SIERRA
C Department of Applied Physics, Faculty of Sciences
C University of Granada, 18071 Granada (SPAIN)
C http://www.ugr.es/local/aquiran/codigos.htm
C aquiran@ugr.es
C
C Last update: 19 May 2.003
C
C This is HOMER, a computer code to calculate cross-sections
C of homogeneous spheres, using the Mie theory
C As an added bonus, it can be adapted for magnetic particles
C Plus, you can use it for a range of size parameters
C
C The equations in Homer are taken from the book:
C "Absorption and Scattering of Light by Small Particles"
C by C.F. Bohren and R.D. Huffman, Wiley & Sons, NY, 1983
C (Bohren-Huffman henceforward)
C
C This code is made available on the following conditions:
C
C - You can use and modify it, as long as acknowledgment is given
C to the author (that's me). I'd appreciate a reprint or e-copy
C of any paper where my code has been used.
C
C - You can not sell it or use it in any way to get money out of it
C (save for the research grants you might get thanks to it :)
C
C - I'm making this code available to the best of my knowledge.
C It has been tested as throughly as possible, but bugs may remain
C and overflow errors, particularly for large sizes and/or extreme
C shapes. It work for Windows and Linux platforms. However, it's
C your responsibility to check for its accuracy in the environment
C and size/shape ranges you use.
C
C - Disclaimer. The usual stuff. I claim no responsibility for any
C misuse or damage arising from its use. If this program erases your
C hard disk, makes your computer cry out for Linux, has you fired,
C makes you lose your boy/girlfriend, melts your TV set, spends out
C your lifesaving or forces you to get to decaf coffee, it's your
C problem. For your tranquility, no such side effects has been
C reported ... yet :-)
c
C - And a final message. Don't you think this is the ultimate code.
C I've been improving and speeding it up for a long time, and I think
C there's still room for improvement. So, as you grow confidence,
C I encourage you to strip it, analize it, and improve it at will.
C Improvements, comments, additions: all welcome at aquiran@ugr.es
C
C If it´s your first time with Mie equations, you might like to have
C some benchmark dat to compare with
C All right, here we go:
C
C From Bohren-Huffman, page 114:
C m=1.33+i*10^-8, kr=3
C n An Bn
C 1 5.1631e-1 - i4.9973e-1 7.3767e-1 - i4.3390e-1
C 2 3.4192e-1 - i4.7435e-1 4.0079e-1 - i4.9006e-1
C 3 4.8467e-2 - i2.1475e-1 9.3553e-3 - i9.6269e-2
C 4 1.0347e-3 - i3.2148e-2 6.8810e-5 - i8.2949e-3
C 5 9.0375e-6 - i3.0062e-3 2.8309e-7 - i5.3204e-4
C
C From Wang, van de Hulst, Appl. Opt. 30, 106-117 (1991):
C
C kr m Qext Qabs
C
C 0.0001 1.50 2.3068e-17 0.0
C 1.50 + i0.1 1.9925e-5 1.9925e-5
C 5.21282 1.55 3.10543 0.0
C 1.55 + i0.1 2.86165 1.19740
C 100 1.50 2.0944 0.0
C 1.50 + i0.1 2.0898 0.95769
C 1570.7963 1.342 2.01294 0.0
C 1.342 + i0.1 2.01442 0.93356
C 25000 1.50 2.00235 0.0
C 1.50 + i0.1 2.00232 0.90641
C
C
INTEGER T,N1,NMAX,MSIZE
PARAMETER (MSIZE=50000)
DOUBLE PRECISION KR1,KR2,KR3,KR,QEXT,QSCA,QABS,MUR,PHI(MSIZE),m1,m2
COMPLEX*16 ZETA(MSIZE),AN(MSIZE),BN(MSIZE),M,QQQ
EXTERNAL HANKEL,ABMIE
C The above allows to calculate for size parameters of up to 5000
C You want push the limit up (to allow for larger particles)
C or down (to save memory) by changing the MSIZE value above
C
C WARNING!!!!
C If the particle size is too high, overflow might happen
C Please see remarks on the HANKEL subroutine
C
C M is the particle index of refraction, relative to medium
C For the equations used here, the imaginary part is >=0
M=(1.333D0,0.0d0)
C MUR is the particle's magnetic permeability (relative to vacuum)
C For non-magnetic particles, just set MUR=1
MUR=1.0d0
c This is the file where efficiencies will be stored
OPEN (1,FILE='q.dat')
DO I=1,5
WRITE (*,*)
END DO
WRITE (*,*) ' Welcome to Arturo Quirantes`'
WRITE (*,*)
WRITE (*,*) ' ==================================='
WRITE (*,*) ' =========== H O M E R ============'
WRITE (*,*) ' ==================================='
WRITE (*,*)
WRITE (*,*)
WRITE (*,*) 'This code calculates cross-efficiencies for'
WRITE (*,*) 'sets of spherical, homogeneous (Mie) spheres'
WRITE (*,*)
WRITE (*,*) 'Please enter the (dimensionless) size parameters:'
WRITE (*,*)
10 WRITE (*,*) 'Enter the minimum value of kr'
READ (*,*) kr1
20 WRITE (*,*) 'Enter the maximum value of kr'
WRITE (*,*) '(For a single particle size, make it equal to the minimum value)'
READ (*,*) kr2
C Just in case Homer cannot handle it...
IF((KR+4.0D0*KR**0.33333+2).gt.MSIZE) THEN
WRITE (*,*) 'Sorry, this value is too high'
WRITE (*,*) 'Please enter a value smaller than', int(msize-4.0D0*msize**0.33333+2)
WRITE (*,*) 'Or modify the value of MSIZE and re-complile the code'
WRITE (*,*)
END IF
30 WRITE (*,*) 'Enter the step size for kr'
READ (*,*) kr3
WRITE (*,*) 'Calculation in progress'
C Size parameter loop
DO KR=KR1,KR2,KR3
NMAX=INT(KR+4.0D0*KR**0.33333+2)
C Calculation of Hankel-related Phi, Zeta functions
CALL HANKEL(MSIZE,KR,NMAX,PHI,ZETA)
C Calculation of An, Bn coefficients (n=1 to NMAX)
CALL ABMIE(MSIZE,MUR,KR,M,PHI,ZETA,AN,BN)
QEXT=0.0D0
QSCA=0.0D0
DO T=1,NMAX
QEXT=QEXT+DFLOAT(T+T+1)*REAL(AN(T)+BN(T))
QSCA=QSCA+DFLOAT(T+T+1)*(ABS(AN(T))**2+ABS(BN(T))**2)
END DO
C Qext, Qsca are the extinction and scattering efficiencies
C In order to obtain cross sections, just multiply by pi*r*r
C Where r is the particle radius
QEXT=QEXT*2.0D0/(KR*KR)
QSCA=QSCA*2.0D0/(KR*KR)
QABS=QEXT-QSCA
C And now, let´s save it into disk
WRITE (1,*) KR,real(QEXT),real(QSCA),real(QABS)
WRITE (*,*) KR,QEXT,QSCA
C End of the size parameter loop
END DO
WRITE (*,*)
WRITE (*,*) 'Calculation process is finished'
WRITE (*,*) 'Your efficiency parameters are saved in the "q.dat" file'
WRITE (*,*) 'With the following format: kr,Qext,Qsca,Qabs'
WRITE (*,*)
WRITE (*,*) 'Now HOMER is leaving for a Duff. Hummmmm.... Duff......'
END
C
C
C
C
C Subroutine to calculate the real-argument
c Hankel-related Phi, Zeta functions
C by upward and downward recurrence.
SUBROUTINE HANKEL(MSIZE,KR,NMAX,PHI,ZETA)
INTEGER MSIZE,T,N1,NMAX
DOUBLE PRECISION QA,PHI(MSIZE),YR(MSIZE),JJ(MSIZE),KR,JZ
DOUBLE COMPLEX ZETA(MSIZE)
N1=NMAX
IF(N1.LE.KR) N1=INT(KR)
IF(N1.GT.150) THEN
N1=INT(1.1D0*N1)
ELSE
N1=N1+15
END IF
C First, we´ll calculate JJ by downward recurrence
C W A R N I N G
C It is generally agreed that downward recurrence must be used
C However, for large (kr>10.000) particles, it can lead to
C overflow problems.
C Should this happen, you can still use upward recurrence
C This might be less accurate, so be warned
C
C Downward recurrence
JJ(N1)=0.0D0
JJ(N1-1)=1.0D-35
DO T=N1-2,1,-1
JJ(T)=(DFLOAT(T+T)+3.0D0)*JJ(T+1)/KR-JJ(T+2)
END DO
JZ=3.0D0*JJ(1)/KR-JJ(2)
QA=(SIN(KR)/KR)/JZ
PHI(1)=SIN(KR)/KR
PHI(2)=JJ(1)*QA
C Upward recurrence (next two lines)
PHI(1)=DSIN(kr)/kr
PHI(2)=(phi(1)-cos(kr))/kr
C Next, we´ll calculate the real part of Ha by upward recurrence
YR(1)=-COS(KR)/KR
YR(2)=(YR(1)-SIN(KR))/kr
ZETA(1)=DCMPLX(PHI(1),YR(1))
ZETA(2)=DCMPLX(PHI(2),YR(2))
DO T=2,NMAX
PHI(T+1)=JJ(T)*QA
C Upward recurrence (next line only)
c PHI(T+1)=(DFLOAT(T+T)-1.0D0)*PHI(T)/KR-PHI(T-1)
YR(T+1)=(DFLOAT(T+T)-1.0D0)*YR(T)/KR-YR(T-1)
ZETA(T+1)=DCMPLX(PHI(T+1),YR(T+1))
END DO
CONTINUE
C And back to the main routine
RETURN
END
C
C
C
C
C Subroutine to calculate An, Bn coefficients
SUBROUTINE ABMIE(MSIZE,MUR,KR,M,PHI,ZETA,AN,BN)
INTEGER NT,T,N1,NMAX,RECU,MSIZE
COMPLEX*16 AN(MSIZE),BN(MSIZE),M,D(MSIZE),ZZ,ZETA(MSIZE),ZQ
DOUBLE PRECISION KR,PHI(50000),MUR
NT=INT(KR+4.0D0*KR**0.33333+2)
C Calculation of the logarithmic derivative by downward recurrence
C (Bohren-Huffman, eq. 4.89)
N1=INT(ABS(M*KR))
IF(NT.GT.N1) N1=NT
IF(N1.GT.150) THEN
N1=INT(1.1D0*N1)
ELSE
N1=N1+15
END IF
D(N1)=(0.0D0,0.0D0)
DO T=N1,2,-1
D(T-1)=DFLOAT(T)/(M*KR)-1.0D0/(D(T)+DFLOAT(T)/(M*KR))
END DO
C Calculation of the An, Bn coefficients
C as in Bohren-Huffman, 4.88
DO T=1,NT
ZZ=D(T)*MUR/M+DFLOAT(T)/KR
AN(T)=(ZZ*PHI(T+1)-PHI(T))/(ZZ*ZETA(T+1)-ZETA(T))
ZZ=D(T)*M/MUR+DFLOAT(T)/KR
BN(T)=(ZZ*PHI(T+1)-PHI(T))/(ZZ*ZETA(T+1)-ZETA(T))
END DO
100 CONTINUE
RETURN
END