mcmlio.c

经典的蒙特卡罗程序

/***********************************************************
 *  Copyright Univ. of Texas M.D. Anderson Cancer Center
 *  1992.
 *
 *	Input/output of data.
 ****/

#include "mcml.h"

/***********************************************************
 *	Structure used to check against duplicated file names.
 ****/
struct NameList {
  char name[STRLEN];
  struct NameList * next;
};

typedef struct NameList NameNode;
typedef NameNode * NameLink;


/***********************************************************
 *	Center a string according to the column width.
 ****/
char * CenterStr(short  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);
}

/**************************************************************************
 *	Center a string according to the column width.
 ****/
#define COLWIDTH 80
void
CtrPuts(char *InStr)
{
  short       nspaces;		/* number of spaces to be left-filled. */
  char        outstr[STRLEN];

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

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

  strcat(outstr, InStr);

  puts(outstr);
}

/***********************************************************
 *	Print some messages before starting simulation.  
 *	e.g. author, address, program version, year.
 ****/
#define COLWIDTH 80
void ShowVersion(char *version)
{
  char str[STRLEN];
  
  CtrPuts(" ");
  CtrPuts(
  "Monte Carlo Simulation of Light Transport in Multi-layered Turbid Media");
  CtrPuts(version);
  CtrPuts(" ");

  CtrPuts("Lihong Wang, Ph.D.");
  CtrPuts("Biomedical Engineering Program, Texas A&M University, 3120 TAMU");
  CtrPuts("College Station, Texas 77843-3120, USA");
  CtrPuts("Email: [e-mail]");

  CtrPuts(" ");
  CtrPuts("Steven L. Jacques, Ph.D.");
  CtrPuts("Oregon Medical Laser Center");
  CtrPuts("Providence/St. Vincent Hospital");
  CtrPuts("9205 SW Barnes Rd., Portland, OR 97225, USA");
  CtrPuts("Email: [e-mail]");

  CtrPuts(" ");
  CtrPuts("The program can be obtained from http://oilab.tamu.edu");
  CtrPuts("Please cite the following article in your publications:");
  printf("\tL.-H. Wang, S. L. Jacques, and L.-Q. Zheng, MCML - Monte \n");
  printf("\tCarlo modeling of photon transport in multi-layered\n");
  printf("\ttissues, Computer Methods and Programs in Biomedicine, 47,\n");
  printf("\t131-146 (1995)\n");

  puts("\n");	
}
#undef COLWIDTH

/***********************************************************
 *	Get a filename and open it for reading, retry until 
 *	the file can be opened.  '.' terminates the program.
 *      
 *	If Fname != NULL, try Fname first.
 ****/
FILE *GetFile(char *Fname)
{
  FILE * file=NULL;
  Boolean firsttime=1;
  
  do {
    if(firsttime && Fname[0]!='\0') { 
	  /* use the filename from command line */
      firsttime = 0;
    }
    else {
      printf("Input filename(or . to exit):");
      scanf("%s", Fname);
      firsttime = 0;
    }

    if(strlen(Fname) == 1 && Fname[0] == '.') 
      exit(1);			/* exit if no filename entered. */
    
    file = fopen(Fname, "r");
  }  while(file == NULL);
  
  return(file);
}

/***********************************************************
 *	Kill the ith char (counting from 0), push the following 
 *	chars forward by one.
 ****/
void KillChar(size_t i, char * Str)
{
  size_t sl = strlen(Str);
  
  for(;i<sl;i++) Str[i] = Str[i+1];
}

/***********************************************************
 *	Eliminate the chars in a string which are not printing 
 *	chars or spaces.
 *
 *	Spaces include ' ', '\f', '\t' etc.
 *
 *	Return 1 if no nonprinting chars found, otherwise 
 *	return 0.
 ****/
Boolean CheckChar(char * Str)
{
  Boolean found = 0;	/* found bad char. */
  size_t sl = strlen(Str);
  size_t i=0;
  
  while(i<sl) 
    if (Str[i]<0 || Str[i]>255)
      nrerror("Non-ASCII file\n");
    else if(isprint(Str[i]) || isspace(Str[i])) 
      i++;
    else {
      found = 1;
      KillChar(i, Str);
      sl--;
    }
  
  return(found);	
}

