This code was not written to be user friendly to third parties. I don't have time, or the inclination, to dress it up.
On the other hand it is useable and hence potentially useful.
User beware, that's all.
The code originated for the study described in Rogers and Martin 1979: ApJ, 228, 450.
I later added many subroutines to treat the case of complex refractive indices, and also some checks in some limiting cases (small particles?). Reference to this version could be made to:
Kim and Martin 1995: The Size Distribution of Interstellar Dust Particles as Determined from Polarization: Spheroids. ApJ, 444, 293-305.
The code was originally on cards for an IBM. Various version have since been on a VAX, CRAY XMP (COS), and now several UNIX machines (SUN, SGI, DEC ALPHA). The present version uses the NAG library in one spot (2 routines), but substitutions of library routines should be straightforward.
The required files are available in electronic form in this web archive.
There is a makefile called: Makefileg (also others)
and a driver program gd.f that you can modify to suit, and a gsub.f that you normally would not touch. Also some test input and output files.
To build the executable "g" yourself:
make -f Makefileg
It needs the NAG library. If you don't have that you will get an error: Can't locate file for: -lnag
make -f Makefileg_nonag
It still needs the NAG library, and will produce the error:
Error: Undefined:
f02ajf_
m01daf_
This tells you what routines to substitute in the code, using modules from your favorite libraries.
f02ajf_ : all eigenvalues of a complex matrix
m01daf_ : ranks a vector of real numbers in ascending order
If you profile the program you should see a lot of time spent solving simultaneous equations; subroutine cmateq(a,b,iii,jjj,id) I believe. You might substitute for this equation solver.
-rw-r--r-- 1 pgmartin 300 13895 Jun 17 13:56 sx20o1121out.dat
-rw-r--r-- 1 pgmartin 300 27671 Jun 17 13:56 sx20o1121out.datcheck
-rw-r--r-- 1 pgmartin 300 196883 Jun 17 13:56 x20o1121.out
-rw-r--r-- 1 pgmartin 300 490 Jun 17 11:41 x20o1121.in
Probably the values for the highest x are not converging -- I did not check but you should look at the files to see.
A run closer to what one might want is
-rw-r--r-- 1 pgmartin 300 9959 Apr 24 1993 sx20o1700out.dat
-rw-r--r-- 1 pgmartin 300 478 Apr 24 1993 x20o1700.in
For interstellar applications I think one should think about starting with edge-on oblates, rather than spinning prolates. That would keep the computational burden down.
The x*.out file has lots of verbiage, including diagnostics/warnings. The "s" files summarize the results; their names are built on the fly (supposed to be informative). That is why the *.orig files are there -- the other output files you ftp'd get overwritten when you run the test program as above.
You might want to modify the driver program so that instead of a fixed wavelength and variable size, the setup produces whatever is needed for a fixed size at a set of wavelengths (e.g, for studying the shape of the polarization curve, doing convergence checks). Look carefully at the code so that if you rearrange things the DEPENDENT variables are still computed correctly.
3 39 21 ! ntry,kmaxpup,mtopup
ntry relates to the extra info in the *.datcheck file. The latter two govern the double summation. This is presently set up in a wasteful manner**, since the numbers 39 and 21 are the maximum that can be held by the program with present array sizes, and imply the greatest amount of computing effort. Unlike Mie programs or cylinder programs, you have to know kmax in advance; the program actually works upward to mtop, and so saves a bit if it detects convergence; but mtop has to be set big enough in the first place so that appropriate functions are computed. BUT having kmax as large as 39 is not necessarily needed; it depends on the size parameter and the refractive index. You might as well keep them large until you get some experience.
**To get some idea of this, do some trial computations for several sizes at a range of wavelengths and refractive indices. Look for warning statements in the x*.out output files. Then with this knowledge you might want to tailor kmax and mtop for each category to save production time.