convo.c

卷积程序

/****************************************************************
 *	Functions for file output.
 ****/

#include "conv.h"

/****************************************************************
 *	Center a string according to the column width.
 ****/
char       *
CenterStr(short int Wid,
	  char *InStr,
	  char *OutStr)
{
  size_t      nspaces;
  /* number of spaces to be filled before InStr. */

  nspaces = (Wid - strlen(InStr)) / 2;
  if (nspaces < 0)
    nspaces = 0;

  strcpy(OutStr, "");
  while (nspaces--)
    strcat(OutStr, " ");

  strcat(OutStr, InStr);

  return (OutStr);
}

/****************************************************************
 *	Print some messages before starting simulation.  e.g. author,
 *	address, program version, year.
 ****/
#define COLWIDTH 80
void 
ShowVersion(void)
{
  char        str[STRLEN];

  puts("");
  CenterStr(COLWIDTH, "Convolution of MCML Simulation Data", str); puts(str);
  CenterStr(COLWIDTH, "Version 1.1, 1994", str); puts(str);
  puts("");

  CenterStr(COLWIDTH, "Lihong Wang, Ph.D.", str); puts(str);
  CenterStr(COLWIDTH, "Biomedical Engineering Program", str); puts(str);
  CenterStr(COLWIDTH, "3120 TAMU, Texas A&M University", str); puts(str);
  CenterStr(COLWIDTH, "College Station, TX 77843-3120", str); puts(str);
  CenterStr(COLWIDTH, "Email: [e-mail]", str); puts(str);
  CenterStr(COLWIDTH, "URL: http://oilab.tamu.edu", str); puts(str);
  puts("");

  CenterStr(COLWIDTH, "Steven L. Jacques, Ph.D.", str); puts(str);
  CenterStr(COLWIDTH, "Oregon Medical Laser Center", str); puts(str);
  CenterStr(COLWIDTH, "Providence/St. Vincent Hospital", str); puts(str);
  CenterStr(COLWIDTH, "9205 SW Barnes Rd., Portland, OR 97225", str); puts(str);
  CenterStr(COLWIDTH, "Email: [e-mail]", str); puts(str);
  CenterStr(COLWIDTH, "URL: http://omlc.ogi.edu/staff/jacques.html", str); 
  puts(str);
  puts("");

  CenterStr(COLWIDTH, 
  "Please cite the following article in your publications:", str); puts(str);
  CenterStr(COLWIDTH, 
  "L.-H. Wang, S. L. Jacques, and L.-Q. Zheng, ", str); puts(str);
  CenterStr(COLWIDTH, 
  "CONV - Convolution for responses to a finite diameter photon beam", str); 
  puts(str);
  CenterStr(COLWIDTH, 
  "incident on multi-layered tissues, Computer Methods and Programs in", str); 
  puts(str);
  CenterStr(COLWIDTH, 
  "Biomedicine 54, 141-150 (1997).", str); puts(str);
}
#undef COLWIDTH

/****************************************************************
 *	Open a file for output with extension Ext.  If file exists,
 *	ask whether to overwrite or append or change filename.
 *
 *	Return file pointer, which could be NULL.
 *	Return the full filename as Ext.
 ****/
FILE       *
GetWriteFile(char *Ext)
{
  FILE       *file;
  char        fname[STRLEN], fmode[STRLEN];

  do {
    printf("Enter output filename with extension .%s (or . to quit): ", Ext);
    gets(fname);
    if (strlen(fname) == 1 && fname[0] == '.') {
      fmode[0] = 'q';
      break;
    } else
      fmode[0] = 'w';

    if ((file = fopen(fname, "r")) != NULL) {	/* file exists. */
      printf("File %s exists, %s",
	     fname, "w=overwrite, a=append, n=new filename, q=quit: ");
      do
	gets(fmode);
      while (!strlen(fmode));	/* avoid null line. */
      fclose(file);
    }
  } while (fmode[0] != 'w' && fmode[0] != 'a' && fmode[0] != 'q');

  if (fmode[0] != 'q')
    file = fopen(fname, fmode);
  else
    file = NULL;

  strcpy(Ext, fname);
  return (file);		/* could be NULL. */
}

