Updated April 14, 2005 by Bill Coles, bcoles@ucsd.edu, 858 534 2703
These programs simulate small angle forward scattering such as occurs in the neutral atmosphere at optical frequencies, and in the ionosphere, the solar wind, and the interstellar plasma at radio frequencies. Related scattering occurs in acoustic propagation in the ocean but it is usually severely modified by the stratified medium. The primary simulation code originated with Joe Filice at UCSD and was used as part of his PhD thesis. It was extended by Rod Frehlich and M. Yadlowski at Boulder, and subsequently refined further by Rod Frehlich. The present code was modified heavily from Rod's last version by Bill Coles and Barney Rickett at UCSD.
There are two separate simulation programs, one for spherical waves and one for plane waves. The spherical wave program is required for a pulsar embedded in an extended scattering region. The plane wave program is required for an extra- galactic source scintillating in the interstellar plasma. It can also be used for a pulsar problem where the scattering comes from a thin region, although the effective distance and velocity must be mapped appropriately.
The canonical problem is a "thin phase screen." In this case the incident field is first phase modulated by the screen. Then that screen is propagated forward by fourier transforming the electric field into an angular spectrum of plane waves. Each plane wave is then propagated to the plane of interest with its propagation constant. Then the angular spectrum is converted back into a spatial electric field by inverse fourier transformation. An extended scattering region is modelled by a series of thin screens. This approach is discussed in detail in by Coles, Filice, Frehlich and Yadlowsky, Applied Optics 34, 2089, 1995.
The plane wave program is a little more efficient than the spherical wave program because each screen has the same statistics. The spherical wave is modelled in a spherical coordinate system, which requires that each screen be scaled separately. This takes twice as long to generate the screens but this is not usually the main computation time factor.
Both the plane wave and the spherical wave can be "beamed", that is the intensity of the incident field can be reduced to zero smoothly at the edge of the analysis window. This reduces effects of the discrete fourier transform, which implies that the window is periodic. These effects are not important if the phase screen is periodic, however if a phase gradient is introduced the phase screen will not be periodic. In this case a "beam" must be used or severe errors will occur. In any case use of a phase gradient requires careful checking that the screen sampling and window size are adequate.
In many cases the observations are time series of a pattern convected past the observer with some velocity. This can often be modelled by simulating a spatial pattern and taking a scaled cut through that pattern. However velocity shear in the scattering medium cannot be modelled with either of these codes.
The screen phase fluctuations are modelled by creating a two-dimensional gaussian random process and weighting it by the square root of the spatial spectrum of phase, then inverse transforming it to obtain a phase realization. This provides a periodic phase screen. To this the user can add a phase gradient (with great care). The spectral model is a power law with a gaussian cutoff at an "inner scale". It may be anisotropic as parameterized by an axial ratio and orientation angle. In this description the Kolmogorov exponent is 5/3.
The key to an accurate simulation is getting the sampling of the screen fine enough and the window size large enough, without having so many points that the simulation is very slow. There are two primary parameters that set the basic scales. First is the strength of scattering Mb^2. This is the intensity variance under the weak scattering approximation. It is linearly proportional to the phase variance, but it also depends on the wavelength and the distance from the screen. The dividing line between weak and strong scattering is Mb^2 = 1. Second is the Fresnel scale Rf = sqrt(lamdba Z/2pi) where lambda is the wavelength and Z the distance.
In weak scattering the spatial scale of intensity is Rf. You need to sample at least 10 samples per Rf, and you need a window at least 10 Rf. So Rf should be the geometric mean of the sample interval and the window size. We normally set Rf = 1, and set dx = Rf/sqrt(n) as a first guess. This can be adjusted a bit depending on the inner scale, axial ratio, and spectral exponent.
In strong scattering the spatial scales of intensity bifurcate into So and Sr. The geometric mean of So and Sr is Rf, so we continue to keep Rf at the geometric mean of dx and ndx. However So gets smaller and Sr gets larger so n must increase as the scattering gets stronger. When the spectrum is anisotropic Sr gets larger even faster and one must allow a larger window.
The only reliable way to get the scales right is to run a couple of tests first.
1. Do a weak, single-frequency, beamed simulation with the spectral parameters you want and a modest window size, say 512 square. Examine this run first with iplot. You should see a square beam in the middle of the window. If the beam is small or misshapen you probably need a larger window, i.e. keep Rf=1, but increase dx. When iplot looks reasonable for the beamed simulation, you can do an unbeamed simulation. Check iplot, then pdfm, then corr2dn. Corr2dn should give you a clean correlation with a spatial scale roughly Rf. Keep track of the actual M^2, because the Mb^2 is just a rough estimate for ideal spectra. You can use the actual M^2/Mb^2 ratio to correct Mb^2 in strong scattering.
2. Then increase the strength of scattering to the Mb^2 you have in mind and increase n to what seems reasonable. Again check iplot, pdfm, and corr2dn. In this case you want to be sure that the correlation has decreased to zero by half the window size. Increase n by a factor of two and see if the key parameters change significantly. If not then your window size was probably OK and you can run your dynamic spectrum with this scaling.
Simulation of the intensity time series must be done carefully. The simulations provide the electric field in a window around the observer. If the observer is moving at a velocity V with respect to the pattern, then we can simply take a cut through the spatial pattern where t = x dot V/V^2. This is often the case in plane wave simulations. Then the apparent velocity is just Vapp = Vobs - Vmedium. The velocity of the source does not matter since it is at infinite distance. It is important that the medium velocity is constant, if there is velocity shear within the medium then a different type of simulation will be required.
The spherical wave problem is more complex because the velocity of the source is important. The actual velocity seen by each slab of medium is the velocity of the line of sight from the source to the observer with respect to the medium. If the source is fixed with respect to the medium then we see no velocity shear, just the velocity of the observer. If the observer is fixed with respect to the medium then we see just the velocity of the source. We can prove this using reciprocity. However the effective location of the scattering medium along the line of sight must be inverted because the reciprocal problem is to put a source at the observer and calculate the field in the plane the source is traversing. Then the source moves through this field with its own velocity and there is no shear. If neither source nor observer is fixed with respect to the medium then a different simulation is required.
It is obviously interesting to consider a "screen" with a phase gradient. Such deterministic components have been invoked to explain extreme scattering events in the interstellar plasma and peculiar features in dynamic spectra. However the simulation is a non-trivial extension of this code. We have found that introducing a phase gradient makes the phase non-periodic and caused discontinuities at the edge of the window. These cause high-angle scattering and this interfers with the simulation. We introduced a beam to eliminate the discontinuities. The beamed simulation requires a much larger window size. This is because the simulation is unusually efficient for homogeneous media. Waves scattered at high angles cannot leave the simulation box because it is periodic. However with a beam highly scattered waves can leave the beam. Thus the beam width B must be > 2Z Thetamax. This reduces to B > 2 pi Rf^2/dx. If dx = 0.05 and Rf=1, we need B > 120 Rf = 2400 samples. So a window with n=4096 is barely adequate. By comparison the case with dx=0.05 can be simulated without a beam using a window n=256! The basic minimum test is to run the same simulation with and without a beam but with no phase gradient. The secondary spectra must be the same.
We create dynamic spectra by running a complete simulation at each frequency and writing that file to disc. This can be a very large file (30 GB for a 4096x4096x256 simulation). This file is then read and a dynamic spectrum file created in a second step. The dynamic spectrum file is very much smaller (by a factor of 8192 in this example). If you are limited in storage you could combine the two steps, but 30 GB is not too bad if you run the code from a script and remove the large file immediately after creating the dynamic spectrum file. To choose dlam and nf efficiently one should try to get the parabolic arc to run through the corner of the secondary spectrum. This requires that dlam/nf = (2/pi)*(dx/Rf)^2.
(a) simpw filename does the plane wave simulation and creates "filename.out" (b) simsw filename does the spherical wave simulation and creates "filename.out" (c) iplot filename plots the intensity in the image plane (d) pdfm filename plots the probability density of intensity (e) bplot filename plots the apparent brightness distribution as viewed (f) corr2dm filename plots a 2-d correlation of intensity and two orthogonal cuts through it. It also plots cuts of the power spectra (g) getdyn filename reads the output file and creates a dynamic spectrum file (assuming that nf > 1). The header is duplicated on the dynamic spectrum file "filename.spc" (h) plotdyn filename plots the dynamic spectrum file and the secondary spectrum
The first line is a comment that will not be retained in the subsequent files. The file is ascii line oriented with the parameters first. There are one or two parameters spaced by blanks. Everything after the last parameter is a comment. The simulation program reads the first 31 characters in each line and puts them in the output file with a new line (32 bytes per line). Thus you can read the first few lines of the output file with any ascii editor or display. The first two lines of the output file are the name of the simulation program and the date/time. After that are the lines from the input parameter file. These files must be called "filename.in", this one is called "weak.in".
comment line 0.5 MB2 BORN VARIANCE 0.1 DX SAMPLE SPACING 1.0 RF FRESNEL SCALE 1.66667 ALPHA STRUCTURE FUNCTION EXPONENT 1.0 0.0 AXIRAT ANGLE 0.01 L0 INNER SCALE 256 N NUMBER OF SAMPLES ACROSS WINDOW 1 NS NUMBER OF THIN SCREENS ALONG PATH 10 NREAL NUMBER OF REALIZATIONS 999 IDUM SEED FOR RANDOM NUMBER GENERATOR swdsp SPECT function for plasma 1 NF number of frequencies 1.0 DLAM wavelength range (max-min)/mean 1.0 TOVERZ thickness/meandist 0.0 phgrad phase gradient in radians/So n beam beamed y/n y lamscale lambda scaling y/n
The code is vanilla fortran 77 and compiles with any fortran compiler we have tried. We use g77, though, so if you have trouble with your compiler try it with g77 before you waste a lot of time. The plotting programs rely on pgplot which is documented in http://www.astro.caltech.edu/~tjp/pgplot/ and is available for macintosh osx via fink.
An input file suitable for testing weak scattering is included "weak.in". The key parameters are mb2=0.5, dx=0.1, Rf=1, n=256, nreal=10, nscreen=1, nf=1, beam='n'. Simulate this using "simpw weak" and check that the variance is close to 0.44. Examine the output with "iplot weak" and observe that the window is large enough to display the intensity fluctuations well. Then look at the pdf with "pdfm weak". You should see a quasi gaussian with substantial skewness. Finally look at the spatial correlations and spectra of intensity with "corr2dn weak". This will show the isotropic correlation with a spatial scale of Rf and no significant large scale component. The spatial spectra will show a clearly defined low frequency region and a power law high frequency slope of -8/3.
Test the dynamic spectrum by increasing nf to 128 and reducing nreal to 1. Then create a dynamic spectrum with "getdyn weak" and plot it with "plotdyn weak". You will see a clear secondary spectrum. You might want to replot with contrast = 70 (dB). The default astronomical parameters are set for the solar wind with V = 400 km/s. You can overplot an arc with this velocity as a test. This is probably the minimum sampling you can get away with in computing a secondary spectrum. You can halve the sample interval to dx=0.05 and keep the window the same using n=512. Then you need to cut dlam/nf by a factor of 4. Say nf=256, and dlam=0.5. This will give a much sharper secondary spectrum.
Test the same dynamic spectrum with equal frequency steps, i.e. set lamscale='n'. The resulting secondary spectrum will not have as sharp an edge as before. The difference becomes less noticeable when dlam is smaller.
Test the thick screen simulation with a dynamic spectrum by setting ns=10 and toverz=2.0. This will provide a uniform medium all the way to the observer. Remember to set lamscale = 'y' again. Resimulate the dynamic spectrum. You will be able to see each of the 10 screens contributing to the secondary spectrum. This demonstrates the ability of the secondary spectrum to resolve thick vs thin screens. You will not be able to see much difference in the shape of the correlation function or the pdf.
Test the spherical wave simulator with the same dynamic spectrum. You will again see each of the screens in the secondary spectrum, and you will observe that they are better resolved than in the case of the plane wave. This dynamic spectrum is a simulation of what you would observe if the pulsar velocity is much larger than the Earth's velocity.
An input file suitable for testing stronger scattering is included "strong.in". Here the key parameters are changed to mb2=5.0, dx=0.02, nreal=10, nf=1, and n=1024. You will see with iplot that the intensity fluctuation scale is much smaller, but the fluctuations cluster in larger groups. The pdf looks almost exponential. The spatial correlations show two scales. The smaller scale is about 0.3 Rf and there is a larger scale carrying about 20% of the total variance with a scale of about 3 Rf. The larger scale is more variable. This is just estimation error. There are fewer degrees of freedom in sampling the larger (refractive) scale. The spatial spectra are no longer flat at low frequencies due to the larger scale component. We would need a larger window to resolve the large component fully.
Test the dynamic spectrum for this case. Since dx has been reduced by a factor of 2.5 to sample the diffractive scale adequately you need to reduce dlam/nf by roughly a factor of 6 to get a good sampling of the secondary spectrum. Since the calculation will run 16 times slower you should reduce dlam to 0.1 and leave nf=256. The output file will also be much larger at 8 MB/window*256 windows = 2 GB. Check your free disc space. The secondary spectrum in strong scattering lacks the razor edge seen in weak scattering, but it is still well defined at mb2 = 5.0. You will see the alias of the parabolic arc reflected in the right wall of the secondary spectrum box.
The output file has an ascii header which can be read using a command like "head -20 filename.out". The first line is the name of the simulator "simpw". The second line is the date and time of the simulation. The remainder should be the same as filename.in. The length of the header is equal to 8n bytes, the length of one row of the simulation. If n < 128 this will cause a problem because there are 19 lines of length 32 at present. We will feel free to increase this.
2. The spherical wave code "simsw" breaks the medium into ns screens which completely fill the space between the source and the observer. Each slab has the same level of turbulence and thickness. The simulation is done in scaled coordinate system in which the transverse variables are scaled by the distance from the source. The individual screen then are all scaled by different factors. The random screen generator produces a pair of screens. Both are used in simpw, but only one can be used in simsw. This makes simsw a bit slower. It uses an input file and creates an output file which are the same as for simpw except that the program name is different.
3. The raw intensity viewer iplot simply read the output file and extracts the simulated electric field windows one after another. It plots the squared magnitude of the electric field for each window. The mean intensity is unity, so one can get an idea of the intensity of the peaks using the gray scale wedge attached to the plot. Like all the analysis programs it does not need the filename.in parameter file, it gets the header information from the filename.out file.
4. The pdf code reads the first nreal windows and computes a histogram of intensity normalized to make a pdf. It also computes the moments of the pdf and error estimates. If both nreal and nf are greater than one it fails because the first nreal windows will not have the same frequency.
5. The corr2 code computes a spatial correlation of intensity for each of the first nreal windows and averages them together. It does this by computing the intensity in each window, then a 2-d fft, then taking the squared magnitude, then taking an inverse 2-d fft. These correlations are then averaged. The squared magnitudes, which are the spectra, could have been averaged and only one inverse fft formed. The way it is done was convenient during testing and its efficiency doesn't seem to be a problem. It also fails if both nreal and nf are greater than 1. The 2-d correlation is plotted first. Then a pair of cuts in the x and y directions are plotted with the corresponding spectra.
6. Getdyn reads filename.out and creates filename.spc which contains the header and binary dynamic spectrum. The header may occupy multiple lines because the dynamic spectrum is stored by frequency. There are only nf*4 bytes per line, where as the header is n*8 bytes long. The entire header is prepended to the binary data. The middle line of the window is taken. Getdyn does not distinguish between beamed and unbeamed simulations, that must be done in the secondary analysis code.
7. Plotdyn reads the header and data from filename.spc, then plots the dynamic spectrum and the secondary spectrum. The user is then asked to input various options. If the data are replotted the user can change the dynamic range (in dB) of the secondary spectrum. If the velocity is called for an arc with that velocity will be plotted. The user can also plot vertical cuts through the secondary spectrum. These are averaged over a range of delay taking into account the parabolic curvature in that range.
The code is not public domain software, however it is freely available for non-commercial use. We would prefer that you not distribute modified code. If you make a modification that could be generally useful, let us know and we will try to add it to the code. Please let us know about errors, omissions, and desired features.