
Revised 8-20-96

The program dave is a modification of the excellent Mie Scattering program
by JV Dave written at IBM in 1967.
 
Reference: FORTRAN code listed (and explained) in J. V. Dave, "Subroutines
for Computing the Parameters of the Electromagnetic Radiation Scattered by
Spheres" Report No. 320-3237, May 1968 still available from IBM Palo Alto
Scientific Center, 1530 Page Mill Road Palo Alto, California 94304, USA

The best reference for a quick bootstrap on the mathematics of the Mie
scattering code is  "Electromagnetic Scattering on Spherical
Polydispersions, D. Deirmendjian, American Elsevier Publishing Company, NY,
1969. The instability of the backward Miller recursion in the Dave program
and in many other programs (such as Bohren) is discussed.  

The method of calculating spherical Bessel functions in the Dave program
has been improved to provide stability as well as using memory more
efficiently. The An(y)=-n/y+j{n-1}(y)/j{n}(y) are calculated directly from 
the ratio j{n-1}/j{n} using continued fractions for the highest order 
allowed by the array size and needed for convergence.  If the An array is not 
large enough, another block of An is calculated until the series converges. 
Therefore, ANY size parameter can be calculated without worry of running 
out of memory.  For large size paramenters, the array sizes should be 
increased to minimize the number of times the continued fraction ratio must 
be computed. 

For stability, one should never recur forward in jn(z) or An. For accuracy,
one should recur backwards on j{n-1}/j{n} rather than the full An. Forward
recurrence in y{n}(x), x real, is always stable, but forward recurrence in
y{n}(z),  where z is complex, is NOT stable.  This last point is not well
known.


Reference: W.J. Lentz, 'Generating Bessel functions in Mie scattering 
calculations using continued fractions', Applied Optics, vol 15, 
no 3, pp668-671, March 1976

ERRORS

T1 GREATER THAN NDIM. POSSIBLE ERROR' is an error message that appears
for large size parameters greater than the array size.  It is a WARNING
for people using compilers that limit the integer word length.  The integer 
must be long enough to represent the number of terms, eg 1000000.  A usual 
symptom is this integer will become negative because it is too short and 
the series will NEVER CONVERGE.  Once you know your word length is OK, you
can IGNORE the message as it gives a convenient indication of the status of 
the computations for each cycle of An. 

Most problems with the programs are due to compiler options designed to 
impress buyers with their speed (and inaccuracy). Use the high accuracy 
options in all cases.