/****************************************************************
 *	Return the index to the layer, where iz is in.
 *
 *	Use the center of box.
 ****/
short 
IzToLayer(short Iz,
	  InputStruct * In_Ptr)
{
  short       i = 1;		/* index to layer. */
  short       num_layers = In_Ptr->num_layers;
  double      dz = In_Ptr->dz;

  while ((Iz + 0.5) * dz >= In_Ptr->layerspecs[i].z1
	 && i < num_layers)
    i++;

  return (i);
}

/****************************************************************
 *	Write the input parameter for Monte Carlo simulation program
 *	in such a format that it can be read directly back.
 ****/
void 
WriteInParm(InputStruct * In_Ptr)
{
  short       i;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "InP");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "1.1\t# file version.\n");
  fprintf(file, "1\t# number of runs.\n\n");
  fprintf(file, "temp.out\tA\t\t#output filename.\n");

  fprintf(file, "%ld \t\t\t# No. of photons.\n",
	  In_Ptr->num_photons);

  fprintf(file, "%G\t%G\t\t# dz, dr.\n", In_Ptr->dz, In_Ptr->dr);
  fprintf(file, "%hd\t%hd\t%hd\t# No. of dz, dr, da.\n\n",
	  In_Ptr->nz, In_Ptr->nr, In_Ptr->na);

  fprintf(file, "%hd\t\t\t\t\t# Number of layers.\n",
	  In_Ptr->num_layers);
  fprintf(file,
	  "#n\tmua\tmus\tg\td\t# One line for each layer.\n");
  fprintf(file, "%G\t\t\t\t\t# n for medium above.\n",
	  In_Ptr->layerspecs[0].n);
  for (i = 1; i <= In_Ptr->num_layers; i++) {
    LayerStruct s;
    s = In_Ptr->layerspecs[i];
    fprintf(file, "%G\t%G\t%G\t%G\t%G\t# layer %hd\n",
	    s.n, s.mua, s.mus, s.g, s.z1 - s.z0, i);
  }
  fprintf(file, "%G\t\t\t\t\t# n for medium below.\n\n",
	  In_Ptr->layerspecs[i].n);

  fclose(file);
}

/****************************************************************
 *	Write reflectance, absorption, transmission.
 ****/
void 
WriteRAT(OutStruct * Out_Ptr)
{
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "RAT");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file,			/* flag. */
	  "RAT #Reflectance, absorption, transmission.\n");
  fprintf(file, "%-12.4G \t#Specular reflectance.\n",
	  Out_Ptr->Rsp);
  fprintf(file, "%-12.4G \t#Diffuse reflectance.\n",
	  Out_Ptr->Rd);
  fprintf(file, "%-12.4G \t#Absorption.\n",
	  Out_Ptr->A);
  fprintf(file, "%-12.4G \t#Transmission.\n",
	  Out_Ptr->Tt);

  fprintf(file, "\n");

  fclose(file);
}

/****************************************************************
 *	Write absorption as a function of layer.
 *	2 numbers each line: layer, A[layer].
 ****/
void 
WriteA_layer(short Num_Layers,
	     double *A_l)
{
  short       i;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "Al");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "layer\t%-s[-]\n", fname);
  for (i = 1; i <= Num_Layers; i++)
    fprintf(file, "%-4hd\t%-12.4G\n", i, A_l[i]);

  fclose(file);
}

/****************************************************************
 *	2 numbers each line: z, A[z].
 ****/
void 
WriteA_z(InputStruct * In_Ptr, double *A_z)
{
  short       nz = In_Ptr->nz;
  double      dz = In_Ptr->dz;
  short       iz;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "Az");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-s[1/cm]\n", "z[cm]", fname);
  for (iz = 0; iz < nz; iz++)
    fprintf(file, "%-12.4E\t%-12.4E\n", (iz + 0.5) * dz, A_z[iz]);

  fclose(file);
}

/****************************************************************
 *	3 numbers each line: r, z, A[r][z].
 ****/
