Wed, 10 Sep 2014
James A Lock <j.lock@csuohio.edu> 

Dear Christopher,

The attachment is a scan of 4 pages of Fortran code to calculate scattering by a finely-
stratified multi-layer sphere.  The basic ideas behind the code are described in Journ. Opt. 
Soc. Am. A25, 2991-3000 (2008), section 3.  Without reading this section first, the code 
probably won't make much sense.  On p.1 of the code, notice the boxed lines around statement 
102 and between 200 and 201.  These lines give the refractive index profile of a generalized 
Luneburg lens.  If your refractive index profile is different, you will have to change theses 
lines. On p.3, AR and AI are the real and imaginary parts of the final partial wave scattering 
amplitudes a-sub-n, and BR and BI are the real and imaginary parts of b-sub-n.  Similarly, 
PT(1,N) and PT(2,N) are the scattering intensities /S1/**2 and /S2/**2.

The program will on rare occasions develop overflow or underflow errors.  But these appear 
to be many fewer than for the progressive iteration approach.  A year ago I gave this code 
to Philip Laven for the project we are working on, and he has reorganized a little bit of 
it and rewritten it in Visual Basic order to decrease the overflow/underflow errors to 
virtually none.  Yesterday he told me that within a month he will be issuing a new version 
of his MIEPLOT code that uses this procedure for a finely-stratified multi-layer sphere.  
So then it will be in the public domain.

If you have any questions about the code, please let me know, and I might be able to answer them.

Sincerely yours,
Jim Lock 



Tue, 30 Sep 2014
James A Lock <j.lock@csuohio.edu> 
Dear Thomas,

Thanks for putting in great effort on my Fortran code for a multilayer sphere.  This code 
uses my parallel iteration method which starts with amplitudes for 256 layers, for example, 
combines them to get amplitudes for 128 composite layers, combines them to get amplitudes 
for 64 doubly-combined layers, etc.  This is based on the successive doubling strategy in 
the FFT algorithm.  The method most of the rest of the world uses is the successive iteration 
method which starts with the amplitudes for the core, adds a layer and updates the combined 
amplitudes, adds another layer and updates the amplitudes again, etc., all the way to the 
sphere surface.  Both strategies are potentially prone to overflow/underflow errors since 
one is calculating the amplitudes for very large partial waves for very small size layers.  
However, since I do many fewer iterations of the combining, I have the gut feeling that my 
parallel iteration method accumulates less round-off error.  I hope it helps you out in 
your project.

About a year ago I gave my code to Philip Laven for the project we are working on together. 
 He made a few organizational changes in the code so that it is less prone to overflows and 
underflows for some of the extreme refractive indices we have been using, and he has been 
using his improved version on the project ever since.  He told me he was going to be putting 
it in the public domain on his web site any week now.  So he may have beaten you to getting 
it in out to the world-at-large.  But certainly feel free to put it on your web site after 
you suitably test its accuracy and stability.  I am glad that I can help.

By the way, thanks for the photos you sent me this morning.  Looking at them brings back 
many good memories.

Take care,
Jim 