mcmlio.c

经典的蒙特卡罗程序

    nrerror("Wrong grid parameters.\n");
  
  /* Init pure numbers. */
  Out_Ptr->Rsp = 0.0;
  Out_Ptr->Rd  = 0.0;
  Out_Ptr->A   = 0.0;
  Out_Ptr->Tt  = 0.0;
  
  /* Allocate the arrays and the matrices. */
  Out_Ptr->Rd_ra = AllocMatrix(0,nr-1,0,na-1);
  Out_Ptr->Rd_r  = AllocVector(0,nr-1);
  Out_Ptr->Rd_a  = AllocVector(0,na-1);
  
  Out_Ptr->A_rz  = AllocMatrix(0,nr-1,0,nz-1);
  Out_Ptr->A_z   = AllocVector(0,nz-1);
  Out_Ptr->A_l   = AllocVector(0,nl+1);
  
  Out_Ptr->Tt_ra = AllocMatrix(0,nr-1,0,na-1);
  Out_Ptr->Tt_r  = AllocVector(0,nr-1);
  Out_Ptr->Tt_a  = AllocVector(0,na-1);
}

/***********************************************************
 *	Undo what InitOutputData did.
 *  i.e. free the data allocations.
 ****/
void FreeData(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. */
  
  free(In_Parm.layerspecs);
  
  FreeMatrix(Out_Ptr->Rd_ra, 0,nr-1,0,na-1);
  FreeVector(Out_Ptr->Rd_r, 0,nr-1);
  FreeVector(Out_Ptr->Rd_a, 0,na-1);
  
  FreeMatrix(Out_Ptr->A_rz, 0, nr-1, 0,nz-1);
  FreeVector(Out_Ptr->A_z, 0, nz-1);
  FreeVector(Out_Ptr->A_l, 0,nl+1);
  
  FreeMatrix(Out_Ptr->Tt_ra, 0,nr-1,0,na-1);
  FreeVector(Out_Ptr->Tt_r, 0,nr-1);
  FreeVector(Out_Ptr->Tt_a, 0,na-1);
}

/***********************************************************
 *	Get 1D array elements by summing the 2D array elements.
 ****/
void Sum2DRd(InputStruct In_Parm, OutStruct * Out_Ptr)
{
  short nr = In_Parm.nr;
  short na = In_Parm.na;
  short ir,ia;
  double sum;
  
  for(ir=0; ir<nr; ir++)  {
    sum = 0.0;
    for(ia=0; ia<na; ia++) sum += Out_Ptr->Rd_ra[ir][ia];
    Out_Ptr->Rd_r[ir] = sum;
  }
  
  for(ia=0; ia<na; ia++)  {
    sum = 0.0;
    for(ir=0; ir<nr; ir++) sum += Out_Ptr->Rd_ra[ir][ia];
    Out_Ptr->Rd_a[ia] = sum;
  }
  
  sum = 0.0;
  for(ir=0; ir<nr; ir++) sum += Out_Ptr->Rd_r[ir];
  Out_Ptr->Rd = sum;
}

/***********************************************************
 *	Return the index to the layer according to the index
 *	to the grid line system in z direction (Iz).
 *
 *	Use the center of box.
 ****/
short IzToLayer(short Iz, InputStruct In_Parm)
{
  short i=1;	/* index to layer. */
  short num_layers = In_Parm.num_layers;
  double dz = In_Parm.dz;
  
  while( (Iz+0.5)*dz >= In_Parm.layerspecs[i].z1 
	&& i<num_layers) i++;
  
  return(i);
}

/***********************************************************
 *	Get 1D array elements by summing the 2D array elements.
 ****/
