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