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];

  CenterStr(COLWIDTH,
	    "Convolution of mcml Simulation Data", str);
  puts("");
  puts(str);
  puts("");

  CenterStr(COLWIDTH, "Lihong Wang, Ph.D.", str);
  puts(str);

  CenterStr(COLWIDTH, "Steven L. Jacques, Ph.D.", str);
  puts(str);

  CenterStr(COLWIDTH, "Laser Biology Research Laboratory - 017", str);
  puts(str);

  CenterStr(COLWIDTH, "University of Texas M.D. Anderson Cancer Center", str);
  puts(str);

  CenterStr(COLWIDTH, "Houston, Texas 77030, USA", str);
  puts(str);

  CenterStr(COLWIDTH, "Voice: [phone], Fax: [phone]", str);
  puts(str);

  puts("");

  CenterStr(COLWIDTH, "Version 1.1, 1994", str);
  puts(str);
  puts("\n\n");
}
#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]);

  fclose(file);
}

/****************************************************************
 *	2 numbers each line: a, Rd[a].
 ****/
void 
WriteRd_a(InputStruct * In_Ptr,
	  double *Rd_a)
{
  short       ia, na = In_Ptr->na;
  double      da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

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

  fprintf(file, "%-12s\t%-s[1/sr]\n", "a[rad]", fname);
  for (ia = 0; ia < na; ia++)
    fprintf(file, "%-12.4E\t%-12.4E\n", (ia + 0.5) * da,
	    Rd_a[ia]);

  fclose(file);
}

/****************************************************************
 *	3 numbers each line:r, a, Tt[r][a]. a = theta.
 ****/
void 
WriteTt_ra(InputStruct * In_Ptr,
	   double **Tt_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, "Tra");
  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, Tt_ra[ir][ia]);
    }
  }

  fclose(file);
}

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

  strcpy(fname, "Tr");
  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,
	    Tt_r[ir]);

  fclose(file);
}

/****************************************************************
 *	2 numbers each line: theta, Tt[theta].
 ****/
void 
WriteTt_a(InputStruct * In_Ptr,
	  double *Tt_a)
{
  short       ia, na = In_Ptr->na;
  double      da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

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

  fprintf(file, "%-12s\t%-s[1/sr]\n", "a[rad]", fname);
  for (ia = 0; ia < na; ia++)
    fprintf(file, "%-12.4E\t%-12.4E\n", (ia + 0.5) * da,
	    Tt_a[ia]);

  fclose(file);
}

/****************************************************************
 *	Write output in M. Keijzer's format so that the file can be
 *	read by the convolution program written by Keijzer in Pascal.
 ****/
void 
WriteKFormat(InputStruct * In_Ptr,
	     OutStruct * Out_Ptr)
{
  short       i, j;
  double      dz, dr;
  FILE       *file;
  char        fname[STRLEN];

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

  fputs("output.filename\n", file);
  fprintf(file, "%12hd layers\n", In_Ptr->num_layers);
  fprintf(file, "%12s %12s %12s %12s %12s %12s\n",
	  "layer", "mua", "mus", "g", "nt", "thickness");

  for (i = 1; i <= In_Ptr->num_layers; i++)
    fprintf(file,
	    "%12hd %12.6lf %12.6lf %12.6lf %12.6lf %12.6lf\n",
	    i, In_Ptr->layerspecs[i].mua,
	    In_Ptr->layerspecs[i].mus,
	    In_Ptr->layerspecs[i].g, In_Ptr->layerspecs[i].n,
	    In_Ptr->layerspecs[i].z1 -
	    In_Ptr->layerspecs[i].z0);
  fprintf(file, "%12.6lf index of refraction above\n",
	  In_Ptr->layerspecs[0].n);
  fprintf(file, "%12.6lf index of refraction below\n",
	  In_Ptr->layerspecs[i].n);
  fprintf(file, "\n");

  fprintf(file, "%12ld photons\n", In_Ptr->num_photons);
  fprintf(file, "%12.6lf critical weight\n", In_Ptr->Wth);
  fprintf(file, "%12.6lf depth of boxes micron\n",
	  In_Ptr->dz * 1e4);
  fprintf(file, "%12.6lf width of boxes micron\n",
	  In_Ptr->dr * 1e4);
  fprintf(file, "%12hd number of boxes in z \n", In_Ptr->nz);
  fprintf(file, "%12hd number of boxes in r \n", In_Ptr->nr);
  fprintf(file, "\n");

  fprintf(file,
	  "%12.6lf Total reflection (including direct R)\n",
	  Out_Ptr->Rsp + Out_Ptr->Rd);
  for (i = 1; i <= In_Ptr->num_layers; i++)
    fprintf(file, "%12.6lf Absorbed in layer %12hd\n",
	    Out_Ptr->A_l[i], i);
  fprintf(file, "%12.6lf Total transmission\n", Out_Ptr->Tt);
  fprintf(file, "\n");

  fprintf(file, "Reflectance and Transmission in [cm-2]\n");
  fprintf(file, "Absorption in z-layers in [cm-1]\n");
  fprintf(file, "Absorption in z/r-boxes in [cm-3]\n");

  fprintf(file, "z/r [cm] Layer");
  dr = In_Ptr->dr;
  for (i = 0; i < In_Ptr->nr; i++)
    fprintf(file, "%12.6lf ", i * dr);
  fprintf(file, "\n");

  fprintf(file, "Refl. ");
  for (i = 0; i < In_Ptr->nr; i++)
    fprintf(file, "%12.6lf ", Out_Ptr->Rd_r[i]);
  fprintf(file, "\n");

  dz = In_Ptr->dz;
  for (i = 0; i < In_Ptr->nz; i++) {
    fprintf(file, "%12.6lf %12.6lf", i * dz, Out_Ptr->A_z[i]);
    for (j = 0; j < In_Ptr->nr; j++)
      fprintf(file, "%12.6lf ", Out_Ptr->A_rz[j][i]);
    fprintf(file, "\n");
  }

  fprintf(file, "Transm.");
  for (i = 0; i < In_Ptr->nr; i++)
    fprintf(file, "%12.6lf ", Out_Ptr->Tt_r[i]);
  fprintf(file, "\n");

  fclose(file);
}

/****************************************************************
 ****/
void 
ShowOutMenu(char *in_fname)
{
  printf("I   = Input parameters of mcml\n");
  printf("3   = reflectance, absorption, and transmittance\n");

  printf("AL  = absorption vs layer [-]\n");
  printf("Az  = absorption vs z [1/cm]\n");
  printf("Arz = absorption vs r & z [1/cm3]\n");