void Sum2DA(InputStruct In_Parm, OutStruct * Out_Ptr)
{
  short nz = In_Parm.nz;
  short nr = In_Parm.nr;
  short iz,ir;
  double sum;
  
  for(iz=0; iz<nz; iz++)  {
    sum = 0.0;
    for(ir=0; ir<nr; ir++) sum += Out_Ptr->A_rz[ir][iz];
    Out_Ptr->A_z[iz] = sum;
  }
  
  sum = 0.0;
  for(iz=0; iz<nz; iz++) {
    sum += Out_Ptr->A_z[iz];
    Out_Ptr->A_l[IzToLayer(iz, In_Parm)] 
	  += Out_Ptr->A_z[iz];
  }
  Out_Ptr->A = sum;
}

/***********************************************************
 *	Get 1D array elements by summing the 2D array elements.
 ****/
void Sum2DTt(InputStruct In_Parm, OutStruct * Out_Ptr)
{
  short nr = In_Parm.nr;
  short na = In_Parm.na;
  short ir,ia;
  double sum;
  
  for(ir=0; ir<nr; ir++)  {
    sum = 0.0;
    for(ia=0; ia<na; ia++) sum += Out_Ptr->Tt_ra[ir][ia];
    Out_Ptr->Tt_r[ir] = sum;
  }
  
  for(ia=0; ia<na; ia++)  {
    sum = 0.0;
    for(ir=0; ir<nr; ir++) sum += Out_Ptr->Tt_ra[ir][ia];
    Out_Ptr->Tt_a[ia] = sum;
  }
  
  sum = 0.0;
  for(ir=0; ir<nr; ir++) sum += Out_Ptr->Tt_r[ir];
  Out_Ptr->Tt = sum;
}

/***********************************************************
 *	Scale Rd and Tt properly.
 *
 *	"a" stands for angle alpha.
 ****
 *	Scale Rd(r,a) and Tt(r,a) by
 *      (area perpendicular to photon direction)
 *		x(solid angle)x(No. of photons).
 *	or
 *		[2*PI*r*dr*cos(a)]x[2*PI*sin(a)*da]x[No. of photons]
 *	or
 *		[2*PI*PI*dr*da*r*sin(2a)]x[No. of photons]
 ****
 *	Scale Rd(r) and Tt(r) by
 *		(area on the surface)x(No. of photons).
 ****
 *	Scale Rd(a) and Tt(a) by
 *		(solid angle)x(No. of photons).
 ****/
void ScaleRdTt(InputStruct In_Parm, OutStruct *	Out_Ptr)
{
  short nr = In_Parm.nr;
  short na = In_Parm.na;
  double dr = In_Parm.dr;
  double da = In_Parm.da;
  short ir,ia;
  double scale1, scale2;
  
  scale1 = 4.0*PI*PI*dr*sin(da/2)*dr*In_Parm.num_photons;
	/* The factor (ir+0.5)*sin(2a) to be added. */

  for(ir=0; ir<nr; ir++)  
    for(ia=0; ia<na; ia++) {
      scale2 = 1.0/((ir+0.5)*sin(2.0*(ia+0.5)*da)*scale1);
      Out_Ptr->Rd_ra[ir][ia] *= scale2;
      Out_Ptr->Tt_ra[ir][ia] *= scale2;
    }
  
  scale1 = 2.0*PI*dr*dr*In_Parm.num_photons;  
	/* area is 2*PI*[(ir+0.5)*dr]*dr.*/ 
	/* ir+0.5 to be added. */

  for(ir=0; ir<nr; ir++) {
    scale2 = 1.0/((ir+0.5)*scale1);
    Out_Ptr->Rd_r[ir] *= scale2;
    Out_Ptr->Tt_r[ir] *= scale2;
  }
  
  scale1  = 2.0*PI*da*In_Parm.num_photons;
    /* solid angle is 2*PI*sin(a)*da. sin(a) to be added. */

  for(ia=0; ia<na; ia++) {
    scale2 = 1.0/(sin((ia+0.5)*da)*scale1);
    Out_Ptr->Rd_a[ia] *= scale2;
    Out_Ptr->Tt_a[ia] *= scale2;
  }
  
  scale2 = 1.0/(double)In_Parm.num_photons;
  Out_Ptr->Rd *= scale2;
  Out_Ptr->Tt *= scale2;
}

