@@ -75,10 +75,10 @@ SHARE
7575%{
7676 %include "read_table-lib"
7777 %include "interoff-lib"
78-
78+
7979 #ifndef SHAPES_T
8080 #define SHAPES_T
81- enum shapes_t {NONE= -1,SPHERE, CYLINDER, CUBE, ELLIPSOID, ANY};
81+ enum shapes_t { NONE = -1, SPHERE, CYLINDER, CUBE, ELLIPSOID, ANY };
8282 #endif
8383%}
8484
@@ -106,178 +106,182 @@ INITIALIZE
106106%{
107107 int status;
108108 V_i = 1;
109-
109+
110110 /* Checking volumes */
111- if (geometry_o && strlen(geometry_o) && strcmp(geometry_o, "NULL") && strcmp(geometry_o, "0")) {
112- if (!off_init(geometry_o, xwidth_o, yheight_o, zdepth_o, 0, &offdata_o))
113- exit(printf("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n",
114- NAME_CURRENT_COMP, geometry_o));
111+ if (geometry_o && strlen (geometry_o) && strcmp (geometry_o, "NULL") && strcmp (geometry_o, "0")) {
112+ if (!off_init (geometry_o, xwidth_o, yheight_o, zdepth_o, 0, &offdata_o))
113+ exit (printf ("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n", NAME_CURRENT_COMP, geometry_o));
115114 shape_o = ANY;
116115 } else if (!yheight_o && !xwidth_o && !zdepth_o && radius_o)
117116 shape_o = SPHERE;
118117 else if (yheight_o && radius_o)
119118 shape_o = CYLINDER;
120119 else if (yheight_o && !radius_o) {
121- if (!zdepth_o) zdepth_o=yheight_o;
122- if (!xwidth_o) xwidth_o=yheight_o;
120+ if (!zdepth_o)
121+ zdepth_o = yheight_o;
122+ if (!xwidth_o)
123+ xwidth_o = yheight_o;
123124 shape_o = CUBE;
124125 } else {
125- exit(fprintf(stderr,"ERROR: %s: Unmeaningful outer volume description.\n",NAME_CURRENT_COMP));
126+ exit (fprintf (stderr, "ERROR: %s: Unmeaningful outer volume description.\n", NAME_CURRENT_COMP));
126127 }
127-
128-
129- if (geometry_i && strlen(geometry_i) && strcmp(geometry_i, "NULL") && strcmp(geometry_i, "0")) {
130- if (!off_init(geometry_i, xwidth_i, yheight_i, zdepth_i, 0, &offdata_i))
131- exit(printf("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n",
132- NAME_CURRENT_COMP, geometry_i));
128+
129+ if (geometry_i && strlen (geometry_i) && strcmp (geometry_i, "NULL") && strcmp (geometry_i, "0")) {
130+ if (!off_init (geometry_i, xwidth_i, yheight_i, zdepth_i, 0, &offdata_i))
131+ exit (printf ("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n", NAME_CURRENT_COMP, geometry_i));
133132 shape_i = ANY;
134133 } else if (!yheight_i && !xwidth_i && !zdepth_i && radius_i)
135134 shape_i = SPHERE;
136135 else if (yheight_i && radius_i)
137136 shape_i = CYLINDER;
138137 else if (yheight_i && !radius_i) {
139- if (!zdepth_i) zdepth_i=yheight_i;
140- if (!xwidth_i) xwidth_i=yheight_i;
138+ if (!zdepth_i)
139+ zdepth_i = yheight_i;
140+ if (!xwidth_i)
141+ xwidth_i = yheight_i;
141142 shape_i = CUBE;
142143 } else {
143- fprintf(stderr,"Warning: %s: Invalid inner volume specification. Ignoring.\n",NAME_CURRENT_COMP);
144- V_i= 0;
145- shape_i= NONE;
144+ fprintf (stderr, "Warning: %s: Invalid inner volume specification. Ignoring.\n", NAME_CURRENT_COMP);
145+ V_i = 0;
146+ shape_i = NONE;
146147 }
147-
148+
148149 /* Loading datafiles */
149- if ( material_datafile_o && (status= Table_Read(&t_o,material_datafile_o,0))== -1){
150- fprintf(stderr,"Error: (%s) Could not parse file \"%s\".\n",NAME_CURRENT_COMP,material_datafile_o);
151- exit(-1);
150+ if (material_datafile_o && (status = Table_Read (&t_o, material_datafile_o, 0)) == -1) {
151+ fprintf (stderr, "Error: (%s) Could not parse file \"%s\".\n", NAME_CURRENT_COMP, material_datafile_o);
152+ exit (-1);
152153 }
153- if (t_o.columns==3) { /*which column is the energy in and which holds mu*/
154- E_c_o=0;mu_c_o=1;
155- }else{
156- E_c_o=0;mu_c_o=5;
154+ if (t_o.columns == 3) { /*which column is the energy in and which holds mu*/
155+ E_c_o = 0;
156+ mu_c_o = 1;
157+ } else {
158+ E_c_o = 0;
159+ mu_c_o = 5;
157160 }
158- if(rho_o==0) {
159- char ** header_parsed;
160- header_parsed= Table_ParseHeader(t_o.header,"Z","A[r]","rho","Z/A",NULL);
161- if(header_parsed[3]){ /*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
162- Z_A_o= strtod(header_parsed[3],NULL);
163- }else if ( (strlen(header_parsed[0])) && (strlen(header_parsed[1])) ) {
164- Z_A_o= strtod(header_parsed[0],NULL)/ strtod(header_parsed[1],NULL);
161+ if (rho_o == 0) {
162+ char** header_parsed;
163+ header_parsed = Table_ParseHeader (t_o.header, "Z", "A[r]", "rho", "Z/A", NULL);
164+ if (header_parsed[3]) { /*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
165+ Z_A_o = strtod (header_parsed[3], NULL);
166+ } else if ((strlen (header_parsed[0])) && (strlen (header_parsed[1]))) {
167+ Z_A_o = strtod (header_parsed[0], NULL) / strtod (header_parsed[1], NULL);
165168 }
166- if(strlen(header_parsed[2])){
167- rho_o= strtod(header_parsed[2],NULL);
168- printf("INFO: %s: Setting rho_o=%g from %s\n", NAME_CURRENT_COMP, rho_o, material_datafile_o);
169+ if (strlen (header_parsed[2])) {
170+ rho_o = strtod (header_parsed[2], NULL);
171+ printf ("INFO: %s: Setting rho_o=%g from %s\n", NAME_CURRENT_COMP, rho_o, material_datafile_o);
169172 }
170173 }
171- if (V_i){ /*if volume is zero - don't bother to read the file*/
172- if ( material_datafile_i && (status= Table_Read(&t_i,material_datafile_i,0))== -1){
173- fprintf(stderr,"Error: Could not parse file \"%s\" in COMP %s\n",material_datafile_i,NAME_CURRENT_COMP);
174- exit(-1);
174+ if (V_i) { /*if volume is zero - don't bother to read the file*/
175+ if (material_datafile_i && (status = Table_Read (&t_i, material_datafile_i, 0)) == -1) {
176+ fprintf (stderr, "Error: Could not parse file \"%s\" in COMP %s\n", material_datafile_i, NAME_CURRENT_COMP);
177+ exit (-1);
175178 }
176- if (t_i.columns==3) {
177- E_c_i=0;mu_c_i=1;
178- }else{
179- E_c_i=0;mu_c_i=5;
179+ if (t_i.columns == 3) {
180+ E_c_i = 0;
181+ mu_c_i = 1;
182+ } else {
183+ E_c_i = 0;
184+ mu_c_i = 5;
180185 }
181- if(rho_i==0) { /*when density is not from input, try to read from tables */
182- char ** header_parsed;
183- header_parsed= Table_ParseHeader(t_i.header,"Z","A[r]","rho","Z/A",NULL);
184- if(header_parsed[3]){ /*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
185- Z_A_i= strtod(header_parsed[3],NULL);
186- }else if ( (strlen(header_parsed[0])) && (strlen(header_parsed[1])) ) {
187- Z_A_i= strtod(header_parsed[0],NULL)/ strtod(header_parsed[1],NULL);
186+ if (rho_i == 0) { /*when density is not from input, try to read from tables */
187+ char** header_parsed;
188+ header_parsed = Table_ParseHeader (t_i.header, "Z", "A[r]", "rho", "Z/A", NULL);
189+ if (header_parsed[3]) { /*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
190+ Z_A_i = strtod (header_parsed[3], NULL);
191+ } else if ((strlen (header_parsed[0])) && (strlen (header_parsed[1]))) {
192+ Z_A_i = strtod (header_parsed[0], NULL) / strtod (header_parsed[1], NULL);
188193 }
189- if(strlen(header_parsed[2])){
190- rho_i= strtod(header_parsed[2],NULL);
191- printf("INFO: %s: Setting rho_i=%g from %s\n", NAME_CURRENT_COMP, rho_i, material_datafile_i);
194+ if (strlen (header_parsed[2])) {
195+ rho_i = strtod (header_parsed[2], NULL);
196+ printf ("INFO: %s: Setting rho_i=%g from %s\n", NAME_CURRENT_COMP, rho_i, material_datafile_i);
192197 }
193198 }
194199 }
195200%}
196201
197202TRACE
198203%{
199- double l0o=0, l1o=0, l0i=0, l1i=0, mu_o=0, mu_i= 0;
200- int status,status_i;
201-
202- if (shape_o== ANY) {
203- status= off_x_intersect (&l0o, &l1o, NULL, NULL, x,y, z, kx,ky,kz, offdata_o);
204- } else if (shape_o== CYLINDER) {
205- status= cylinder_intersect (&l0o,&l1o,x,y,z, kx,ky,kz,radius_o,yheight_o);
206- } else if (shape_o== CUBE) {
207- status= box_intersect (&l0o,&l1o,x,y,z, kx,ky,kz,xwidth_o,yheight_o,zdepth_o);
208- } else if (shape_o== SPHERE) {
209- status= sphere_intersect (&l0o,&l1o,x,y,z, kx,ky,kz,radius_o);
204+ double l0o = 0, l1o = 0, l0i = 0, l1i = 0, mu_o = 0, mu_i = 0;
205+ int status, status_i;
206+
207+ if (shape_o == ANY) {
208+ status = off_x_intersect (&l0o, &l1o, NULL, NULL, x, y, z, kx, ky, kz, offdata_o);
209+ } else if (shape_o == CYLINDER) {
210+ status = cylinder_intersect (&l0o, &l1o, x, y, z, kx, ky, kz, radius_o, yheight_o);
211+ } else if (shape_o == CUBE) {
212+ status = box_intersect (&l0o, &l1o, x, y, z, kx, ky, kz, xwidth_o, yheight_o, zdepth_o);
213+ } else if (shape_o == SPHERE) {
214+ status = sphere_intersect (&l0o, &l1o, x, y, z, kx, ky, kz, radius_o);
210215 }
211-
212- if (status ) { /*rays intersects the enclosing material*/
213- PROP_DL(l0o);
216+
217+ if (status) { /*rays intersects the enclosing material*/
218+ PROP_DL (l0o);
214219 SCATTER;
215- status_i= 0;
216- double k= sqrt(scalar_prod(kx,ky,kz,kx,ky,kz));
217-
218- if (shape_i== ANY) {
219- status_i= off_x_intersect (&l0i, &l1i, NULL, NULL, (x- x_i),(y- y_i),(z- z_i), kx,ky,kz, offdata_i);
220- } else if (shape_i== CYLINDER) {
221- status_i= cylinder_intersect (&l0i,&l1i,(x- x_i),(y- y_i),(z- z_i),kx,ky,kz,radius_i,yheight_i);
222- } else if (shape_i== CUBE) {
223- status_i= box_intersect (&l0i,&l1i,(x- x_i),(y- y_i),(z- z_i),kx,ky,kz,xwidth_i,yheight_i,zdepth_i);
224- } else if (shape_i== SPHERE) {
225- status_i= sphere_intersect (&l0i,&l1i,(x- x_i),(y- y_i),(z- z_i),kx,ky,kz,radius_i);
220+ status_i = 0;
221+ double k = sqrt (scalar_prod (kx, ky, kz, kx, ky, kz));
222+
223+ if (shape_i == ANY) {
224+ status_i = off_x_intersect (&l0i, &l1i, NULL, NULL, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, offdata_i);
225+ } else if (shape_i == CYLINDER) {
226+ status_i = cylinder_intersect (&l0i, &l1i, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, radius_i, yheight_i);
227+ } else if (shape_i == CUBE) {
228+ status_i = box_intersect (&l0i, &l1i, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, xwidth_i, yheight_i, zdepth_i);
229+ } else if (shape_i == SPHERE) {
230+ status_i = sphere_intersect (&l0i, &l1i, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, radius_i);
226231 }
227-
228- if(status_i){ /*rays intersect the inclusion*/
229- PROP_DL(l0i);
232+
233+ if (status_i) { /*rays intersect the inclusion*/
234+ PROP_DL (l0i);
230235 SCATTER;
231- PROP_DL(l1i- l0i);
236+ PROP_DL (l1i - l0i);
232237 SCATTER;
233238 /*now calculate the mu*/
234- mu_i= Table_Value(t_i,k* K2E,mu_c_i)* rho_i;
235- p*= exp(-(l1i- l0i)* mu_i* 1e2); /* 1e2 to have in m unit */
239+ mu_i = Table_Value (t_i, k * K2E, mu_c_i) * rho_i;
240+ p *= exp (-(l1i - l0i) * mu_i * 1e2); /* 1e2 to have in m unit */
236241 }
237- PROP_DL(l1o- l0o- l1i);
242+ PROP_DL (l1o - l0o - l1i);
238243 SCATTER;
239- mu_o= Table_Value(t_o,k* K2E,mu_c_o)* rho_o;
240- p*= exp(-(l1o- l0o- (l1i- l0i))* mu_o* 1e2);
244+ mu_o = Table_Value (t_o, k * K2E, mu_c_o) * rho_o;
245+ p *= exp (-(l1o - l0o - (l1i - l0i)) * mu_o * 1e2);
241246 }
242247%}
243248
244249MCDISPLAY
245250%{
246-
247- if (shape_o== CUBE){
248- box(0,0,0, xwidth_o,yheight_o,zdepth_o,0, 0, 1, 0);
249- } else if (shape_o== CYLINDER){
250- circle("xy",0, yheight_o/ 2.0,0, radius_o);
251- circle("xy",0, -yheight_o/ 2.0,0, radius_o);
252- line( radius_o,yheight_o/ 2.0,0, radius_o,-yheight_o/ 2.0,0);
253- line(-radius_o,yheight_o/ 2.0,0, -radius_o,-yheight_o/ 2.0,0);
254- line(0,yheight_o/ 2.0, radius_o, 0,-yheight_o/ 2.0, radius_o);
255- line(0,yheight_o/ 2.0,-radius_o, 0,-yheight_o/ 2.0,-radius_o);
256- } else if (shape_o== SPHERE){
257- circle("xy",0,0,0, radius_o);
258- circle("xz",0,0,0, radius_o);
259- circle("yz",0,0,0, radius_o);
260- } else if (shape_o == ANY) { /* OFF file */
261- off_display(offdata_o);
251+
252+ if (shape_o == CUBE) {
253+ box (0, 0, 0, xwidth_o, yheight_o, zdepth_o, 0, 0, 1, 0);
254+ } else if (shape_o == CYLINDER) {
255+ circle ("xy", 0, yheight_o / 2.0, 0, radius_o);
256+ circle ("xy", 0, -yheight_o / 2.0, 0, radius_o);
257+ line ( radius_o, yheight_o / 2.0, 0, radius_o, -yheight_o / 2.0, 0);
258+ line (-radius_o, yheight_o / 2.0, 0, -radius_o, -yheight_o / 2.0, 0);
259+ line (0, yheight_o / 2.0, radius_o, 0, -yheight_o / 2.0, radius_o);
260+ line (0, yheight_o / 2.0, -radius_o, 0, -yheight_o / 2.0, -radius_o);
261+ } else if (shape_o == SPHERE) {
262+ circle ("xy", 0, 0, 0, radius_o);
263+ circle ("xz", 0, 0, 0, radius_o);
264+ circle ("yz", 0, 0, 0, radius_o);
265+ } else if (shape_o == ANY) { /* OFF file */
266+ off_display (offdata_o);
262267 }
263-
264- if (shape_i== CUBE){
265- box(0,0,0, xwidth_i,yheight_i,zdepth_i,0, 0, 1, 0);
266- } else if (shape_i== CYLINDER){
267- circle("xy",0, yheight_i/ 2.0,0, radius_i);
268- circle("xy",0, -yheight_i/ 2.0,0, radius_i);
269- line( radius_i,yheight_i/ 2.0,0, radius_i,-yheight_i/ 2.0,0);
270- line(-radius_i,yheight_i/ 2.0,0, -radius_i,-yheight_i/ 2.0,0);
271- line(0,yheight_i/ 2.0, radius_i, 0,-yheight_i/ 2.0, radius_i);
272- line(0,yheight_i/ 2.0,-radius_i, 0,-yheight_i/ 2.0,-radius_i);
273- } else if (shape_i== SPHERE){
274- circle("xy",0,0,0, radius_i);
275- circle("xz",0,0,0, radius_i);
276- circle("yz",0,0,0, radius_i);
277- } else if (shape_i == ANY) { /* OFF file */
278- off_display(offdata_i);
268+
269+ if (shape_i == CUBE) {
270+ box (0, 0, 0, xwidth_i, yheight_i, zdepth_i, 0, 0, 1, 0);
271+ } else if (shape_i == CYLINDER) {
272+ circle ("xy", 0, yheight_i / 2.0, 0, radius_i);
273+ circle ("xy", 0, -yheight_i / 2.0, 0, radius_i);
274+ line ( radius_i, yheight_i / 2.0, 0, radius_i, -yheight_i / 2.0, 0);
275+ line (-radius_i, yheight_i / 2.0, 0, -radius_i, -yheight_i / 2.0, 0);
276+ line (0, yheight_i / 2.0, radius_i, 0, -yheight_i / 2.0, radius_i);
277+ line (0, yheight_i / 2.0, -radius_i, 0, -yheight_i / 2.0, -radius_i);
278+ } else if (shape_i == SPHERE) {
279+ circle ("xy", 0, 0, 0, radius_i);
280+ circle ("xz", 0, 0, 0, radius_i);
281+ circle ("yz", 0, 0, 0, radius_i);
282+ } else if (shape_i == ANY) { /* OFF file */
283+ off_display (offdata_i);
279284 }
280-
281285%}
282286
283287END
0 commit comments