convconv.c

蒙特卡罗算法。用盟特卡罗算法模拟光在组织钟的传播

/****************************************************************
 *	Specify photon beam profile, and convolute over the beam.
 *
 ****/

#include "conv.h"

/****************************************************************
 *	If IBMPC is 1, file extensions are limited to 3 letters
 *	at most.
 ****/
#define IBMPC 0

/****************************************************************
 *	Used by convolution over Gaussian beam.  Ignore the value
 *	beyond GAUSSLIMIT radius.
 ****/
#define GAUSSLIMIT 4


/***********************Declarations****************************/
FILE       *GetWriteFile(char *);
void        AllocConvData(InputStruct *, OutStruct *);
void        FreeConvData(InputStruct *, OutStruct *);
float       qtrap(float (*) (float), float, float, float);
double      BessI0(double);

/****************************************************************
 *	Data structures for the binary tree used to store part of
 *	the integrand evaluation.
 ****/
struct Node {
  float       x, y;
  struct Node *left, *right;
};

typedef struct Node *LINK;
typedef LINK TREE;

/****************************************************************
 *	A global structure to pass the current coordinate of the
 * 	physical quantities being evaluated and the pointers of the
 *	input and output parameters to the integration function.
 ****/
struct {
  double      r;
  short       iz, ia;
  InputStruct *in_ptr;
  OutStruct  *out_ptr;
  TREE        tree;		/* A tree to store ITheta() &*
				 * ExpBessI0(). */
}           ConvVar;


/***********************Convolution*****************************/
/****************************************************************
 *	Specify the parameters for a flat beam.
 ****/
void
GetFlatBeam(BeamStruct * Beam_Ptr)
{
  Beam_Ptr->type = flat;

  printf("Total energy of the flat beam [J]: ");
  Beam_Ptr->P = GetFloat(FLT_MIN, FLT_MAX);
  printf("Radius of the flat beam [cm]: ");
  Beam_Ptr->R = GetFloat(FLT_MIN, FLT_MAX);
  printf("Total power: %8.2lg J, and radius: %8.2lg cm.\n",
	 Beam_Ptr->P, Beam_Ptr->R);
}

/****************************************************************
 *	Specify the parameters for a Gaussian beam.
 ****/
void
GetGaussianBeam(BeamStruct * Beam_Ptr)
{
  Beam_Ptr->type = gaussian;

  printf("Total energy of the Gaussian beam [J]: ");
  Beam_Ptr->P = GetFloat(FLT_MIN, FLT_MAX);
  printf("1/e2 Radius of the Gaussian beam [cm]: ");
  Beam_Ptr->R = GetFloat(FLT_MIN, FLT_MAX);
  printf("Total power: %8.2lg J, and radius: %8.2lg cm.\n",
	 Beam_Ptr->P, Beam_Ptr->R);
}

/****************************************************************
 *	Specify the beam profile.
 ****/
void
LaserBeam(BeamStruct * Beam_Ptr, OutStruct * Out_Ptr)
{
  ConvStruct  null_conved = NULLCONVSTRUCT;
  char        cmd_str[STRLEN];

  printf("Beam profile:f=flat, g=Gaussian. q=quit: ");
  do
    gets(cmd_str);
  while (!strlen(cmd_str));

  switch (toupper(cmd_str[0])) {
  case 'F':
    Out_Ptr->conved = null_conved;
    GetFlatBeam(Beam_Ptr);
    break;
  case 'G':
    Out_Ptr->conved = null_conved;
    GetGaussianBeam(Beam_Ptr);
    break;
  default:;
  }
}


/****************************************************************
 *	Convolution.
 ***************************************************************/

/****************************************************************
 *	Specify the resolution and the number of points in r
 *	direction.
 *	Set the Out_Ptr->conved to null.
 *	Reallocate the arrays for the convolution arrays.
 ****/