/***********************************************************
 *	Scale absorption arrays properly.
 ****/
void ScaleA(InputStruct In_Parm, OutStruct * Out_Ptr)
{
  short nz = In_Parm.nz;
  short nr = In_Parm.nr;
  double dz = In_Parm.dz;
  double dr = In_Parm.dr;
  short nl = In_Parm.num_layers;
  short iz,ir;
  short il;
  double scale1;
  
  /* Scale A_rz. */
  scale1 = 2.0*PI*dr*dr*dz*In_Parm.num_photons;	
	/* volume is 2*pi*(ir+0.5)*dr*dr*dz.*/ 
	/* ir+0.5 to be added. */
  for(iz=0; iz<nz; iz++) 
    for(ir=0; ir<nr; ir++) 
      Out_Ptr->A_rz[ir][iz] /= (ir+0.5)*scale1;
  
  /* Scale A_z. */
  scale1 = 1.0/(dz*In_Parm.num_photons);
  for(iz=0; iz<nz; iz++) 
    Out_Ptr->A_z[iz] *= scale1;
  
  /* Scale A_l. Avoid int/int. */
  scale1 = 1.0/(double)In_Parm.num_photons;	
  for(il=0; il<=nl+1; il++)
    Out_Ptr->A_l[il] *= scale1;
  
  Out_Ptr->A *=scale1;
}

/***********************************************************
 *	Sum and scale results of current run.
 ****/
void SumScaleResult(InputStruct In_Parm, 
					OutStruct * Out_Ptr)
{
  /* Get 1D & 0D results. */
  Sum2DRd(In_Parm, Out_Ptr);
  Sum2DA(In_Parm,  Out_Ptr);
  Sum2DTt(In_Parm, Out_Ptr);
  
  ScaleRdTt(In_Parm, Out_Ptr);
  ScaleA(In_Parm, Out_Ptr);
}

/***********************************************************
 *	Write the version number as the first string in the 
 *	file.
 *	Use chars only so that they can be read as either 
 *	ASCII or binary.
 ****/
void WriteVersion(FILE *file, char *Version)
{
  fprintf(file, 
	"%s \t# Version number of the file format.\n\n", 
	Version);
  fprintf(file, "####\n# Data categories include: \n");
  fprintf(file, "# InParm, RAT, \n");
  fprintf(file, "# A_l, A_z, Rd_r, Rd_a, Tt_r, Tt_a, \n");
  fprintf(file, "# A_rz, Rd_ra, Tt_ra \n####\n\n");
}

/***********************************************************
 *	Write the input parameters to the file.
 ****/
void WriteInParm(FILE *file, InputStruct In_Parm)
{
  short i;
  
  fprintf(file, 
	"InParm \t\t\t# Input parameters. cm is used.\n");
  
  fprintf(file, 
	"%s \tA\t\t# output file name, ASCII.\n", 
	In_Parm.out_fname);
  fprintf(file, 
	"%ld \t\t\t# No. of photons\n", In_Parm.num_photons);
  
  fprintf(file, 
	"%G\t%G\t\t# dz, dr [cm]\n", In_Parm.dz,In_Parm.dr);
  fprintf(file, "%hd\t%hd\t%hd\t# No. of dz, dr, da.\n\n", 
	  In_Parm.nz, In_Parm.nr, In_Parm.na);
  
  fprintf(file, 
	"%hd\t\t\t\t\t# Number of layers\n", 
	In_Parm.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_Parm.layerspecs[0].n); 
  for(i=1; i<=In_Parm.num_layers; i++)  {
    LayerStruct s;
    s = In_Parm.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_Parm.layerspecs[i].n); 
}

/***********************************************************
 *	Write reflectance, absorption, transmission. 
 ****/