void 
WriteA_rz(InputStruct * In_Ptr,
	  double **A_rz)
{
  short       ir, iz, nz = In_Ptr->nz, nr = In_Ptr->nr;
  double      r, z, dr = In_Ptr->dr, dz = In_Ptr->dz;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "Arz");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-12s\t%-s[1/cm3]\n", "r[cm]", "z[cm]", fname);
  for (ir = 0; ir < nr; ir++) {
    r = (ir + 0.5) * dr;
    for (iz = 0; iz < nz; iz++) {
      z = (iz + 0.5) * dz;
      fprintf(file, "%-12.4E\t%-12.4E\t%-12.4E\n",
	      r, z, A_rz[ir][iz]);
    }
  }

  fclose(file);
}

/****************************************************************
 *	2 numbers each line: z, F[z].
 ****/
void 
WriteF_z(InputStruct * In_Ptr,
	 double *A_z)
{
  FILE       *file;
  short       iz, nz = In_Ptr->nz;
  double      mua, dz = In_Ptr->dz;
  char        fname[STRLEN];

  strcpy(fname, "Fz");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-s[-]\n", "z[cm]", fname);
  for (iz = 0; iz < nz; iz++) {
    mua = In_Ptr->layerspecs[IzToLayer(iz, In_Ptr)].mua;
    if (mua > 0.0)
      fprintf(file, "%-12.4E\t%-12.4E\n", (iz + 0.5) * dz,
	      A_z[iz] / mua);
    else
      fprintf(file, "%-12.4E\t%-12.4E\n", (iz + 0.5) * dz, 0.0);
  }

  fclose(file);
}

/****************************************************************
 *	3 numbers each line: r, z, F[r][z].
 ****/
void 
WriteF_rz(InputStruct * In_Ptr,
	  double **A_rz)
{
  FILE       *file;
  short       ir, iz, nz = In_Ptr->nz, nr = In_Ptr->nr;
  double      mua, r, z, dr = In_Ptr->dr, dz = In_Ptr->dz;
  char        fname[STRLEN];

  strcpy(fname, "Frz");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-12s\t%-s[1/cm2]\n",
	  "r[cm]", "z[cm]", fname);
  for (ir = 0; ir < nr; ir++) {
    r = (ir + 0.5) * dr;
    for (iz = 0; iz < nz; iz++) {
      z = (iz + 0.5) * dz;
      mua = In_Ptr->layerspecs[IzToLayer(iz, In_Ptr)].mua;
      if (mua > 0.0)
	fprintf(file, "%-12.4E\t%-12.4E\t%-12.4E\n",
		r, z, A_rz[ir][iz] / mua);
      else
	fprintf(file, "%-12.4E\t%-12.4E\t%-12.4E\n", r, z, 0.0);
    }
  }

  fclose(file);
}

/****************************************************************
 *	3 numbers each line: r, a, Rd[r][a].
 ****/
void 
WriteRd_ra(InputStruct * In_Ptr,
	   double **Rd_ra)
{
  short       ir, ia, nr = In_Ptr->nr, na = In_Ptr->na;
  double      r, a, dr = In_Ptr->dr, da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "Rra");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-12s\t%-s[1/(cm2sr)]\n",
	  "r[cm]", "a[rad]", fname);
  for (ir = 0; ir < nr; ir++) {
    r = (ir + 0.5) * dr;
    for (ia = 0; ia < na; ia++) {
      a = (ia + 0.5) * da;
      fprintf(file, "%-12.4E\t%-12.4E\t%-12.4E\n",
	      r, a, Rd_ra[ir][ia]);
    }
  }

  fclose(file);
}

/****
 *	2 numbers each line: r, Rd[r]
 ****/
void 
WriteRd_r(InputStruct * In_Ptr,
	  double *Rd_r)
{
  short       ir, nr = In_Ptr->nr;
  double      dr = In_Ptr->dr;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "Rr");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-s[1/cm2]\n", "r[cm]", fname);
  for (ir = 0; ir < nr; ir++)
    fprintf(file, "%-12.4E\t%-12.4E\n", (ir + 0.5) * dr,
	    Rd_r[ir]);