/***********************************************************
 *	Return 1 if this line is a comment line in which the 
 *	first non-space character is "#".
 *
 *	Also return 1 if this line is space line.
 ****/
Boolean CommentLine(char *Buf)
{
  size_t spn, cspn;
  
  spn = strspn(Buf, " \t");	
  /* length spanned by space or tab chars. */

  cspn = strcspn(Buf, "#\n");
  /* length before the 1st # or return. */

  if(spn == cspn) 	/* comment line or space line. */
	return(1);
  else				/* the line has data. */	 
	return(0);		
}

/***********************************************************
 *	Skip space or comment lines and return a data line only.
 ****/
char * FindDataLine(FILE *File_Ptr)
{
  static char buf[STRLEN];	/* LW 1/11/2000. Added static. */
  
  buf[0] = '\0';
  do {	/* skip space or comment lines. */
    if(fgets(buf, 255, File_Ptr) == NULL)  {
      printf("Incomplete data\n");
      buf[0]='\0';
      break;
    }
    else
      CheckChar(buf);
  } while(CommentLine(buf));
  
  return(buf);
}

/***********************************************************
 *	Skip file version, then read number of runs.
 ****/
short ReadNumRuns(FILE* File_Ptr)
{
  char buf[STRLEN];
  short n=0;
  
  FindDataLine(File_Ptr); /* skip file version. */

  strcpy(buf, FindDataLine(File_Ptr));
  if(buf[0]=='\0') nrerror("Reading number of runs\n");
  sscanf(buf, "%hd",&n);	        
  return(n);
}

  
/***********************************************************
 *	Read the file name and the file format.
 *
 *	The file format can be either A for ASCII or B for
 *	binary.
 ****/
void ReadFnameFormat(FILE *File_Ptr, InputStruct *In_Ptr)
{
  char buf[STRLEN];

  /** read in file name and format. **/
  strcpy(buf, FindDataLine(File_Ptr));
  if(buf[0]=='\0') 
	nrerror("Reading file name and format.\n");
  sscanf(buf, "%s %c", 
	In_Ptr->out_fname, &(In_Ptr->out_fformat) );
  if(toupper(In_Ptr->out_fformat) != 'B') 
	In_Ptr->out_fformat = 'A';
}


/***********************************************************
 *	Read the number of photons.
 ****/
void ReadNumPhotons(FILE *File_Ptr, InputStruct *In_Ptr)
{
  char buf[STRLEN];

  /** read in number of photons. **/
  strcpy(buf, FindDataLine(File_Ptr));
  if(buf[0]=='\0') 
	nrerror("Reading number of photons.\n");
  sscanf(buf, "%ld", &In_Ptr->num_photons);
  if(In_Ptr->num_photons<=0) 
	nrerror("Nonpositive number of photons.\n");
}


/***********************************************************
 *	Read the members dz and dr.
 ****/
void ReadDzDr(FILE *File_Ptr, InputStruct *In_Ptr)
{
  char buf[STRLEN];

  /** read in dz, dr. **/
  strcpy(buf, FindDataLine(File_Ptr));
  if(buf[0]=='\0') nrerror("Reading dz, dr.\n");
  sscanf(buf, "%lf%lf", &In_Ptr->dz, &In_Ptr->dr);
  if(In_Ptr->dz<=0) nrerror("Nonpositive dz.\n");
  if(In_Ptr->dr<=0) nrerror("Nonpositive dr.\n");
}


/***********************************************************
 *	Read the members nz, nr, na.
 ****/
void ReadNzNrNa(FILE *File_Ptr, InputStruct *In_Ptr)
{
  char buf[STRLEN];

  /** read in number of dz, dr, da. **/
  strcpy(buf, FindDataLine(File_Ptr));
  if(buf[0]=='\0') 
	nrerror("Reading number of dz, dr, da's.\n");
  sscanf(buf, "%hd%hd%hd", 
	&In_Ptr->nz, &In_Ptr->nr, &In_Ptr->na);
  if(In_Ptr->nz<=0) 
	nrerror("Nonpositive number of dz's.\n");
  if(In_Ptr->nr<=0) 
	nrerror("Nonpositive number of dr's.\n");
  if(In_Ptr->na<=0) 
	nrerror("Nonpositive number of da's.\n");
  In_Ptr->da = 0.5*PI/In_Ptr->na;
}