void
ConvResolution(InputStruct * In_Ptr, OutStruct * Out_Ptr)
{
  char        in_str[STRLEN];
  ConvStruct  null_conved = NULLCONVSTRUCT;

  if (!Out_Ptr->allocated) {
    puts("...No mcml output data to work with");
    return;
  }
  printf("Current resolution: %8.2lg cm and number of points: %hd\n",
	 In_Ptr->drc, In_Ptr->nrc);
  printf("Input resolution in r direction [cm]: ");
  In_Ptr->drc = GetFloat(FLT_MIN, FLT_MAX);
  printf("Input number of points in r direction: ");
  In_Ptr->nrc = GetShort(1, SHRT_MAX);

  printf("Resolution: %8.2lg cm and number of points: %hd\n",
	 In_Ptr->drc, In_Ptr->nrc);

  Out_Ptr->conved = null_conved;
  FreeConvData(In_Ptr, Out_Ptr);
  AllocConvData(In_Ptr, Out_Ptr);
}

/****************************************************************
 *	Specify the relative convolution error. 0.001-0.1 recomended.
 *	Set the Out_Ptr->conved to null.
 ****/
void
ConvError(InputStruct * In_Ptr, OutStruct * Out_Ptr)
{
  char        in_str[STRLEN];
  ConvStruct  null_conved = NULLCONVSTRUCT;
  float       eps;

  printf("Relative convolution error\n");
  printf("Current value is %8.2g (0.001-0.1 recommended): ",
	 In_Ptr->eps);
  In_Ptr->eps = GetFloat(FLT_MIN, 1);
  Out_Ptr->conved = null_conved;
}

/***********************Binary Tree*****************************/
LINK
FillNode(float x, float y)
{
  LINK        l;

  l = (LINK) malloc(sizeof(struct Node));
  if (l) {
    l->x = x;
    l->y = y;
    l->left = l->right = NULL;
  }
  return (l);
}

/****************************************************************
 *	Assume the (x, y) is not in the tree.
 ****/
void
InsertNode(TREE * TreePtr, float x, float y)
{
  LINK        l1, l2;		/* l1 buffers l2. */

  l1 = NULL;
  l2 = *TreePtr;
  while (l2 != NULL) {
    l1 = l2;
    if (x < l2->x)
      l2 = l2->left;
    else
      l2 = l2->right;
  }

  if (l1 == NULL)		/* Empty tree. */
    *TreePtr = FillNode(x, y);
  else if (x < l1->x)
    l1->left = FillNode(x, y);
  else
    l1->right = FillNode(x, y);
}

LINK
SearchNode(TREE Tree, float x)
{
  LINK        l = Tree;
  Boolean     found = 0;

  while (l != NULL && !found)
    if (x < l->x)
      l = l->left;
    else if (x > l->x)
      l = l->right;
    else
      found = 1;

  return (l);
}

void
FreeTree(TREE Tree)
{
  if (Tree) {
    FreeTree(Tree->left);
    FreeTree(Tree->right);
    free(Tree);
  }
}

/***********************Integration*****************************/
/****************************************************************
 *	Compute Itheta shown in the manual.
 ****/
double
ITheta(double r, double r2, double R)
{
  double      temp;

  if (R >= r + r2)
    temp = 1;
  else if (fabs(r - r2) <= R) {
    temp = (r * r + r2 * r2 - R * R) / (2 * r * r2);
    if (fabs(temp) > 1)
      temp = SIGN(temp);
    temp = acos(temp) / PI;
  } else			/* R < fabs(r-r2) */
    temp = 0;

  return (temp);
}

/****************************************************************
 ****/
double
ExpBessI0(double r, double r2, double R)
{
  double      expbess;
  double      _RR = 1 / (R * R);
  double      x = 4 * r * r2 * _RR;
  double      y = 2 * (r2 * r2 + r * r) * _RR;

  expbess = exp(-y + x) * BessI0(x);
  return (expbess);
}

/****************************************************************
 *	Interpolate for the arrays A_rz[].
 ****/