void WriteRAT(FILE * file, OutStruct Out_Parm)
{
  fprintf(file, 
	"RAT #Reflectance, absorption, transmission. \n");
	/* flag. */

  fprintf(file, 
	"%-14.6G \t#Specular reflectance [-]\n", Out_Parm.Rsp);
  fprintf(file, 
	"%-14.6G \t#Diffuse reflectance [-]\n", Out_Parm.Rd);
  fprintf(file, 
	"%-14.6G \t#Absorbed fraction [-]\n", Out_Parm.A);
  fprintf(file, 
	"%-14.6G \t#Transmittance [-]\n", Out_Parm.Tt);
  
  fprintf(file, "\n");
}

/***********************************************************
 *	Write absorption as a function of layer. 
 ****/
void WriteA_layer(FILE * file, 
				  short Num_Layers,
				  OutStruct Out_Parm)
{
  short i;
  
  fprintf(file, 
	"A_l #Absorption as a function of layer. [-]\n");
	/* flag. */

  for(i=1; i<=Num_Layers; i++)
    fprintf(file, "%12.4G\n", Out_Parm.A_l[i]);
  fprintf(file, "\n");
}

/***********************************************************
 *	5 numbers each line.
 ****/
void WriteRd_ra(FILE * file, 
				short Nr,
				short Na,
				OutStruct Out_Parm)
{
  short ir, ia;
  
  fprintf(file, 
	  "%s\n%s\n%s\n%s\n%s\n%s\n",	/* flag. */
	  "# Rd[r][angle]. [1/(cm2sr)].",
	  "# Rd[0][0], [0][1],..[0][na-1]",
	  "# Rd[1][0], [1][1],..[1][na-1]",
	  "# ...",
	  "# Rd[nr-1][0], [nr-1][1],..[nr-1][na-1]",
	  "Rd_ra");
  
  for(ir=0;ir<Nr;ir++)
    for(ia=0;ia<Na;ia++)  {
      fprintf(file, "%12.4E ", Out_Parm.Rd_ra[ir][ia]);
      if( (ir*Na + ia + 1)%5 == 0) fprintf(file, "\n");
    }
  
  fprintf(file, "\n");
}

/***********************************************************
 *	1 number each line.
 ****/
void WriteRd_r(FILE * file, 
			   short Nr,
			   OutStruct Out_Parm)
{
  short ir;
  
  fprintf(file, 
	"Rd_r #Rd[0], [1],..Rd[nr-1]. [1/cm2]\n");	/* flag. */
  
  for(ir=0;ir<Nr;ir++) {
    fprintf(file, "%12.4E\n", Out_Parm.Rd_r[ir]);
  }
  
  fprintf(file, "\n");
}

/***********************************************************
 *	1 number each line.
 ****/
void WriteRd_a(FILE * file, 
			   short Na,
			   OutStruct Out_Parm)
{
  short ia;
  
  fprintf(file, 
	"Rd_a #Rd[0], [1],..Rd[na-1]. [sr-1]\n");	/* flag. */
  
  for(ia=0;ia<Na;ia++) {
    fprintf(file, "%12.4E\n", Out_Parm.Rd_a[ia]);
  }
  
  fprintf(file, "\n");
}

/***********************************************************
 *	5 numbers each line.
 ****/
void WriteTt_ra(FILE * file, 
				short Nr,
				short Na,
				OutStruct Out_Parm)
{
  short ir, ia;
  
  fprintf(file, 
	  "%s\n%s\n%s\n%s\n%s\n%s\n",	/* flag. */
	  "# Tt[r][angle]. [1/(cm2sr)].",
	  "# Tt[0][0], [0][1],..[0][na-1]",
	  "# Tt[1][0], [1][1],..[1][na-1]",
	  "# ...",
	  "# Tt[nr-1][0], [nr-1][1],..[nr-1][na-1]",
	  "Tt_ra");
  
  for(ir=0;ir<Nr;ir++)
    for(ia=0;ia<Na;ia++)  {
      fprintf(file, "%12.4E ", Out_Parm.Tt_ra[ir][ia]);
      if( (ir*Na + ia + 1)%5 == 0) fprintf(file, "\n");
    }
  
  fprintf(file, "\n");
}