/***********************************************************
 *	Read the number of layers.
 ****/
void ReadNumLayers(FILE *File_Ptr, InputStruct *In_Ptr)
{
  char buf[STRLEN];

  /** read in number of layers. **/
  strcpy(buf, FindDataLine(File_Ptr));
  if(buf[0]=='\0') 
	nrerror("Reading number of layers.\n");
  sscanf(buf, "%hd", &In_Ptr->num_layers);
  if(In_Ptr->num_layers<=0) 
	nrerror("Nonpositive number of layers.\n");
}


/***********************************************************
 *	Read the refractive index n of the ambient.
 ****/
void ReadAmbient(FILE *File_Ptr, 
				 LayerStruct * Layer_Ptr,
				 char *side)
{
  char buf[STRLEN], msg[STRLEN];
  double n;

  strcpy(buf, FindDataLine(File_Ptr));
  if(buf[0]=='\0') {
    sprintf(msg, "Rading n of %s ambient.\n", side);
	nrerror(msg);
  }

  sscanf(buf, "%lf", &n );
  if(n<=0) nrerror("Wrong n.\n");
  Layer_Ptr->n = n;
}


/***********************************************************
 *	Read the parameters of one layer.
 *
 *	Return 1 if error detected.
 *	Return 0 otherwise.
 *
 *	*Z_Ptr is the z coordinate of the current layer, which
 *	is used to convert thickness of layer to z coordinates
 *	of the two boundaries of the layer.
 ****/
Boolean ReadOneLayer(FILE *File_Ptr, 
					 LayerStruct * Layer_Ptr,
					 double *Z_Ptr)
{
  char buf[STRLEN], msg[STRLEN];
  double d, n, mua, mus, g;	/* d is thickness. */

  strcpy(buf, FindDataLine(File_Ptr));
  if(buf[0]=='\0') return(1);	/* error. */

  sscanf(buf, "%lf%lf%lf%lf%lf", &n, &mua, &mus, &g, &d);
  if(d<0 || n<=0 || mua<0 || mus<0 || g<0 || g>1) 
    return(1);			/* error. */
    
  Layer_Ptr->n	= n;
  Layer_Ptr->mua = mua;	
  Layer_Ptr->mus = mus;	
  Layer_Ptr->g   = g;
  Layer_Ptr->z0	= *Z_Ptr;
  *Z_Ptr += d;
  Layer_Ptr->z1	= *Z_Ptr;

  return(0);
}

/***********************************************************
 *	Read the parameters of one layer at a time.
 ****/
void ReadLayerSpecs(FILE *File_Ptr, 
					short Num_Layers,
					LayerStruct ** Layerspecs_PP)
{
  char msg[STRLEN];
  short i=0;
  double z = 0.0;	/* z coordinate of the current layer. */
  
  /* Allocate an array for the layer parameters. */
  /* layer 0 and layer Num_Layers + 1 are for ambient. */
  *Layerspecs_PP = (LayerStruct *)
    malloc((unsigned) (Num_Layers+2)*sizeof(LayerStruct));
  if (!(*Layerspecs_PP)) 
    nrerror("allocation failure in ReadLayerSpecs()");
  
  ReadAmbient(File_Ptr, &((*Layerspecs_PP)[i]), "top"); 
  for(i=1; i<=Num_Layers; i++)  
    if(ReadOneLayer(File_Ptr, &((*Layerspecs_PP)[i]), &z)) {
      sprintf(msg, "Error reading %hd of %hd layers\n", 
			  i, Num_Layers);
	  nrerror(msg);
    }
  ReadAmbient(File_Ptr, &((*Layerspecs_PP)[i]), "bottom"); 
}

/***********************************************************
 *	Compute the critical angles for total internal
 *	reflection according to the relative refractive index
 *	of the layer.
 *	All layers are processed.
 ****/