double
A_rzInterp(double **A_rz, double r2)
{
  double      ir2, A_lo, A_hi, A_at_r2;
  short       iz, ir2lo, nr = ConvVar.in_ptr->nr;

  ir2 = r2 / ConvVar.in_ptr->dr;
  iz = ConvVar.iz;
  if (nr < 3)
    A_at_r2 = A_rz[0][iz];
  else if (ir2 < nr - 1.5) {	/* interpolation. */
    ir2lo = MAX(0, (short) (ir2 - 0.5));	/* truncation. */
    A_lo = A_rz[ir2lo][iz];
    A_hi = A_rz[ir2lo + 1][iz];
    A_at_r2 = A_lo + (A_hi - A_lo) * (ir2 - ir2lo - 0.5);
  } else {			/* extrapolation. */
    ir2lo = nr - 3;
    A_lo = A_rz[ir2lo][iz];
    A_hi = A_rz[ir2lo + 1][iz];
    if (A_lo >= A_hi)		/* Noise test. */
      A_at_r2 = A_lo + (A_hi - A_lo) * (ir2 - ir2lo - 0.5);
    else
      A_at_r2 = 0.0;
  }

  return (MAX(0, A_at_r2));
}

/****************************************************************
 *	Interpolate for the arrays Rd_ra[] or Tt_ra[].
 ****/
double
RT_raInterp(double **RT_ra, double r2)
{
  double      ir2, RT_lo, RT_hi, RT_at_r2;
  short       ia, ir2lo, nr = ConvVar.in_ptr->nr;

  ir2 = r2 / ConvVar.in_ptr->dr;
  ia = ConvVar.ia;
  if (nr < 3)
    RT_at_r2 = RT_ra[0][ia];
  else if (ir2 < nr - 1.5) {	/* interpolation. */
    ir2lo = MAX(0, (short) (ir2 - 0.5));	/* truncation. */
    RT_lo = RT_ra[ir2lo][ia];
    RT_hi = RT_ra[ir2lo + 1][ia];
    RT_at_r2 = RT_lo + (RT_hi - RT_lo) * (ir2 - ir2lo - 0.5);
  } else {			/* extrapolation. */
    ir2lo = nr - 3;
    RT_lo = RT_ra[ir2lo][ia];
    RT_hi = RT_ra[ir2lo + 1][ia];
    if (RT_lo >= RT_hi)		/* Noise test. */
      RT_at_r2 = RT_lo + (RT_hi - RT_lo) * (ir2 - ir2lo - 0.5);
    else
      RT_at_r2 = 0.0;
  }

  return (MAX(0, RT_at_r2));
}

/****************************************************************
 *	Interpolate for the arrays Rd_r[] or Tt_r[].
 ****/
double
RT_rInterp(double *RT_r, double r2)
{
  double      ir2, RT_lo, RT_hi, RT_at_r2;
  short       ir2lo, nr = ConvVar.in_ptr->nr;

  ir2 = r2 / ConvVar.in_ptr->dr;
  if (nr < 3)
    RT_at_r2 = RT_r[0];
  else if (ir2 < nr - 1.5) {	/* interpolation. */
    ir2lo = MAX(0, (short) (ir2 - 0.5));	/* truncation. */
    RT_lo = RT_r[ir2lo];
    RT_hi = RT_r[ir2lo + 1];
    RT_at_r2 = RT_lo + (RT_hi - RT_lo) * (ir2 - ir2lo - 0.5);
  } else {			/* extrapolation. */
    ir2lo = nr - 3;
    RT_lo = RT_r[ir2lo];
    RT_hi = RT_r[ir2lo + 1];
    if (RT_lo >= RT_hi)		/* Noise test. */
      RT_at_r2 = RT_lo + (RT_hi - RT_lo) * (ir2 - ir2lo - 0.5);
    else
      RT_at_r2 = 0.0;
  }

  return (MAX(0, RT_at_r2));
}

/****************************************************************
 *	Convolution integrand for either flat or gaussian beams.
 *	Return the integrand for the convolution integral.
 *	r2 is the r" in the formula shown in the manual.
 *	When r2 is in the range of recorded array, interpolation
 *	is used to evaluate the diffuse reflectance at r2.
 *	Note that since the last grid elements collect all the
 *	photon weight that falls beyond the grid system, we should
 *	avoid using them in the convolution.
 ****/
