make_poly.cpp 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122
  1. /* make_poly.cpp */
  2. /* Creates a poly-energetic kernel given the energy-binned beam fluence and
  3. kernels for all of the beam energies. */
  4. #include "defs.h"
  5. //fillers for these entries in kernel structure
  6. #define UNCERT 1.0
  7. #define MEAN_RADIUS 0.0
  8. #define MEAN_ANGLE 0.0
  9. extern char errstr[200]; // error string that all routines have access to
  10. int make_poly_kernel(MONO_KERNELS *mono, KERNEL *poly)
  11. {
  12. // There is a problem with the first block of commented statements, likely a memory leak
  13. KERNEL_CATEGORIES category;
  14. int i, j, e;
  15. float sum;
  16. poly->nradii = N_KERNEL_RADII;
  17. poly->ntheta = N_KERNEL_ANGLES;
  18. poly->radial_boundary = (float *)malloc(poly->nradii*sizeof(float));
  19. poly->angular_boundary = (float *)malloc(poly->ntheta*sizeof(float));
  20. //copy radial boundaries from first mono kernel
  21. for (i=0;i<poly->nradii;i++)
  22. poly->radial_boundary[i] = mono->kernel[0].radial_boundary[i];
  23. //copy angular boundaries from first mono kernel
  24. for (i=0;i<poly->ntheta;i++)
  25. poly->angular_boundary[i] = mono->kernel[0].angular_boundary[i];
  26. for (i=0;i<N_KERNEL_CATEGORIES;i++)
  27. if ( (poly->matrix[i] =
  28. (float *) malloc(poly->ntheta*poly->nradii*sizeof(float))) == NULL)
  29. {
  30. sprintf(errstr,"Cannot allocate space for matrix %d\n",i);
  31. return(FAILURE);
  32. }
  33. if ( (poly->total_matrix =
  34. (float *) malloc(poly->ntheta*poly->nradii*sizeof(float))) == NULL)
  35. {
  36. sprintf(errstr,"Cannot allocate space for total matrix\n");
  37. return(FAILURE);
  38. }
  39. for (j=0;j<poly->ntheta;j++)
  40. for (i=0;i<poly->nradii;i++)
  41. {
  42. KERNEL_TOTAL_VALUE(poly,i,j) = 0.0;
  43. //weight of each mono kernel value in sum is fluence*energy*mu
  44. category = primary_;
  45. KERNEL_VALUE(poly,category,i,j) = 0.0;
  46. sum = 0.0;
  47. for (e=0;e<mono->nkernels;e++)
  48. {
  49. KERNEL_VALUE(poly,category,i,j) += mono->fluence[e]*mono->energy[e]*mono->mu[e]
  50. * KERNEL_VALUE(&(mono->kernel[e]),category,i,j);
  51. sum += mono->fluence[e]*mono->energy[e]*mono->mu[e];
  52. }
  53. KERNEL_VALUE(poly,category,i,j) /= sum;
  54. KERNEL_TOTAL_VALUE(poly,i,j) += KERNEL_VALUE(poly,category,i,j);
  55. category = first_scatter_;
  56. KERNEL_VALUE(poly,category,i,j) = 0.0;
  57. sum = 0.0;
  58. for (e=0;e<mono->nkernels;e++)
  59. {
  60. KERNEL_VALUE(poly,category,i,j) += mono->fluence[e]*mono->energy[e]*mono->mu[e]
  61. * KERNEL_VALUE(&(mono->kernel[e]),category,i,j);
  62. sum += mono->fluence[e]*mono->energy[e]*mono->mu[e];
  63. }
  64. KERNEL_VALUE(poly,category,i,j) /= sum;
  65. KERNEL_TOTAL_VALUE(poly,i,j) += KERNEL_VALUE(poly,category,i,j);
  66. category = second_scatter_;
  67. KERNEL_VALUE(poly,category,i,j) = 0.0;
  68. sum = 0.0;
  69. for (e=0;e<mono->nkernels;e++)
  70. {
  71. KERNEL_VALUE(poly,category,i,j) += mono->fluence[e]*mono->energy[e]*mono->mu[e]
  72. * KERNEL_VALUE(&(mono->kernel[e]),category,i,j);
  73. sum += mono->fluence[e]*mono->energy[e]*mono->mu[e];
  74. }
  75. KERNEL_VALUE(poly,category,i,j) /= sum;
  76. KERNEL_TOTAL_VALUE(poly,i,j) += KERNEL_VALUE(poly,category,i,j);
  77. category = multiple_scatter_;
  78. KERNEL_VALUE(poly,category,i,j) = 0.0;
  79. sum = 0.0;
  80. for (e=0;e<mono->nkernels;e++)
  81. {
  82. KERNEL_VALUE(poly,category,i,j) += mono->fluence[e]*mono->energy[e]*mono->mu[e]
  83. * KERNEL_VALUE(&(mono->kernel[e]),category,i,j);
  84. sum += mono->fluence[e]*mono->energy[e]*mono->mu[e];
  85. }
  86. KERNEL_VALUE(poly,category,i,j) /= sum;
  87. KERNEL_TOTAL_VALUE(poly,i,j) += KERNEL_VALUE(poly,category,i,j);
  88. category = brem_annih_;
  89. KERNEL_VALUE(poly,category,i,j) = 0.0;
  90. sum = 0.0;
  91. for (e=0;e<mono->nkernels;e++)
  92. {
  93. KERNEL_VALUE(poly,category,i,j) += mono->fluence[e]*mono->energy[e]*mono->mu[e]
  94. * KERNEL_VALUE(&(mono->kernel[e]),category,i,j);
  95. sum += mono->fluence[e]*mono->energy[e]*mono->mu[e];
  96. }
  97. KERNEL_VALUE(poly,category,i,j) /= sum;
  98. KERNEL_TOTAL_VALUE(poly,i,j) += KERNEL_VALUE(poly,category,i,j);
  99. }
  100. return(SUCCESS);
  101. }