by David Boas

The following set of programs will calculate the distortion of diffuse photon density waves (DPDW) by spherical or cylindrical objects embedded in either an infinite homogeneous medium or a homogeneous medium with planar free-space boundaries. The theory is described in a paper by Boas, O'Leary, Chance, and Yodh, which is to be published in the Proceedings of the National Academy of Sciences in 1994. Please contact me if you have any questions or comments, or if you uncover any bugs.

Download Source Code

FILES

p1.c
Main program. Inputs all the necessary parameters and initializes the calculation.
p1.h
Header file containing all global declerations, etc.
dpdw.c
Calculates the DPDW for general source-arrays, detector arrays, planar boundaries, and object distributions.
helmholtz.c
Contains functions which are solutions of the 3D Helmholtz equation in spherical or cylindrical coordinates.
cbess.c
Spherical Bessel functions and Modified bessel Functions with complex arguments.
compfunc.c
Some useful functions with complex arguments.
dcomplex.c
More useful functions with complex arguments.
quad.c
Integration routine by quadrature (from Numerical Recipes).
nrutil.c
Some useful utilities from Numerical Recipes.
nrutil.h
Numerical Recipes header file.
makefile
This file is needed for compiling the code into an executable.

If one of these files is missing, then you will not be able to create the executable.

INSTALLATION

I describe here the method for installing this program on a unix machine. Installation on other machines should be no problem as long as the c-compiler is ANSI standard.

  1. Copy all the files into a common directory
  2. Compile the executable with the make command, i.e type 'make' at the command line.
  3. If all works well, this will create an executable called 'p1' which can do many calculations related to the scattering of DPDW's from spherical and cylindrical objects.

Feel free to contact me if this installation has not worked or if you are interested in installing this program on another system.

USAGE

To calculate the distortion/scattering of DPDW's in a model system, the system has to be characterized for the program. We characterize a system by telling the computer all of the simulation parameters including:

  • number of objects
  • number of sources
  • number of detectors
  • number of free space boundaries
  • reduced scattering and absorption coef. of the homogeneous medium
  • index of refraction of the homogeneous medium
  • modulation frequency
  • reduced scattering and absorption coef. of each object
  • position and size of each object
  • position of sources
  • position of detectors

There are two different methods of entering the simulation parameters into the program. One is to answer a series of questions asked by the program at run time. The other is to prepare two input files of a format that is described below. The run time input method is best for calculating on a 2 or 3 dimensional grid, the distortion/scattering of a DPDW generated by a fixed source, detected by one detector and scattered by any number of objects and/or boundaries. On the other hand, the input files are best when doing several calculations for which the source is not necessarily fixed.

Below I will run through examples of using either input method.

RUN TIME INPUT

Let's calculate the DPDW generated by a phased-array of two sources and scattered by one absorbing sphere. Output from the computer is enclosed in quotes while your input is indented.

To start the program, type

    p1 s

If we were doing the calculation for a cylinder, we would instead type p1 c

We are first asked

    "Number of sources, objects, and free-space planar boundaries:"

Respond with

    2 1 0

Next, we are asked about the homogeneous medium

"Homogeneous Medium
      mu_s' and mu_a (units of inverse centimeters):"

Respond with

10 .02

All length units are in centimeters.

Next,

"Index of Refraction of the turbid medium:"
    1.333

Next

"Modulation Frequency in MHz"  

This is the frequency at which the intensity of your source is modulated, respond in MHz, e.g.

    220

We are then asked about the object.

    "Object 1: mu_s' mu_a  
               Position Xo Yo Zo
               Size a (units of centimeters):" 
respond
    10 5 0 4 0 .5

This object has a reduced scattering coef. of 10 inverse centimeters, an absorption coef. of 5 inverse centimeters. It is positioned at Xo=0.0 cm, Yo=4.0 cm, Zo=0.0 cm, and it has a radius of 0.5 cm.

Next, the sources.

    "Source 1: Position  Xs Ys Zs 
               Initial Amplitude and Phase A_o Phi_o:"
    -.5 0 0 1 0
    "Source 2: Position  Xs Ys Zs 
               Initial Amplitude and Phase A_o Phi_o:"
    .5 0 0 1 180

The first source is at Xs=-0.5 cm, Ys=0.0 cm, Zs=0.0 cm, and has an initial amplitude of 1 and phase of 0 degrees. The second source is at Xs=0.5 cm, Ys=0.0 cm, Zs=0.0 cm, and has an initial amplitude of 1 and phase of 180 degrees.