float
A_rzFGIntegrand(float r2)
{				/* r" in the integration. */
  float       f;
  short       nr = ConvVar.in_ptr->nr;
  double      R, r, A_at_r2;
  LINK        link;

  A_at_r2 = A_rzInterp(ConvVar.out_ptr->A_rz, r2);

  R = ConvVar.in_ptr->beam.R;
  r = ConvVar.r;
  if ((link = SearchNode(ConvVar.tree, r2)))	/* f in tree. */
    f = link->y;
  else {
    if (ConvVar.in_ptr->beam.type == flat)
      f = ITheta(r, r2, R);
    else			/* Gaussian. */
      f = ExpBessI0(r, r2, R);
    InsertNode(&ConvVar.tree, r2, f);
  }

  f *= A_at_r2 * r2;
  return (f);
}

/****************************************************************
 *	Convolution integrand for either flat or gaussian beams.
 *	See comments for A_rzFGIntegrand().
 ****/
float
Rd_raFGIntegrand(float r2)
{				/* r" in the integration. */
  float       f;
  short       nr = ConvVar.in_ptr->nr;
  double      R, r, Rd_at_r2;
  LINK        link;

  Rd_at_r2 = RT_raInterp(ConvVar.out_ptr->Rd_ra, r2);

  R = ConvVar.in_ptr->beam.R;
  r = ConvVar.r;
  if ((link = SearchNode(ConvVar.tree, r2)))	/* f in tree. */
    f = link->y;
  else {
    if (ConvVar.in_ptr->beam.type == flat)
      f = ITheta(r, r2, R);
    else			/* Gaussian. */
      f = ExpBessI0(r, r2, R);
    InsertNode(&ConvVar.tree, r2, f);
  }

  f *= Rd_at_r2 * r2;
  return (f);
}

/****************************************************************
 *	Convolution integrand for either flat or gaussian beams.
 *	See comments for A_rzFGIntegrand().
 ****/
float
Rd_rFGIntegrand(float r2)
{				/* r" in the integration. */
  float       f;
  short       nr = ConvVar.in_ptr->nr;
  double      R, r, Rd_at_r2;

  Rd_at_r2 = RT_rInterp(ConvVar.out_ptr->Rd_r, r2);

  R = ConvVar.in_ptr->beam.R;
  r = ConvVar.r;
  if (ConvVar.in_ptr->beam.type == flat)
    f = Rd_at_r2 * ITheta(r, r2, R) * r2;
  else				/* Gaussian. */
    f = Rd_at_r2 * ExpBessI0(r, r2, R) * r2;

  return (f);
}

/****************************************************************
 *	Convolution integrand for either flat or gaussian beams.
 *	See comments for A_rzFGIntegrand().
 ****/
float
Tt_raFGIntegrand(float r2)
{				/* r" in the integration. */
  float       f;
  short       nr = ConvVar.in_ptr->nr;
  double      R, r, Tt_at_r2;
  LINK        link;

  Tt_at_r2 = RT_raInterp(ConvVar.out_ptr->Tt_ra, r2);

  R = ConvVar.in_ptr->beam.R;
  r = ConvVar.r;
  if ((link = SearchNode(ConvVar.tree, r2)))	/* f in tree. */
    f = link->y;
  else {
    if (ConvVar.in_ptr->beam.type == flat)
      f = ITheta(r, r2, R);
    else			/* Gaussian. */
      f = ExpBessI0(r, r2, R);
    InsertNode(&ConvVar.tree, r2, f);
  }

  f *= Tt_at_r2 * r2;
  return (f);
}

/****************************************************************
 *	Convolution integrand for either flat or gaussian beams.
 *	See comments for A_rzFGIntegrand().
 ****/
float
Tt_rFGIntegrand(float r2)
{				/* r" in the integration. */
  float       f;
  short       nr = ConvVar.in_ptr->nr;
  double      R, r, Tt_at_r2;

  Tt_at_r2 = RT_rInterp(ConvVar.out_ptr->Tt_r, r2);

  R = ConvVar.in_ptr->beam.R;
  r = ConvVar.r;
  if (ConvVar.in_ptr->beam.type == flat)
    f = Tt_at_r2 * ITheta(r, r2, R) * r2;
  else				/* Gaussian. */
    f = Tt_at_r2 * ExpBessI0(r, r2, R) * r2;

  return (f);
}

/****************************************************************
 ****/