TEST PROGRAM
I have included a program, dtable, to calculate spherical bessel functions 
of the first and second kind of complex arguments.  I recently revised the 
program and checked its accuracy against older programs and found it more 
accurate. This can be used to check Mie programs.`  

PROGRAM ASSEMBLY
There are three files to compile. The first ddave is a driver (d)
for the dave subroutine that does all the work.  dave calls 
zratio which calculates the ratios of spherical Bessel functions
of complex argument.  

When you run ddave, it will ask 
INPUT NTHETA,X,RFR,RFI 
ntheta = number of angles to calculate.  This can be zero.
x = size parameter = pi * diameter / lambda (all in the same units)
	Values less than 0.01 may cause overflow unless the convergence 
	criteria is reduced.
rfr = positive real index. This number has been tested from 0.1 to
  about 40. In the infrared, real indices can be less than 1.0 
rfi = positive imaginary index. This number can be much larger than
the real index.  See Peterson and Weinman, Infrared indices of Quartz

Sample run
91,10,1.33,0

This will calculated all the angles from 0 to 180 degrees for the four
stokes parameters as well as the extinction,scattering, and absorption
efficiency for a size parmeter 10 for water.  The size parameter times
pi r squared gives the cross section.  Please see the comments in dave.
The program then asks for stopping angle up to 90 and starting angle.
input 0 for no angles (and a big time saving)
The program calculates the complementary angle automatically. ie
180-theta. You can calculate just 60 degrees if you want.
The following is a sample run done with a Sun.  Compiled versions in 
4 different fortrans have been run.  Note the time call in ddave. 
If you have a get the time subroutine, you may want to add it in ddave.

SAMPLE RUN
Ntheta=number of angles to 2001
x=size parameter
rfr=real index, rfi=(positive)imaginary index
ztheta=starting angle or its complement 180-ztheta
(dtheta + 1)*NTHETA + ztheta must be LE 90 DEGREES
INPUT NTHETA,X,RFR,RFI   
Input angle increment in degrees OR final angle <_ 90 
Input starting angle 
Adjusting dtheta to 
 Angles   0.000000E+00   0.900000E+02   0.100000E+01
 Complementary Angles   0.180000E+03   0.900000E+02   0.100000E+01
 RFR= 0.133000E+01  RFI= 0.000000E+00  X= 0.100000E+02
DTHETA=    1.0000000000000
Working
 QEXT= 0.220655E+01  QSCAT= 0.220655E+01  QABS= 0.000000E+00
   0.00000   0.357396E+04   0.357396E+04   0.357396E+04   0.000000E+00
   1.00000   0.353213E+04   0.352918E+04   0.353066E+04   0.269408E+01
   2.00000   0.340898E+04   0.339742E+04   0.340318E+04   0.104665E+02
   3.00000   0.321139E+04   0.318626E+04   0.319872E+04   0.224212E+02
   4.00000   0.295027E+04   0.290771E+04   0.292868E+04   0.371721E+02
   5.00000   0.263972E+04   0.257731E+04   0.260779E+04   0.529909E+02
   6.00000   0.229609E+04   0.221306E+04   0.225316E+04   0.679934E+02
   7.00000   0.193675E+04   0.183410E+04   0.188301E+04   0.803396E+02
   8.00000   0.157897E+04   0.145938E+04   0.151542E+04   0.884311E+02
   9.00000   0.123875E+04   0.110645E+04   0.116718E+04   0.910813E+02
  10.00000   0.929772E+03   0.790224E+03   0.852671E+03   0.876405E+02
  11.00000   0.662742E+03   0.522242E+03   0.583110E+03   0.780669E+02
  12.00000   0.444791E+03   0.309995E+03   0.365955E+03   0.629302E+02
  13.00000   0.279305E+03   0.156717E+03   0.204676E+03   0.433527E+02
  14.00000   0.166016E+03   0.614742E+02   0.988388E+02   0.208949E+02
  15.00000   0.101359E+03   0.195687E+02   0.444597E+02  -0.260617E+01
  16.00000   0.790570E+02   0.232080E+02   0.346206E+02  -0.252225E+02
  17.00000   0.908707E+02   0.623683E+02   0.602568E+02  -0.451282E+02
  18.00000   0.127428E+03   0.125767E+03   0.111050E+03  -0.607793E+02
  19.00000   0.179064E+03   0.201858E+03   0.176340E+03  -0.710597E+02
  20.00000   0.236602E+03   0.279767E+03   0.245990E+03  -0.753821E+02
  21.00000   0.292009E+03   0.350089E+03   0.311115E+03  -0.737307E+02
  22.00000   0.338881E+03   0.405488E+03   0.364651E+03  -0.666461E+02
  23.00000   0.372740E+03   0.441074E+03   0.401701E+03  -0.551521E+02
  24.00000   0.391130E+03   0.454527E+03   0.419676E+03  -0.406373E+02
  25.00000   0.393524E+03   0.445996E+03   0.418210E+03  -0.247031E+02
  26.00000   0.381072E+03   0.417775E+03   0.398900E+03  -0.899357E+01
  27.00000   0.356238E+03   0.373840E+03   0.364899E+03   0.496894E+01
  28.00000   0.322361E+03   0.319268E+03   0.320415E+03   0.159322E+02
  29.00000   0.283194E+03   0.259623E+03   0.270173E+03   0.230261E+02
  30.00000   0.242472E+03   0.200356E+03   0.218891E+03   0.258318E+02
  31.00000   0.203542E+03   0.146290E+03   0.170824E+03   0.243989E+02
  32.00000   0.169083E+03   0.101206E+03   0.129395E+03   0.192146E+02
  33.00000   0.140930E+03   0.675899E+02   0.969620E+02   0.111277E+02
  34.00000   0.120017E+03   0.465202E+02   0.747107E+02   0.123893E+01
  35.00000   0.106410E+03   0.377174E+02   0.626762E+02  -0.923038E+01
  36.00000   0.994327E+02   0.397183E+02   0.598787E+02  -0.190746E+02
  37.00000   0.978500E+02   0.501546E+02   0.645451E+02  -0.272315E+02
  38.00000   0.100084E+03   0.660940E+02   0.743864E+02  -0.328874E+02
  39.00000   0.104433E+03   0.844042E+02   0.868953E+02  -0.355503E+02
  40.00000   0.109281E+03   0.102100E+03   0.996332E+02  -0.350836E+02
  41.00000   0.113254E+03   0.116641E+03   0.110477E+03  -0.316995E+02
  42.00000   0.115337E+03   0.126149E+03   0.117805E+03  -0.259155E+02
  43.00000   0.114935E+03   0.129538E+03   0.120611E+03  -0.184796E+02
  44.00000   0.111867E+03   0.126543E+03   0.118535E+03  -0.102743E+02
  45.00000   0.106327E+03   0.117663E+03   0.111830E+03  -0.221212E+01
  46.00000   0.987978E+02   0.104021E+03   0.101259E+03   0.486629E+01
  47.00000   0.899495E+02   0.871799E+02   0.879551E+02   0.102807E+02
  48.00000   0.805300E+02   0.689178E+02   0.732505E+02   0.135761E+02
  49.00000   0.712612E+02   0.510124E+02   0.585084E+02   0.145592E+02
  50.00000   0.627553E+02   0.350410E+02   0.449669E+02   0.133035E+02
  51.00000   0.554554E+02   0.222235E+02   0.336141E+02   0.101244E+02
  52.00000   0.496073E+02   0.133214E+02   0.251052E+02   0.552879E+01
  53.00000   0.452585E+02   0.859679E+01   0.197245E+02   0.146066E+00
  54.00000   0.422832E+02   0.783231E+01   0.173941E+02  -0.534997E+01
  55.00000   0.404257E+02   0.104024E+02   0.177209E+02  -0.103197E+02
  56.00000   0.393535E+02   0.153848E+02   0.200755E+02  -0.142274E+02
  57.00000   0.387140E+02   0.216940E+02   0.236883E+02  -0.166951E+02
  58.00000   0.381840E+02   0.282222E+02   0.277519E+02  -0.175348E+02
  59.00000   0.375099E+02   0.339676E+02   0.315169E+02  -0.167573E+02
  60.00000   0.365310E+02   0.381394E+02   0.343707E+02  -0.145576E+02
  61.00000   0.351872E+02   0.402277E+02   0.358921E+02  -0.112807E+02
  62.00000   0.335114E+02   0.400343E+02   0.358785E+02  -0.737134E+01
  63.00000   0.316089E+02   0.376638E+02   0.343440E+02  -0.331671E+01
  64.00000   0.296283E+02   0.334800E+02   0.314926E+02   0.412888E+00
  65.00000   0.277302E+02   0.280356E+02   0.276723E+02   0.341724E+01
  66.00000   0.260564E+02   0.219864E+02   0.233161E+02   0.540784E+01
  67.00000   0.247061E+02   0.160012E+02   0.188804E+02   0.623380E+01
  68.00000   0.237207E+02   0.106788E+02   0.147857E+02   0.589003E+01
  69.00000   0.230789E+02   0.648143E+01   0.113695E+02   0.450774E+01
  70.00000   0.227028E+02   0.369137E+01   0.885310E+01   0.232960E+01
  71.00000   0.224724E+02   0.239388E+01   0.732732E+01  -0.326548E+00
  72.00000   0.222468E+02   0.248613E+01   0.675525E+01  -0.311047E+01
  73.00000   0.218873E+02   0.370874E+01   0.699107E+01  -0.568324E+01
  74.00000   0.212804E+02   0.569398E+01   0.781063E+01  -0.775658E+01
  75.00000   0.203561E+02   0.802315E+01   0.894879E+01  -0.912355E+01
  76.00000   0.190991E+02   0.102848E+02   0.101381E+02  -0.967726E+01
  77.00000   0.175522E+02   0.121266E+02   0.111437E+02  -0.941630E+01
  78.00000   0.158111E+02   0.132947E+02   0.117907E+02  -0.843697E+01
  79.00000   0.140106E+02   0.136559E+02   0.119800E+02  -0.691431E+01
  80.00000   0.123064E+02   0.132027E+02   0.116930E+02  -0.507469E+01
  81.00000   0.108526E+02   0.120407E+02   0.109845E+02  -0.316426E+01
  82.00000   0.978063E+01   0.103617E+02   0.996672E+01  -0.141713E+01
  83.00000   0.917990E+01   0.840861E+01   0.878575E+01  -0.273386E-01
  84.00000   0.908508E+01   0.643601E+01   0.759680E+01   0.871957E+00
  85.00000   0.947020E+01   0.467319E+01   0.653965E+01   0.122026E+01
  86.00000   0.102510E+02   0.329422E+01   0.571885E+01   0.103134E+01
  87.00000   0.112947E+02   0.239882E+01   0.519092E+01   0.385229E+00
  88.00000   0.124368E+02   0.200571E+01   0.495977E+01  -0.587666E+00
  89.00000   0.135012E+02   0.205849E+01   0.498064E+01  -0.172783E+01
  90.00000   0.143223E+02   0.244228E+01   0.517123E+01  -0.287009E+01
  91.00000   0.147657E+02   0.300772E+01   0.542811E+01  -0.386608E+01
  92.00000   0.147447E+02   0.359834E+01   0.564552E+01  -0.460268E+01
  93.00000   0.142312E+02   0.407683E+01   0.573367E+01  -0.501431E+01
  94.00000   0.132584E+02   0.434627E+01   0.563374E+01  -0.508779E+01
  95.00000   0.119163E+02   0.436330E+01   0.532730E+01  -0.485946E+01
  96.00000   0.103405E+02   0.414158E+01   0.483911E+01  -0.440556E+01
  97.00000   0.869461E+01   0.374544E+01   0.423267E+01  -0.382749E+01
  98.00000   0.715037E+01   0.327514E+01   0.359961E+01  -0.323439E+01
  99.00000   0.586624E+01   0.284640E+01   0.304439E+01  -0.272568E+01
 100.00000   0.496853E+01   0.256778E+01   0.266691E+01  -0.237607E+01
 101.00000   0.453660E+01   0.251957E+01   0.254535E+01  -0.222519E+01
 102.00000   0.459425E+01   0.273774E+01   0.272223E+01  -0.227317E+01
 103.00000   0.510813E+01   0.320573E+01   0.319550E+01  -0.248275E+01
 104.00000   0.599323E+01   0.385523E+01   0.391622E+01  -0.278721E+01
 105.00000   0.712479E+01   0.457628E+01   0.479328E+01  -0.310314E+01
 106.00000   0.835468E+01   0.523497E+01   0.570463E+01  -0.334569E+01
 107.00000   0.953035E+01   0.569624E+01   0.651344E+01  -0.344418E+01
 108.00000   0.105138E+02   0.584819E+01   0.708710E+01  -0.335555E+01
 109.00000   0.111981E+02   0.562393E+01   0.731644E+01  -0.307357E+01
 110.00000   0.115193E+02   0.501761E+01   0.713232E+01  -0.263236E+01
 111.00000   0.114629E+02   0.409157E+01   0.651744E+01  -0.210344E+01
 112.00000   0.110633E+02   0.297306E+01   0.551130E+01  -0.158666E+01
 113.00000   0.103969E+02   0.184063E+01   0.420785E+01  -0.119616E+01
 114.00000   0.957017E+01   0.901353E+00   0.274552E+01  -0.104317E+01
 115.00000   0.870419E+01   0.361979E+00   0.129109E+01  -0.121812E+01
 116.00000   0.791730E+01   0.397837E+00   0.192177E-01  -0.177466E+01
 117.00000   0.730901E+01   0.112389E+01  -0.909660E+00  -0.271792E+01
 118.00000   0.694685E+01   0.257232E+01  -0.137029E+01  -0.399897E+01
 119.00000   0.685806E+01   0.468011E+01  -0.128959E+01  -0.551665E+01
 120.00000   0.702714E+01   0.728928E+01  -0.657441E+00  -0.712675E+01
 121.00000   0.739953E+01   0.101601E+02   0.470302E+00  -0.865788E+01
 122.00000   0.789098E+01   0.129968E+02   0.197701E+01  -0.993222E+01
 123.00000   0.840115E+01   0.154822E+02   0.369826E+01  -0.107885E+02
 124.00000   0.882961E+01   0.173182E+02   0.544093E+01  -0.111045E+02
 125.00000   0.909200E+01   0.182658E+02   0.700628E+01  -0.108160E+02
 126.00000   0.913403E+01   0.181804E+02   0.821381E+01  -0.992943E+01
 127.00000   0.894137E+01   0.170367E+02   0.892333E+01  -0.852676E+01
 128.00000   0.854386E+01   0.149399E+02   0.905203E+01  -0.676058E+01
 129.00000   0.801339E+01   0.121210E+02   0.858471E+01  -0.484076E+01
 130.00000   0.745550E+01   0.891513E+01   0.757553E+01  -0.301299E+01
 131.00000   0.699567E+01   0.572688E+01   0.614140E+01  -0.153184E+01
 132.00000   0.676196E+01   0.298408E+01   0.444743E+01  -0.631331E+00
 133.00000   0.686638E+01   0.108702E+01   0.268654E+01  -0.496401E+00
 134.00000   0.738715E+01   0.358580E+00   0.105562E+01  -0.123877E+01
 135.00000   0.835461E+01   0.100157E+01  -0.268651E+00  -0.288021E+01
 136.00000   0.974256E+01   0.306904E+01  -0.115141E+01  -0.534551E+01
 137.00000   0.114665E+02   0.645143E+01  -0.151570E+01  -0.846628E+01
 138.00000   0.133893E+02   0.108830E+02  -0.135123E+01  -0.119954E+02
 139.00000   0.153339E+02   0.159670E+02  -0.714558E+00  -0.156309E+02
 140.00000   0.171018E+02   0.212178E+02   0.279323E+00  -0.190469E+02
 141.00000   0.184951E+02   0.261145E+02   0.147293E+01  -0.219276E+02
 142.00000   0.193399E+02   0.301605E+02   0.268767E+01  -0.240016E+02
 143.00000   0.195077E+02   0.329417E+02   0.374810E+01  -0.250713E+02
 144.00000   0.189326E+02   0.341769E+02   0.450589E+01  -0.250351E+02
 145.00000   0.176220E+02   0.337531E+02   0.486035E+01  -0.238993E+02
 146.00000   0.156594E+02   0.317423E+02   0.477312E+01  -0.217780E+02
 147.00000   0.131979E+02   0.283969E+02   0.427487E+01  -0.188813E+02
 148.00000   0.104467E+02   0.241239E+02   0.346326E+01  -0.154926E+02
 149.00000   0.765025E+01   0.194407E+02   0.249217E+01  -0.119380E+02
 150.00000   0.506293E+01   0.149174E+02   0.155363E+01  -0.855057E+01
 151.00000   0.292247E+01   0.111128E+02   0.854435E+00  -0.563443E+01
 152.00000   0.142475E+01   0.850999E+01   0.590379E+00  -0.343162E+01
 153.00000   0.702750E+00   0.746121E+01   0.921176E+00  -0.209638E+01
 154.00000   0.812547E+00   0.814686E+01   0.194910E+01  -0.167950E+01
 155.00000   0.172749E+01   0.105542E+02   0.370395E+01  -0.212437E+01
 156.00000   0.334158E+01   0.144781E+02   0.613609E+01  -0.327537E+01
 157.00000   0.548169E+01   0.195440E+02   0.911862E+01  -0.489745E+01
 158.00000   0.792678E+01   0.252501E+02   0.124579E+02  -0.670460E+01
 159.00000   0.104325E+02   0.310236E+02   0.159121E+02  -0.839393E+01
 160.00000   0.127581E+02   0.362850E+02   0.192145E+02  -0.968150E+01
 161.00000   0.146919E+02   0.405123E+02   0.220994E+02  -0.103352E+02
 162.00000   0.160739E+02   0.432975E+02   0.243287E+02  -0.102019E+02
 163.00000   0.168112E+02   0.443894E+02   0.257129E+02  -0.922400E+01
 164.00000   0.168859E+02   0.437177E+02   0.261298E+02  -0.744622E+01
 165.00000   0.163540E+02   0.413949E+02   0.255321E+02  -0.500859E+01
 166.00000   0.153361E+02   0.376981E+02   0.239501E+02  -0.212940E+01
 167.00000   0.140001E+02   0.330312E+02   0.214846E+02   0.921328E+00
 168.00000   0.125396E+02   0.278741E+02   0.182941E+02   0.385430E+01
 169.00000   0.111498E+02   0.227242E+02   0.145765E+02   0.639517E+01
 170.00000   0.100037E+02   0.180385E+02   0.105489E+02   0.831698E+01
 171.00000   0.923168E+01   0.141824E+02   0.642769E+01   0.946637E+01
 172.00000   0.890749E+01   0.113918E+02   0.241015E+01   0.978079E+01
 173.00000   0.904139E+01   0.975376E+01  -0.133840E+01   0.929496E+01
 174.00000   0.958138E+01   0.920554E+01  -0.469263E+01   0.813517E+01
 175.00000   0.104224E+02   0.955405E+01  -0.756899E+01   0.650283E+01
 176.00000   0.114224E+02   0.105111E+02  -0.992211E+01   0.464906E+01
 177.00000   0.124222E+02   0.117400E+02  -0.117368E+02   0.284326E+01
 178.00000   0.132679E+02   0.129073E+02  -0.130176E+02   0.133985E+01
 179.00000   0.138317E+02   0.137326E+02  -0.137777E+02   0.346775E+00
 180.00000   0.140295E+02   0.140295E+02  -0.140295E+02   0.000000E+00
 RFR= 0.133000E+01  RFI= 0.000000E+00  X= 0.100000E+02
 QEXT= 0.220655E+01  QSCAT= 0.220655E+01  QABS= 0.000000E+00
 Angles   0.000000E+00   0.900000E+02   0.100000E+01
 Complementary Angles   0.180000E+03   0.900000E+02   0.100000E+01
Ntheta=number of angles to 2001
x=size parameter
rfr=real index, rfi=(positive)imaginary index
ztheta=starting angle or its complement 180-ztheta
(dtheta + 1)*NTHETA + ztheta must be LE 90 DEGREES
INPUT NTHETA,X,RFR,RFI   





The following are the program listings
**************************************************************************
DDAVE.FOR
**************************************************************************

C       REVISED 6/6/86 WJL;4-7-91 WJL;10-27-91 angles WJL
C***************************************************************************
        EXTERNAL DAMIE
C	JX MUST be integer*4 or the iteration index will be lost
        INTEGER*4 JX
        REAL*8 ztheta,ztest,THETD(2001),dtheta
        REAL*8 X,RFR,RFI,QEXT,QSCAT,CTBROS,ELTRMX(4,2001,2)
        INTEGER*4 IMIN,ISEC,IHOUR,NMIN,NSEC,NHOUR
C*************************************************************************
2       JX=1
	print *,'Ntheta=number of angles to 2001'
	print *,'x=size parameter'
	print *,'rfr=real index, rfi=(positive)imaginary index'
	print *,'ztheta=starting angle or its complement 180-ztheta'
	print *,'(dtheta + 1)*NTHETA + ztheta must be LE 90 DEGREES'
        PRINT *, 'INPUT NTHETA,X,RFR,RFI   '
        READ *, JX,X,RFR,RFI
        IF(JX.GT.2001) GO TO 2
50      FORMAT(1X,4HRFR=,e13.6,2X,4HRFI=,e13.6,2X,2HX=,e13.6)
        IF(JX.EQ.0) GO TO 78
 	print *,'Input angle increment in degrees OR final angle <_ 90 '
        READ *, DTHETA
 	print *,'Input starting angle '
        READ *, ZTHETA
c       print *,dtheta,ztheta
        if(int(ztheta).ge.90)  ztheta=180.0-ztheta
        ztest=dtheta*dfloat(jx-1)+ztheta
        if(ztest.gt.90) then
        	print *,'Adjusting dtheta to '
        	dtheta=(dtheta-ztheta)/dfloat(jx-1)
        endif
        ztest=dtheta*dfloat(jx-1)+ztheta
        write(6,60) ztheta,ztest,dtheta
        write(6,61) 180-ztheta,180-ztest,dtheta
60	format(1x,6hAngles,3(2x,e13.6))
61	format(1x,20hComplementary Angles,3(2x,e13.6))
        if(ztheta.ge.ztest) then
        	Print *,' Oops. ztheta inconsistent'
        	goto 2
        endif
	DO 77 J=1,JX
        THETD(J)=dfloat(J-1)*DTHETA+ztheta
77      CONTINUE
78      CONTINUE
        WRITE(6,50) RFR,RFI,X
        IF(JX.NE.0) PRINT *, 'DTHETA=',DTHETA
        print *,'Working'
        CALL DAMIE(X,RFR,RFI,THETD,JX,QEXT,QSCAT,CTBROS,ELTRMX)
        QABS=QEXT-QSCAT
        IF(QABS.LT.1.0E-15) QABS=0.0
        WRITE(6,777) QEXT,QSCAT,QABS
777     FORMAT(1X,5HQEXT=,e13.6,2X,6HQSCAT=,e13.6,2X,5HQABS=,e13.6)
        IF(JX.EQ.0) GO TO 301
        J2X=2*JX+1
        IF(JX.EQ.2001) J2X=J2X-1
334     CONTINUE
34      FORMAT(I8)
        DO 200 J=1,JX
35      FORMAT(3(e13.6))
        WRITE(6,100) THETD(J), (ELTRMX(I,J,1),I=1,4)
C100    FORMAT(F10.5,4(2X,e8.2))
100     FORMAT(F10.5,4(2X,e13.6))
200     CONTINUE
        DO 300 J=1,JX
        JJ=JX+1-J
        THETD(JJ)=180.0-THETD(JJ)
        IF(THETD(JJ).NE.90) WRITE(6,100) THETD(JJ), 
     1  (ELTRMX(I,JJ,2),I=1,4)
300     CONTINUE
301     CONTINUE
        WRITE(6,50) RFR,RFI,X
        WRITE(6,777) QEXT,QSCAT,QABS
	print *,'Seconds= ',SECS
        write(6,60)ztheta,ztest,dtheta
        write(6,61) 180-ztheta,180-ztest,dtheta
	goto 2
        END

		******************dave.for***********************

C       EDITED 4-7-91,6-17-91 W.J. Lentz
c	changed convergence criteria for small x
C*****************************************************************************
C       DAVE SUBROUTINE (IBM 1967) MODIFIED TO USE CONTINUED FRACTIONS WJL 
C*****************************************************************************
C       SUBROUTINE FOR COMPUTING THE PARAMETERS OF THE ELECTROMAGNETIC
C       RADIATION SCATTERED BY A SPHERE.  THIS SUBROUTINE CARRIES OUT ALL
C       COMPUTATIONS IN DOUBLE PRECISION ARITHMETIC.
C       THIS SUBROUTINE COMPUTES THE CAPITAL A FUNCTION BY MAKING USE OF
C       DOWNWARD RECURRENCE RELATIONSHIP AS WELL AS CONTINUED FRACTIONS.
C       X  : SIZE PARAMETER OF THE SPHERE,( 2 * PI * RADIUS OF THE SPHERE)/
C       WAVELENGTH OF THE INCIDENT RADIATION).
C       RF: REFRACTIVE INDEX OF THE MATERIAL OF THE SPHERE. COMPLEX
C       QUANTITY..FORM: (RFR -I *RFI)
C       THETD(J): ANGLE IN DEGREES BETWEEN THE DIRECTIONS OF THE INCIDENT
C       AND THE SCATTERED RADIATION.  THETD(J) IS < OR = 90.0D0
C       IF THETD(J) SHOULD HAPPEN TO BE GREATER THAN 90.0, ENTER WITH
C       SUPPLEMENTARY VALUE; SEE COMMENTS BELOW ON ELTRMX..
C       JX:  TOTAL NUMBER OF THETD FOR WHICH THE COMPUTATIONS ARE
C       REQUIRED.  JX SHOULD NOT EXCEED NDIM UNLESS THE DIMENSIONS
C       STATEMENTS ARE APPROPRIATELY MODIFIED.
C       MAIN PROGRAM SHOULD ALSO HAVE REAL * 8 THETD(100),ELTRMX(4,NDIM,2).
C       THE DEFINITIONS FOR THE FOLLOWING SYMBOLS CAN BE FOUND IN 'LIGHT
C       SCATTERING BY SMALL PARTICLES,H.C. VAN DE HULST, JOHN WILEY & SONS
C       INC.,NEW YORK, 1957'.
C       QSCAT: EFFICIENCY FACTOR FOR SCATTERING, VAN DE HULST, P.14 & 127.
C       CTBRQS: AVERAGE(COSINE THETA) * QSCAT, VAN DE HYULST, P.128.
C       ELTRMX(I,J,K):ELEMENTS OF THE TRANSFORMATION MATRIX F,VAN DE
C       HULST,P.34,45 & 125,I=1: ELEMENT M SUB 2..I=2: ELEMENT M SUB 1..
C       I=3: ELEMENT S SUB 21.. I=4: ELEMENT D SUB 21..
C       ELTRMX(I,J,1)REPRESENTS THE ITH ELEMENT OF THE MATRIX FOR
C       THE ANGLE THETD(J).. EKTRNX(I,J,2) REPRESENTS THE ITH ELEMENT
C       OF THE MATRIX  FOR THE ANGLE 91.0 - THETD(J)
C       ELTRMX(2,I,J) IS THE VERTICAL COMPONENT SCATTERING I1
C       ELTRMX(1,I,J) IS THE HORIZONTAL COMPONENT SCATTERING I2
C*****************************************************************************
        SUBROUTINE DAMIE(X,RFR,RFI,THETD,JX,QEXT,QSCAT,CTBROS,ELTRMX)
        EXTERNAL ZRATIO
C	NDIM MUST match A(DIM) below in zratio.for
	COMPLEX*16 A(5000)
        COMPLEX*16 ZRATIO
        REAL*8 X,RX,RFR,RFI,QEXT,QSCAT,T(5),TA(4),TB(2),TC(2)
        REAL*8 TD(2),TE(2),CTBROS,PII
        REAL*8 ELTRMX(4,2001,2),PI(4,2001),TAU(4,2001),CSTHT(2001)
        REAL*8 SI2THT(2001),THETD(2001),MAX,XN
        COMPLEX*16 RF,RRF,RRFX,WM1,FNA,FNB,TC1
        COMPLEX*16 FNAP,FNBP,TC2,WFN(2),ACAPN
        COMMON TA,TB,TC,TD,TE
        EQUIVALENCE (WFN(1),TA(1)), (FNA,TB(1)),(FNB,TC(1))
        EQUIVALENCE (FNAP,TD(1)),(FNBP,TE(1))
C*****************************************************************************
C	CONVERGENCE CRITERIA
	ZCONV=1.0D-10
        MAX=0
C*************************************************************************
C	NDIM MUST match A(DIM) below in zratio.for
	NDIM=5000
C*************************************************************************
        PII=3.14159265358979D0
        pii=datan(1.0d0)*4.0d0
        IF(JX.LE.2001) GO TO 20
        PRINT *, 'PLEASE READ THE COMMENTS; JX MUST BE LESS THAN 2001'
        STOP
20      RF=CMPLX(RFR,-RFI)
        T(1)=(X*(RFR+RFI))+9
        IF(T(1).GT.NDIM) PRINT *, 'T1 GREATER THAN NDIM. POSSIBLE
     1  ERROR',T(1)
        MAX=T(1)
        RRF=1.0D0/RF
        RX=1.0D0/X
        RRFX=RRF*RX
        DO 30 J=1,JX
        IF (THETD(J).LT.0.0) THETD(J)=ABS(THETD(J))
        IF (THETD(J).GT.0.0) GO TO 24
        CSTHT(J)=1.0D0
        SI2THT(J)=0.0D0
        GO TO 30
24      IF (THETD(J).GE.90.0) GO TO 25
        T(1)=PII*THETD(J)/180.0D0
        CSTHT(J)=COS(T(1))
        SI2THT(J)=1.0D0-CSTHT(J)**2
        GO TO 30
25      IF (THETD(J).NE.90.0) GO TO 28
        CSTHT(J)=0.0D0
        SI2THT(J)=1.0D0
        GO TO 30
28      IF(THETD(J).NE.180.0) GO TO 29
        CSTHT(J)=-1.
        SI2THT(J)=0.
        GO TO 30
29      T(1)=PII*(180.0-THETD(J))/180.0D0
        CSTHT(J)=-COS(T(1))
        SI2THT(J)=1.0D0-CSTHT(J)*CSTHT(J)
30      CONTINUE
        DO 35 J=1,JX
        PI(1,J)=0.0D0
        PI(2,J)=1.0D0
        TAU(1,J)=0.0D0
        TAU(2,J)=CSTHT(J)
35      CONTINUE
        T(1)=COS(X)
        T(2)=SIN(X)
        WM1=CMPLX(T(1),-T(2))
        WFN(1)=CMPLX(T(2),T(1))
        WFN(2)=RX*WFN(1)-WM1
        T(1)=RFI*X
        XN=1.
        ACAPN= ZRATIO(XN,X,RF,MAX,A,NDELTA)
        TC1=ACAPN*RRF+RX
        TC2=ACAPN*RF+RX
        FNA=(TC1*TA(3)-TA(1))/(TC1*WFN(2)-WFN(1))
        FNB=(TC2*TA(3)-TA(1))/(TC2*WFN(2)-WFN(1))
        FNAP=FNA
        FNBP=FNB
        T(1)=1.5
        TB(1)=T(1)*TB(1)
        TB(2)=T(1)*TB(2)
        TC(1)=T(1)*TC(1)
        TC(2)=T(1)*TC(2)
        DO 60 J=1,JX
        TAU2J=TAU(2,J)
        ELTRMX(1,J,1)=TB(1)+TC(1)*TAU2J
        ELTRMX(2,J,1)=TB(2)+TC(2)*TAU2J
        ELTRMX(3,J,1)=TC(1)+TB(1)*TAU2J
        ELTRMX(4,J,1)=TC(2)+TB(2)*TAU2J
        ELTRMX(1,J,2)=TB(1)-TC(1)*TAU2J
        ELTRMX(2,J,2)=TB(2)-TC(2)*TAU2J
        ELTRMX(3,J,2)=TC(1)-TB(1)*TAU2J
        ELTRMX(4,J,2)=TC(2)-TB(2)*TAU2J
60      CONTINUE
        QEXT=2.0D0*(TB(1)+TC(1))
        QSCAT=(TB(1)**2+TB(2)**2+TC(1)**2+TC(2)**2)/0.75
        CTBROS=0.0D0
        XN=2
        IFLAG=-1
65      T(1)=2.*XN-1.
        T(2)=XN-1.
        T(3)=2.*XN+1.
        DO 70 J=1,JX
        PI1J=PI(1,J)
        PI2J=PI(2,J)
        CSTHTJ=CSTHT(J)
        PI(3,J)=(T(1)*PI2J*CSTHTJ-XN*PI1J)/T(2)
        TAU(3,J)=CSTHTJ*(PI(3,J)-PI1J)-T(1)*SI2THT(J)*PI2J+TAU(1,J)
70      CONTINUE
        WM1=WFN(1)
        WFN(1)=WFN(2)
        WFN(2)=T(1)*RX*WFN(1)-WM1
        ACAPN= ZRATIO(XN,X,RF,MAX,A,NDELTA)
        TC1=ACAPN*RRF+XN*RX
        TC2=ACAPN*RF+XN*RX
        FNA=(TC1*TA(3)-TA(1))/(TC1*WFN(2)-WFN(1))
        FNB=(TC2*TA(3)-TA(1))/(TC2*WFN(2)-WFN(1))
        T(5)=XN
        T(4)=T(1)/(T(5)*T(2))
        T(2)=(T(2)*(T(5)+1.0D0))/T(5)
        CTBROS=CTBROS+T(2)*(TD(1)*TB(1)+TD(2)*TB(2)+TE(1)*
     1 TC(1)+TE(2)*TC(2))+T(4)*(TD(1)*TE(1)+TD(2)*TE(2))
        QEXT=QEXT+T(3)*(TB(1)+TC(1))
        T(4)=TB(1)**2+TB(2)**2+TC(1)**2+TC(2)**2
        QSCAT=QSCAT+T(3)*T(4)
        T(2)=XN*(XN+1.)
        T(1)=T(3)/T(2)
        DO 80 J=1,JX
        PI3J=PI(3,J)
        TAU3J=TAU(3,J)
        ELTRMX(1,J,1)=ELTRMX(1,J,1)+T(1)*(TB(1)*PI3J+TC(1)*TAU3J)
        ELTRMX(2,J,1)=ELTRMX(2,J,1)+T(1)*(TB(2)*PI3J+TC(2)*TAU3J)
        ELTRMX(3,J,1)=ELTRMX(3,J,1)+T(1)*(TC(1)*PI3J+TB(1)*TAU3J)
        ELTRMX(4,J,1)=ELTRMX(4,J,1)+T(1)*(TC(2)*PI3J+TB(2)*TAU3J)
        IF (IFLAG.EQ.-1) GO TO 75
        ELTRMX(1,J,2)=ELTRMX(1,J,2)+T(1)*(TB(1)*PI3J-TC(1)*TAU3J)
        ELTRMX(2,J,2)=ELTRMX(2,J,2)+T(1)*(TB(2)*PI3J-TC(2)*TAU3J)
        ELTRMX(3,J,2)=ELTRMX(3,J,2)+T(1)*(TC(1)*PI3J-TB(1)*TAU3J)
        ELTRMX(4,J,2)=ELTRMX(4,J,2)+T(1)*(TC(2)*PI3J-TB(2)*TAU3J)
        GO TO 80
75      ELTRMX(1,J,2)=ELTRMX(1,J,2)+T(1)*(-TB(1)*PI3J+TC(1)*TAU3J)
        ELTRMX(2,J,2)=ELTRMX(2,J,2)+T(1)*(-TB(2)*PI3J+TC(2)*TAU3J)
        ELTRMX(3,J,2)=ELTRMX(3,J,2)+T(1)*(-TC(1)*PI3J+TB(1)*TAU3J)
        ELTRMX(4,J,2)=ELTRMX(4,J,2)+T(1)*(-TC(2)*PI3J+TB(2)*TAU3J)
80      CONTINUE
C	The first case says the series must be small and ten terms
C	should be used to eliminate false convergence
C	The second case is for extremely small size parameters to
C	prevent underflow by forcing ten or too many terms
        IF(((T(4).LT.ZCONV).AND.XN.GT.8.0).OR.(T(4).LT.1.0D-30))
     X	GO TO 100
        XN=XN+1.
        IFLAG=-IFLAG
        DO 90 J=1,JX
        PI(1,J)=PI(2,J)
        PI(2,J)=PI(3,J)
        TAU(1,J)=TAU(2,J)
        TAU(2,J)=TAU(3,J)
90      CONTINUE
        FNAP=FNA
        FNBP=FNB
        GO TO 65
100	DO 120 J=1,JX
        DO 120 K=1,2
        DO 115 I=1,4
        T(I)=ELTRMX(I,J,K)
115     CONTINUE
        ELTRMX(1,J,K)=T(3)**2+T(4)**2
        ELTRMX(2,J,K)=T(1)**2+T(2)**2
        ELTRMX(3,J,K)=T(1)*T(3)+T(2)*T(4)
        ELTRMX(4,J,K)=T(2)*T(3)-T(4)*T(1)
C       ELTRMX(2,I,J) IS THE VERTICAL COMPONENT SCATTERING I1
C       ELTRMX(1,I,J) IS THE HORIZONTAL COMPONENT SCATTERING I2
120     CONTINUE
        T(1)=2.0D0*RX**2
        QEXT=QEXT*T(1)
        QSCAT=QSCAT*T(1)
        CTBROS=2.0D0*CTBROS*T(1)
        RETURN
        END
		**************zratio.for**************************
C       REVISED 6/6/86 WJL;12-28-88 WJL;4-7-91 WJL; 
C***********************************************************************
C       PROGRAM RATIO CALCULATES THE AN VALUES BASED ON THE RATIOS OF
C       SPHERICAL BESSEL FUNCTIONS OF COMPLEX ARGUMENT
C       SETUP FOR DOUBLE PRECISION ONLY-PROGRAM CAN RUN IN SINGLE PRECISION
C       NDIM MUST!!!! equal A(NDIM) in zratio and dave A(NDIM) must match
C************************************************************************
        COMPLEX*16 FUNCTION ZRATIO(XN,X,RF,MAX,A,NDELTA)
        COMPLEX*16 RF,Y,ZANP,ZNUM,ZDEN,ZPDT,ZRPDT,ZAN
        REAL*8 X,V,XN,MAX,MIN,J,ZREAL,ZIMAG
        INTEGER*4 JN,JJ,NDELTA,NDIM
        LOGICAL NUMFLAG,DENFLAG
C************************************************************************
C	NDIM MUST MATCH DIMENSION OF A(NDIM)
	COMPLEX*16 A(5000)
C************************************************************************
	NDIM=5000
C       PRINT *, 'MAX=',MAX
        IF(XN.LT.1.5.AND.MAX.LT.NDIM) NDELTA=int4(MAX)
        IF(XN.LT.1.5.AND.MAX.GE.NDIM) NDELTA=NDIM
        IF(XN.LT.1.5) MAX=0.
        IF(XN.LT.(MAX+1)) GO TO 2
C       PRINT *, N,MAX,'   AFTER GO TO 2'
        MAX=MAX+DFLOAT(NDELTA)
        MIN=MAX+1.-DFLOAT(NDELTA)
C       PRINT *, MAX,MIN
C*****************************************************************************
C       CALCULATE RATIO OF BESSEL FNS OF CONSECUTIVE ORDER
C       PRINT *, 'CALCULATING NEW CONTINUED FRACTION RATIO',MAX,NDIM
	NUMFLAG=.FALSE.
	DENFLAG=.FALSE.
C	NUMERATOR AND DENOMINATOR FLAGS FALSE ALLOW NORMAL RECURSION
        V=MAX+0.5
        Y=RF*X
        ZANP=2.0D0/Y
        ZNUM=ZANP*V
        ZPDT=ZNUM
        V=V+1.0
C       PRINT *, 'V=',V
        ZDEN=ZANP*V
        ZNUM=ZDEN-1.0D0/ZNUM
444     ZRPDT=ZNUM/ZDEN
        ZPDT=ZRPDT*ZPDT
        ZREAL=DREAL(ZRPDT)
        ZIMAG=DIMAG(ZRPDT)
        IF (DABS(ZREAL-1.0D0).LT.1.0D-12.AND.DABS(ZIMAG).LT.1.0D-12) 
     1  GO TO 555
C       PRINT *, 'TEST ZRPDT CONVERGENCE'
779     V=V+1.0
C       PRINT *, 'V=',V
        ZAN=ZANP*V
C	If the numerator has been zero at some time, then numflag
C	will have been set true
        IF(NUMFLAG) THEN
        PRINT *,'NUMFLAG SET'
          IF(ZNUM.EQ.1.0D0) THEN
            ZNUM=2.0D0
	  ELSE
	    ZNUM=ZAN
            NUMFLAG=.FALSE.
          ENDIF
        ELSE
          ZNUM=ZAN-1.0D0/ZNUM
          IF(ZNUM.EQ.0) THEN
            NUMFLAG=.TRUE.
            ZNUM=1.0D0
          ENDIF
        ENDIF
        IF(DENFLAG) THEN
        PRINT *,'DENFLAG SET'
          IF(ZDEN.EQ.1.0D0) THEN
            ZDEN=2.0D0
          ELSE
            ZDEN=ZAN
            DENFLAG=.FALSE.
          ENDIF
        ELSE
          ZDEN=ZAN-1.0D0/ZDEN
          IF(ZDEN.EQ.0) THEN
            DENFLAG=.TRUE.
            ZDEN=1.0D0
          ENDIF
        ENDIF
        GO TO 444
555     CONTINUE
d	PRINT *,'CF CONVERGED FOR V=',V,' MAX=',MAX
C**********************************************************************
C	Downward recursion of ratios to form An from cf ratio
        J=MAX
1       JJ=int4(J-MAX+DFLOAT(NDELTA)+.001)
        IF(JJ.LT.1.OR.JJ.GT.NDIM) PRINT *, 'JJ ERROR',JJ
        A(JJ)=-(J)/Y+ZPDT
        J=J-1
        IF(J.LT.MIN) GO TO 2
        ZPDT=(2.*J+1.)/Y-1.0/ZPDT
        GO TO 1
2       CONTINUE
        JN=int4(XN-MAX+DFLOAT(NDELTA)+.001)
        J=DFLOAT(JN)
        IF(JN.LT.1.OR.JN.GT.NDIM) PRINT *, 'J ERROR',JN
        ZRATIO=A(JN)
        RETURN
        END

		**************************************************
		The value of 5000 is small to prevent memory conflicts. It
		should be changed if large size parametes greater than 1000
		are run.  Be sure to change it in dave and in zratio at all
		places.  If you do this, make sure int4 works in your compiler
		or the integer will wrap around at 32768 causing an infinite
		loop. 