Next, we are asked about the name for the output file.

    "Output file name (excluding extension):"  
Respond with your choosen name
sample

The program will automatically place the proper extension on the file name. It does this because it produces 3 output files, one containing the ampitude and phase of the distorted wave, one containing the amplitude and phase of the scattered wave, and the final containing the X Y Z position of the detector. The format of the output files will be described below.

Finally, you are asked about the detector positions.

    "Detector Positions (define a 3d grid):"
    " Xmin, Xmax, Xstep:"
    -3 0 1
    " Ymin, Ymax, Ystep:"
    6 6 1
    " Zmin, Zmax, Zstep:"
    0 0 1
    

Thus, the DPDW will be calculated at (X=-3.0 cm, Y=6.0, Z=0.0), (X=-2.0, Y=6.0, Z=0.0), (X=-1.0, Y=6.0, Z=0.0), (X=0.0, Y=6.0, Z=0.0).

FILE INPUT

Let's do the same calculation as above, but this time we will input the simulation parameters using input files. Output from the computer is enclosed in quotes while your input is indented. To start the program, type

    p1 s source_detector_file phantom_file output_file

The 's' indicates that we are doing the calculation for a sphere. The next three items on the command line are the names of the two input files and the name for the output file (without extensions). For this calculation, the source_detector_file is 'source_detector_file.sdg'

            2   1
            220
            7
            -.5 0 0 1 0     .5 0 0 1 180    -3.0 6 0 1 0
            -.5 0 0 1 0     .5 0 0 1 180    -2.5 6 0 1 0
            -.5 0 0 1 0     .5 0 0 1 180    -2.0 6 0 1 0
            -.5 0 0 1 0     .5 0 0 1 180    -1.5 6 0 1 0
            -.5 0 0 1 0     .5 0 0 1 180    -1.0 6 0 1 0
            -.5 0 0 1 0     .5 0 0 1 180    -0.5 6 0 1 0
            -.5 0 0 1 0     .5 0 0 1 180     0.0 6 0 1 0

The first line contains the number of sources and the number of detectors used for each measurement. The next line is the modulation freq. of the source in MHz. Then the number of measurements followed by the Xs Ys Zs A_o Phi_o for each source and Xd Yd Zd A_o Phi_o for each detector. The phantom_file is 'phantom_file.phm'

            1    0
            10   .02    1.33333
            10   5   0   4   0   .5

The first line indicates the number of objects and number of planar boundaries. The second line has the homogeneous mu_s', mu_a and the index of refraction of the turbid medium. The next line is for the first object, mu_s' mu_a Xo Yo Zo a. If there are more objects, additional lines would be added for each object.

If the format of your input files is correct, then the program will execute properly and return two output files; the output.pdw file containing the amplitude and phase of the distorted wave, abd the output.sc file containing the amplitude and phase of the scattered wave.

OUTPUT FILES

The output files from the above sample seesions are given below.

The 'sample.pdw' file will be

    7
    5.455546e-15 -78.361698
    6.085271e-15 -87.653065
    6.212503e-15 -95.357244
    5.639028e-15 -101.307615
    4.299213e-15 -105.382961
    2.323047e-15 -107.630010
    5.057899e-24 -162.986157

where the first line is the number of measurements. Each additional line is the amplitude and phase for each measurement.

The 'sample.sc' file will be

    7
    4.761850e-17 113.683236
    7.980468e-17 95.978313
    1.251892e-16 79.943429
    1.746208e-16 66.197048
    1.961558e-16 55.511084
    1.403547e-16 48.691068
    3.136315e-24 5.109805

The format is the same as for 'sample.pdw'.

If the simulation parameters where entered at run time, then the additional file 'sample.xyz' will be created

   
    7
    1
    1
    -3.000000 6.000000 0.000000
    -2.500000 6.000000 0.000000
    -2.000000 6.000000 0.000000
    -1.500000 6.000000 0.000000
    -1.000000 6.000000 0.000000
    -0.500000 6.000000 0.000000
     0.000000 6.000000 0.000000

The first three lines indicate the number of grid points in the x, y, and z directions respectively. Following will be (nx*ny*nz) lines giving the position of the detector, Xd Yd Zd, for each measurement. nx, ny, and nz are the number of grid points along each axis.