void CriticalAngle( short Num_Layers, 
					LayerStruct ** Layerspecs_PP)
{
  short i=0;
  double n1, n2;
  
  for(i=1; i<=Num_Layers; i++)  {
    n1 = (*Layerspecs_PP)[i].n;
    n2 = (*Layerspecs_PP)[i-1].n;
    (*Layerspecs_PP)[i].cos_crit0 = n1>n2 ? 
		sqrt(1.0 - n2*n2/(n1*n1)) : 0.0;
    
    n2 = (*Layerspecs_PP)[i+1].n;
    (*Layerspecs_PP)[i].cos_crit1 = n1>n2 ? 
		sqrt(1.0 - n2*n2/(n1*n1)) : 0.0;
  }
}

/***********************************************************
 *	Read in the input parameters for one run.
 ****/
void ReadParm(FILE* File_Ptr, InputStruct * In_Ptr)
{
  char buf[STRLEN];
  
  In_Ptr->Wth = WEIGHT;
  
  ReadFnameFormat(File_Ptr, In_Ptr);
  ReadNumPhotons(File_Ptr, In_Ptr);
  ReadDzDr(File_Ptr, In_Ptr);
  ReadNzNrNa(File_Ptr, In_Ptr);
  ReadNumLayers(File_Ptr, In_Ptr);

  ReadLayerSpecs(File_Ptr, In_Ptr->num_layers, 
				&In_Ptr->layerspecs);
  CriticalAngle(In_Ptr->num_layers, &In_Ptr->layerspecs);
}

/***********************************************************
 *	Return 1, if the name in the name list.
 *	Return 0, otherwise.
 ****/
Boolean NameInList(char *Name, NameLink List)
{
  while (List != NULL) {
	if(strcmp(Name, List->name) == 0) 
	  return(1);
	List = List->next;
  };
  return(0);
}

/***********************************************************
 *	Add the name to the name list.
 ****/
void AddNameToList(char *Name, NameLink * List_Ptr)
{
  NameLink list = *List_Ptr;

  if(list == NULL) {	/* first node. */
	*List_Ptr = list = (NameLink)malloc(sizeof(NameNode));
	strcpy(list->name, Name);
	list->next = NULL;
  }
  else {				/* subsequent nodes. */
	/* Move to the last node. */
	while(list->next != NULL)
	  list = list->next;

	/* Append a node to the list. */
	list->next = (NameLink)malloc(sizeof(NameNode));
	list = list->next;
	strcpy(list->name, Name);
	list->next = NULL;
  }
}

/***********************************************************
 *	Check against duplicated file names.
 *
 *	A linked list is set up to store the file names used
 *	in this input data file.
 ****/
Boolean FnameTaken(char *fname, NameLink * List_Ptr)
{
  if(NameInList(fname, *List_Ptr))
	return(1);
  else {
	AddNameToList(fname, List_Ptr);
	return(0);
  }
}

/***********************************************************
 *	Free each node in the file name list.
 ****/
void FreeFnameList(NameLink List)
{
  NameLink next;

  while(List != NULL) {
	next = List->next;
	free(List);
	List = next;
  }
}

/***********************************************************
 *	Check the input parameters for each run.
 ****/
void CheckParm(FILE* File_Ptr, InputStruct * In_Ptr)
{
  short i_run;
  short num_runs;	/* number of independent runs. */
  NameLink head = NULL;
  Boolean name_taken;/* output files share the same */
					/* file name.*/
  char msg[STRLEN];
  
  num_runs = ReadNumRuns(File_Ptr);
  for(i_run=1; i_run<=num_runs; i_run++)  {
    printf("Checking input data for run %hd\n", i_run);
    ReadParm(File_Ptr, In_Ptr);

	name_taken = FnameTaken(In_Ptr->out_fname, &head);
	if(name_taken) 
	  sprintf(msg, "file name %s duplicated.\n", 
					In_Ptr->out_fname);

    free(In_Ptr->layerspecs);
	if(name_taken) nrerror(msg);
  }
  FreeFnameList(head);
  rewind(File_Ptr);
}


/***********************************************************
 *	Allocate the arrays in OutStruct for one run, and 
 *	array elements are automatically initialized to zeros.
 ****/
void InitOutputData(InputStruct In_Parm, 
					OutStruct * Out_Ptr)
{
  short nz = In_Parm.nz;
  short nr = In_Parm.nr;
  short na = In_Parm.na;
  short nl = In_Parm.num_layers;	
  /* remember to use nl+2 because of 2 for ambient. */
  
  if(nz<=0 || nr<=0 || na<=0 || nl<=0)