Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
325 changes: 168 additions & 157 deletions mcxtrace-comps/samples/Abs_objects.comp

Large diffs are not rendered by default.

252 changes: 128 additions & 124 deletions mcxtrace-comps/samples/Absorption_sample.comp
Original file line number Diff line number Diff line change
Expand Up @@ -75,10 +75,10 @@ SHARE
%{
%include "read_table-lib"
%include "interoff-lib"

#ifndef SHAPES_T
#define SHAPES_T
enum shapes_t {NONE=-1,SPHERE, CYLINDER, CUBE, ELLIPSOID, ANY};
enum shapes_t { NONE = -1, SPHERE, CYLINDER, CUBE, ELLIPSOID, ANY };
#endif
%}

Expand Down Expand Up @@ -106,178 +106,182 @@ INITIALIZE
%{
int status;
V_i = 1;

/* Checking volumes */
if (geometry_o && strlen(geometry_o) && strcmp(geometry_o, "NULL") && strcmp(geometry_o, "0")) {
if (!off_init(geometry_o, xwidth_o, yheight_o, zdepth_o, 0, &offdata_o))
exit(printf("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n",
NAME_CURRENT_COMP, geometry_o));
if (geometry_o && strlen (geometry_o) && strcmp (geometry_o, "NULL") && strcmp (geometry_o, "0")) {
if (!off_init (geometry_o, xwidth_o, yheight_o, zdepth_o, 0, &offdata_o))
exit (printf ("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n", NAME_CURRENT_COMP, geometry_o));
shape_o = ANY;
} else if (!yheight_o && !xwidth_o && !zdepth_o && radius_o)
shape_o = SPHERE;
else if (yheight_o && radius_o)
shape_o = CYLINDER;
else if (yheight_o && !radius_o) {
if (!zdepth_o) zdepth_o=yheight_o;
if (!xwidth_o) xwidth_o=yheight_o;
if (!zdepth_o)
zdepth_o = yheight_o;
if (!xwidth_o)
xwidth_o = yheight_o;
shape_o = CUBE;
} else {
exit(fprintf(stderr,"ERROR: %s: Unmeaningful outer volume description.\n",NAME_CURRENT_COMP));
exit (fprintf (stderr, "ERROR: %s: Unmeaningful outer volume description.\n", NAME_CURRENT_COMP));
}


if (geometry_i && strlen(geometry_i) && strcmp(geometry_i, "NULL") && strcmp(geometry_i, "0")) {
if (!off_init(geometry_i, xwidth_i, yheight_i, zdepth_i, 0, &offdata_i))
exit(printf("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n",
NAME_CURRENT_COMP, geometry_i));

if (geometry_i && strlen (geometry_i) && strcmp (geometry_i, "NULL") && strcmp (geometry_i, "0")) {
if (!off_init (geometry_i, xwidth_i, yheight_i, zdepth_i, 0, &offdata_i))
exit (printf ("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n", NAME_CURRENT_COMP, geometry_i));
shape_i = ANY;
} else if (!yheight_i && !xwidth_i && !zdepth_i && radius_i)
shape_i = SPHERE;
else if (yheight_i && radius_i)
shape_i = CYLINDER;
else if (yheight_i && !radius_i) {
if (!zdepth_i) zdepth_i=yheight_i;
if (!xwidth_i) xwidth_i=yheight_i;
if (!zdepth_i)
zdepth_i = yheight_i;
if (!xwidth_i)
xwidth_i = yheight_i;
shape_i = CUBE;
} else {
fprintf(stderr,"Warning: %s: Invalid inner volume specification. Ignoring.\n",NAME_CURRENT_COMP);
V_i=0;
shape_i=NONE;
fprintf (stderr, "Warning: %s: Invalid inner volume specification. Ignoring.\n", NAME_CURRENT_COMP);
V_i = 0;
shape_i = NONE;
}

/* Loading datafiles */
if ( material_datafile_o && (status=Table_Read(&t_o,material_datafile_o,0))==-1){
fprintf(stderr,"Error: (%s) Could not parse file \"%s\".\n",NAME_CURRENT_COMP,material_datafile_o);
exit(-1);
if (material_datafile_o && (status = Table_Read (&t_o, material_datafile_o, 0)) == -1) {
fprintf (stderr, "Error: (%s) Could not parse file \"%s\".\n", NAME_CURRENT_COMP, material_datafile_o);
exit (-1);
}
if (t_o.columns==3) { /*which column is the energy in and which holds mu*/
E_c_o=0;mu_c_o=1;
}else{
E_c_o=0;mu_c_o=5;
if (t_o.columns == 3) { /*which column is the energy in and which holds mu*/
E_c_o = 0;
mu_c_o = 1;
} else {
E_c_o = 0;
mu_c_o = 5;
}
if(rho_o==0){
char **header_parsed;
header_parsed=Table_ParseHeader(t_o.header,"Z","A[r]","rho","Z/A",NULL);
if(header_parsed[3]){ /*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
Z_A_o=strtod(header_parsed[3],NULL);
}else if ( (strlen(header_parsed[0])) && (strlen(header_parsed[1])) ){
Z_A_o=strtod(header_parsed[0],NULL)/strtod(header_parsed[1],NULL);
if (rho_o == 0) {
char** header_parsed;
header_parsed = Table_ParseHeader (t_o.header, "Z", "A[r]", "rho", "Z/A", NULL);
if (header_parsed[3]) { /*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
Z_A_o = strtod (header_parsed[3], NULL);
} else if ((strlen (header_parsed[0])) && (strlen (header_parsed[1]))) {
Z_A_o = strtod (header_parsed[0], NULL) / strtod (header_parsed[1], NULL);
}
if(strlen(header_parsed[2])){
rho_o=strtod(header_parsed[2],NULL);
printf("INFO: %s: Setting rho_o=%g from %s\n", NAME_CURRENT_COMP, rho_o, material_datafile_o);
if (strlen (header_parsed[2])) {
rho_o = strtod (header_parsed[2], NULL);
printf ("INFO: %s: Setting rho_o=%g from %s\n", NAME_CURRENT_COMP, rho_o, material_datafile_o);
}
}
if (V_i){ /*if volume is zero - don't bother to read the file*/
if ( material_datafile_i && (status=Table_Read(&t_i,material_datafile_i,0))==-1){
fprintf(stderr,"Error: Could not parse file \"%s\" in COMP %s\n",material_datafile_i,NAME_CURRENT_COMP);
exit(-1);
if (V_i) { /*if volume is zero - don't bother to read the file*/
if (material_datafile_i && (status = Table_Read (&t_i, material_datafile_i, 0)) == -1) {
fprintf (stderr, "Error: Could not parse file \"%s\" in COMP %s\n", material_datafile_i, NAME_CURRENT_COMP);
exit (-1);
}
if (t_i.columns==3) {
E_c_i=0;mu_c_i=1;
}else{
E_c_i=0;mu_c_i=5;
if (t_i.columns == 3) {
E_c_i = 0;
mu_c_i = 1;
} else {
E_c_i = 0;
mu_c_i = 5;
}
if(rho_i==0){ /*when density is not from input, try to read from tables */
char **header_parsed;
header_parsed=Table_ParseHeader(t_i.header,"Z","A[r]","rho","Z/A",NULL);
if(header_parsed[3]){/*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
Z_A_i=strtod(header_parsed[3],NULL);
}else if ( (strlen(header_parsed[0])) && (strlen(header_parsed[1])) ){
Z_A_i=strtod(header_parsed[0],NULL)/strtod(header_parsed[1],NULL);
if (rho_i == 0) { /*when density is not from input, try to read from tables */
char** header_parsed;
header_parsed = Table_ParseHeader (t_i.header, "Z", "A[r]", "rho", "Z/A", NULL);
if (header_parsed[3]) { /*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
Z_A_i = strtod (header_parsed[3], NULL);
} else if ((strlen (header_parsed[0])) && (strlen (header_parsed[1]))) {
Z_A_i = strtod (header_parsed[0], NULL) / strtod (header_parsed[1], NULL);
}
if(strlen(header_parsed[2])){
rho_i=strtod(header_parsed[2],NULL);
printf("INFO: %s: Setting rho_i=%g from %s\n", NAME_CURRENT_COMP, rho_i, material_datafile_i);
if (strlen (header_parsed[2])) {
rho_i = strtod (header_parsed[2], NULL);
printf ("INFO: %s: Setting rho_i=%g from %s\n", NAME_CURRENT_COMP, rho_i, material_datafile_i);
}
}
}
%}

TRACE
%{
double l0o=0,l1o=0,l0i=0,l1i=0,mu_o=0,mu_i=0;
int status,status_i;
if (shape_o==ANY) {
status= off_x_intersect (&l0o, &l1o, NULL, NULL, x,y,z, kx,ky,kz, offdata_o);
} else if (shape_o==CYLINDER) {
status=cylinder_intersect (&l0o,&l1o,x,y,z,kx,ky,kz,radius_o,yheight_o);
} else if (shape_o==CUBE) {
status=box_intersect (&l0o,&l1o,x,y,z,kx,ky,kz,xwidth_o,yheight_o,zdepth_o);
} else if (shape_o==SPHERE) {
status=sphere_intersect (&l0o,&l1o,x,y,z,kx,ky,kz,radius_o);
double l0o = 0, l1o = 0, l0i = 0, l1i = 0, mu_o = 0, mu_i = 0;
int status, status_i;

if (shape_o == ANY) {
status = off_x_intersect (&l0o, &l1o, NULL, NULL, x, y, z, kx, ky, kz, offdata_o);
} else if (shape_o == CYLINDER) {
status = cylinder_intersect (&l0o, &l1o, x, y, z, kx, ky, kz, radius_o, yheight_o);
} else if (shape_o == CUBE) {
status = box_intersect (&l0o, &l1o, x, y, z, kx, ky, kz, xwidth_o, yheight_o, zdepth_o);
} else if (shape_o == SPHERE) {
status = sphere_intersect (&l0o, &l1o, x, y, z, kx, ky, kz, radius_o);
}
if (status ){ /*rays intersects the enclosing material*/
PROP_DL(l0o);

if (status) { /*rays intersects the enclosing material*/
PROP_DL (l0o);
SCATTER;
status_i=0;
double k=sqrt(scalar_prod(kx,ky,kz,kx,ky,kz));
if (shape_i==ANY) {
status_i= off_x_intersect (&l0i, &l1i, NULL, NULL, (x-x_i),(y-y_i),(z-z_i), kx,ky,kz, offdata_i);
} else if (shape_i==CYLINDER) {
status_i=cylinder_intersect (&l0i,&l1i,(x-x_i),(y-y_i),(z-z_i),kx,ky,kz,radius_i,yheight_i);
} else if (shape_i==CUBE) {
status_i=box_intersect (&l0i,&l1i,(x-x_i),(y-y_i),(z-z_i),kx,ky,kz,xwidth_i,yheight_i,zdepth_i);
} else if (shape_i==SPHERE) {
status_i=sphere_intersect (&l0i,&l1i,(x-x_i),(y-y_i),(z-z_i),kx,ky,kz,radius_i);
status_i = 0;
double k = sqrt (scalar_prod (kx, ky, kz, kx, ky, kz));

if (shape_i == ANY) {
status_i = off_x_intersect (&l0i, &l1i, NULL, NULL, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, offdata_i);
} else if (shape_i == CYLINDER) {
status_i = cylinder_intersect (&l0i, &l1i, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, radius_i, yheight_i);
} else if (shape_i == CUBE) {
status_i = box_intersect (&l0i, &l1i, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, xwidth_i, yheight_i, zdepth_i);
} else if (shape_i == SPHERE) {
status_i = sphere_intersect (&l0i, &l1i, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, radius_i);
}
if(status_i){ /*rays intersect the inclusion*/
PROP_DL(l0i);

if (status_i) { /*rays intersect the inclusion*/
PROP_DL (l0i);
SCATTER;
PROP_DL(l1i-l0i);
PROP_DL (l1i - l0i);
SCATTER;
/*now calculate the mu*/
mu_i=Table_Value(t_i,k*K2E,mu_c_i)*rho_i;
p*=exp(-(l1i-l0i)*mu_i*1e2); /* 1e2 to have in m unit */
mu_i = Table_Value (t_i, k * K2E, mu_c_i) * rho_i;
p *= exp (-(l1i - l0i) * mu_i * 1e2); /* 1e2 to have in m unit */
}
PROP_DL(l1o-l0o-l1i);
PROP_DL (l1o - l0o - l1i);
SCATTER;
mu_o=Table_Value(t_o,k*K2E,mu_c_o)*rho_o;
p*=exp(-(l1o-l0o-(l1i-l0i))*mu_o*1e2);
mu_o = Table_Value (t_o, k * K2E, mu_c_o) * rho_o;
p *= exp (-(l1o - l0o - (l1i - l0i)) * mu_o * 1e2);
}
%}

MCDISPLAY
%{
if (shape_o==CUBE){
box(0,0,0,xwidth_o,yheight_o,zdepth_o,0, 0, 1, 0);
} else if (shape_o==CYLINDER){
circle("xy",0, yheight_o/2.0,0,radius_o);
circle("xy",0,-yheight_o/2.0,0,radius_o);
line( radius_o,yheight_o/2.0,0, radius_o,-yheight_o/2.0,0);
line(-radius_o,yheight_o/2.0,0,-radius_o,-yheight_o/2.0,0);
line(0,yheight_o/2.0, radius_o, 0,-yheight_o/2.0, radius_o);
line(0,yheight_o/2.0,-radius_o, 0,-yheight_o/2.0,-radius_o);
} else if (shape_o==SPHERE){
circle("xy",0,0,0,radius_o);
circle("xz",0,0,0,radius_o);
circle("yz",0,0,0,radius_o);
} else if (shape_o == ANY) { /* OFF file */
off_display(offdata_o);

if (shape_o == CUBE) {
box (0, 0, 0, xwidth_o, yheight_o, zdepth_o, 0, 0, 1, 0);
} else if (shape_o == CYLINDER) {
circle ("xy", 0, yheight_o / 2.0, 0, radius_o);
circle ("xy", 0, -yheight_o / 2.0, 0, radius_o);
line (radius_o, yheight_o / 2.0, 0, radius_o, -yheight_o / 2.0, 0);
line (-radius_o, yheight_o / 2.0, 0, -radius_o, -yheight_o / 2.0, 0);
line (0, yheight_o / 2.0, radius_o, 0, -yheight_o / 2.0, radius_o);
line (0, yheight_o / 2.0, -radius_o, 0, -yheight_o / 2.0, -radius_o);
} else if (shape_o == SPHERE) {
circle ("xy", 0, 0, 0, radius_o);
circle ("xz", 0, 0, 0, radius_o);
circle ("yz", 0, 0, 0, radius_o);
} else if (shape_o == ANY) { /* OFF file */
off_display (offdata_o);
}
if (shape_i==CUBE){
box(0,0,0,xwidth_i,yheight_i,zdepth_i,0, 0, 1, 0);
} else if (shape_i==CYLINDER){
circle("xy",0, yheight_i/2.0,0,radius_i);
circle("xy",0,-yheight_i/2.0,0,radius_i);
line( radius_i,yheight_i/2.0,0, radius_i,-yheight_i/2.0,0);
line(-radius_i,yheight_i/2.0,0,-radius_i,-yheight_i/2.0,0);
line(0,yheight_i/2.0, radius_i, 0,-yheight_i/2.0, radius_i);
line(0,yheight_i/2.0,-radius_i, 0,-yheight_i/2.0,-radius_i);
} else if (shape_i==SPHERE){
circle("xy",0,0,0,radius_i);
circle("xz",0,0,0,radius_i);
circle("yz",0,0,0,radius_i);
} else if (shape_i == ANY) { /* OFF file */
off_display(offdata_i);

if (shape_i == CUBE) {
box (0, 0, 0, xwidth_i, yheight_i, zdepth_i, 0, 0, 1, 0);
} else if (shape_i == CYLINDER) {
circle ("xy", 0, yheight_i / 2.0, 0, radius_i);
circle ("xy", 0, -yheight_i / 2.0, 0, radius_i);
line (radius_i, yheight_i / 2.0, 0, radius_i, -yheight_i / 2.0, 0);
line (-radius_i, yheight_i / 2.0, 0, -radius_i, -yheight_i / 2.0, 0);
line (0, yheight_i / 2.0, radius_i, 0, -yheight_i / 2.0, radius_i);
line (0, yheight_i / 2.0, -radius_i, 0, -yheight_i / 2.0, -radius_i);
} else if (shape_i == SPHERE) {
circle ("xy", 0, 0, 0, radius_i);
circle ("xz", 0, 0, 0, radius_i);
circle ("yz", 0, 0, 0, radius_i);
} else if (shape_i == ANY) { /* OFF file */
off_display (offdata_i);
}

%}

END
Loading
Loading