/***********************************************************
 *	5 numbers each line.
 ****/
void WriteA_rz(FILE * file, 
			   short Nr,
			   short Nz,
			   OutStruct Out_Parm)
{
  short iz, ir;
  
  fprintf(file, 
	  "%s\n%s\n%s\n%s\n%s\n%s\n", /* flag. */
	  "# A[r][z]. [1/cm3]",
	  "# A[0][0], [0][1],..[0][nz-1]",
	  "# A[1][0], [1][1],..[1][nz-1]",
	  "# ...",
	  "# A[nr-1][0], [nr-1][1],..[nr-1][nz-1]",
	  "A_rz");
  
  for(ir=0;ir<Nr;ir++)
    for(iz=0;iz<Nz;iz++)  {
      fprintf(file, "%12.4E ", Out_Parm.A_rz[ir][iz]);
      if( (ir*Nz + iz + 1)%5 == 0) fprintf(file, "\n");
    }
  
  fprintf(file, "\n");
}

/***********************************************************
 *	1 number each line.
 ****/
void WriteA_z(FILE * file, 
			  short Nz,
			  OutStruct Out_Parm)
{
  short iz;
  
  fprintf(file, 
	"A_z #A[0], [1],..A[nz-1]. [1/cm]\n");	/* flag. */
  
  for(iz=0;iz<Nz;iz++) {
    fprintf(file, "%12.4E\n", Out_Parm.A_z[iz]);
  }
  
  fprintf(file, "\n");
}

/***********************************************************
 *	1 number each line.
 ****/
void WriteTt_r(FILE * file, 
			   short Nr,
			   OutStruct Out_Parm)
{
  short ir;
  
  fprintf(file, 
	"Tt_r #Tt[0], [1],..Tt[nr-1]. [1/cm2]\n"); /* flag. */
  
  for(ir=0;ir<Nr;ir++) {
    fprintf(file, "%12.4E\n", Out_Parm.Tt_r[ir]);
  }
  
  fprintf(file, "\n");
}

/***********************************************************
 *	1 number each line.
 ****/
void WriteTt_a(FILE * file, 
			   short Na,
			   OutStruct Out_Parm)
{
  short ia;
  
  fprintf(file, 
	"Tt_a #Tt[0], [1],..Tt[na-1]. [sr-1]\n"); /* flag. */
  
  for(ia=0;ia<Na;ia++) {
    fprintf(file, "%12.4E\n", Out_Parm.Tt_a[ia]);
  }
  
  fprintf(file, "\n");
}

/***********************************************************
 ****/
void WriteResult(InputStruct In_Parm, 
				 OutStruct Out_Parm,
				 char * TimeReport)
{
  FILE *file;
  
  file = fopen(In_Parm.out_fname, "w");
  if(file == NULL) nrerror("Cannot open file to write.\n");
  
  if(toupper(In_Parm.out_fformat) == 'A') 
	WriteVersion(file, "A1");
  else 
	WriteVersion(file, "B1");

  fprintf(file, "# %s", TimeReport);
  fprintf(file, "\n");
  
  WriteInParm(file, In_Parm);
  WriteRAT(file, Out_Parm);	
	/* reflectance, absorption, transmittance. */
  
  /* 1D arrays. */
  WriteA_layer(file, In_Parm.num_layers, Out_Parm);
  WriteA_z(file, In_Parm.nz, Out_Parm);
  WriteRd_r(file, In_Parm.nr, Out_Parm);
  WriteRd_a(file, In_Parm.na, Out_Parm);
  WriteTt_r(file, In_Parm.nr, Out_Parm);
  WriteTt_a(file, In_Parm.na, Out_Parm);
  
  /* 2D arrays. */
  WriteA_rz(file, In_Parm.nr, In_Parm.nz, Out_Parm);
  WriteRd_ra(file, In_Parm.nr, In_Parm.na, Out_Parm);
  WriteTt_ra(file, In_Parm.nr, In_Parm.na, Out_Parm);
  
  fclose(file);
}