Permalink
Fetching contributors…
Cannot retrieve contributors at this time
1244 lines (1137 sloc) 48.5 KB
/*****************************************************************************
* McXtrace, x-ray tracing package
* Copyright (C) Risoe National Laboratory, Roskilde, Denmark
*
* Component: PowderN
*
* %I
* Written by: P. Willendrup, L. Chapon, K. Lefmann, A.B.Abrahamsen, N.B.Christensen, E.M.Lauridsen.
* Date: 4.2.98
* Version: $Revision$
* Origin: McXtrace release 1.2
* Modified by: KL, KN 20.03.98 (rewrite)
* Modified by: KL, 28.09.01 (two lines)
* Modified by: KL, 22.05.03 (background)
* Modified by: KL, PW 01.05.05 (N lines)
* Modified by: PW, LC 04.10.05 (Merge with Chapon Powder_multi)
* Modified by: PW, KL 05.06.07 (Concentricity)
* Modified by: EF, 17.10.08 (added any shape sample geometry)
* Modified by: EK, 28.10.10 (for X-ray use)
* Modified by: EK, (Added options for incoherent scattering)
*
* General powder sample (N lines, single scattering, incoherent scattering)
*
* %D
* General powder sample with
* many scattering vectors
* possibility for intrinsic line broadening
* incoherent background ratio is computed from material datafile.
* No multiple scattering. No secondary extinction.
*
* Based on Powder1/Powder2/Single_crystal.
* Geometry is a powder filled cylinder, sphere, box or any shape from an OFF file.
* Incoherent scattering is only provided here to account for a background.
* The efficient is highly improved when restricting the vertical scattering
* range on the Debye-Scherrer cone (with 'd_phi' and 'focus_flip').
* The unit cell volume Vc may also be computed when giving the density,
* the atomic/molecular weight and the number of atoms per unit cell.
*
* <b>Sample shape:</b>
* Sample shape may be a cylinder, a sphere, a box or any other shape.
* box/plate: xwidth x yheight x zdepth (thickness=0)
* hollow box/plate:xwidth x yheight x zdepth and thickness>0
* cylinder: radius x yheight (thickness=0)
* hollow cylinder: radius x yheight and thickness>0
* sphere: radius (yheight=0 thickness=0)
* hollow sphere: radius and thickness>0 (yheight=0)
* any shape: geometry=OFF_file
*
* The complex geometry option handles any closed non-convex polyhedra.
* It computes the intersection points of the xray with the object
* transparently, so that it can be used like a regular sample object.
* It supports the OFF and NOFF file format but not COFF (colored faces).
* Such files may be generated from XYZ data using:
* qhull < coordinates.xyz Qx Qv Tv o > geomview.off
* or
* powercrust coordinates.xyz
* and viewed with geomview or java -jar jroff.jar (see below).
* The default size of the object depends of the OFF file data, but its
* bounding box may be resized using xwidth,yheight and zdepth.
*
* If you use this component and produce valuable scientific results, please
* cite authors with references bellow (in <a href="#links">Links</a>).
*
* Example: PowderN(reflections = "c60.lau", d_phi = 15 , radius = 0.01,
* yheight = 0.05, Vc = 1076.89, Delta_d=0, DW=1,
* format=Crystallographica)
*
* <b>Powder definition file format</b>
* Powder structure is specified with an ascii data file 'reflections'.
* The powder data are free-text column based files.
* The reflection list should be ordered by decreasing d-spacing values.
* ... d ... F2
* Lines begining by '#' are read as comments (ignored) but they may contain
* the following keywords (in the header):
* #Vc <value of unit cell volume Vc [Angs^3]>
* #Debye_Waller <value of Debye-Waller factor DW>
* #Delta_d/d <value of Detla_d/d width for all lines>
* These values are not read if entered as component parameters (Vc=...)
*
* The signification of the columns in the numerical block may be
* set using the 'format' parameter. Built-in formats are:
* format=Crystallographica
* format=Fullprof
* format=Lazy
* and these specifications it is important NOT to use quotes, as shown.
*
* An other possibility to define other formats is to set directly
* the signification of the columns as a vector of indexes in the order
* format={j,d,F2,DW,Delta_d/d,1/2d,q,F}
* Signification of the symbols is given below. Indexes start at 1.
* Indices with zero means that the column are not present, so that:
* Crystallographica={ 4,5,7,0,0,0,0,0 }
* Fullprof ={ 4,0,8,0,0,5,0,0 }
* Lazy ={17,6,0,0,0,0,0,13}
* Here again, NO quotes should be around the 'format' value.
*
* At last, the format may be overridden by direct definition of the
* column indexes in the file itself by using the following keywords
* in the header (e.g. '#column_j 4'):
* #column_j <index of the multiplicity 'j' column>
* #column_d <index of the d-spacing 'd' column [Angs]>
* #column_F2 <index of the squared str. factor '|F|^2' column [b]>
* #column_F <index of the structure factor norm '|F|' column>
* #column_DW <index of the Debye-Waller factor 'DW' column>
* #column_Dd <index of the relative line width Delta_d/d 'Dd' column>
* #column_inv2d <index of the 1/2d=sin(theta)/lambda 'inv2d' column>
* #column_q <index of the scattering wavevector 'q' column [Angs-1]>
*
* <b>Concentricity</b>
*
* PowderN assumes 'concentric' shape, i.e. can contain other components inside its
* optional inner hollow. Example, Sample in Al cryostat:
*
*
* COMPONENT Cryo = PowderN(reflections="Al.laz", radius = 0.01, thickness = 0.001,
* concentric = 1, p_interact=0.1)
* AT (0,0,0) RELATIVE Somewhere
*
* COMPONENT Sample = some_other_component(with geometry FULLY enclosed in the hollow)
* AT (0,0,0) RELATIVE Somewhere
*
* COMPONENT Cryo2 = COPY(Cryo)(concentric = 0)
* AT (0,0,0) RELATIVE Somewhere
*
*
* (The second instance of the cryostat component can also be written out completely
* using PowderN(...). In both cases, this second instance needs concentric = 0.)
* The concentric arrangment can not be used with OFF geometry specification.
*
* %P
* INPUT PARAMETERS
* radius: Outer radius of sample in (x,z) plane [m]
* xwidth: Horiz. dimension of sample, as a width [m]
* yheight: Height of sample y direction [m]
* zdepth: Depth of box sample [m]
* thickness:Thickness of hollow sample.
* Negative value extends the hollow volume outside of the box/cylinder. [m]
* reflections: Input file for reflections.
* Use only incoherent scattering if NULL or "" [string]
*
* Optional parameters:
* d_phi: [deg] Angle corresponding to the vertical angular range to focus to, e.g. detector height. 0 for no focusing.
* d_theta: [deg] Longitudinal angle to focus to, i.e. lateral opening angle. This means, the maximum scattering angle, 2*theta, is equal to d_theta/2.0
* focus_flip: [1] Controls the sense of d_phi. If 0 d_phi is measured against the xz-plane. If !=0 d_phi is measured against zy-plane.
* pack: [1] Packing factor
* Delta_d: [AA] Global relative Delta_d/d spreading when the 'w' column is not available. Use 0 if ideal.
* format: [ ] Name of the format, or list of column indexes (see Description). N.b. no qoutes!
* p_inc: [1] Fraction of incoherently scattered rays.
* p_transmit: [1] Fraction of transmitted (only attenuated) rays.
* p_interact: [1] Fraction of events interacting with sample, e.g. 1-p_transmit-p_inc.
* concentric: [1] Indicate that this component has a hollow geometry and may contain other components. It should then be duplicated after the inside part (only for box, cylinder, sphere).
* geometry: [str] Name of an Object File Format (OFF) or PLY file for complex geometry. The OFF/PLY file may be generated from XYZ coordinates using qhull/powercrust.
* barns: [1] Flag to indicate if |F|^2 from 'reflections' is in barns or fm^2, (barns=1 for laz, barns=0 for lau type files).
* Vc: [AA^3] Volume of unit cell=nb atoms per cell/density of atoms.
* DW: [1] Global Debye-Waller factor when the 'DW' column is not available. Use 1 if included in F2.
* weight: [g/mol] Atomic/molecular weight of material [g/mol]
* density: [g/cm^3] Density of material. rho=density/weight/1e24*N_A.
* nb_atoms:[1] Number of sub-unit per unit cell, that is ratio of sigma for chemical formula to sigma per unit cell.
*
* OUTPUT PARAMETERS:
* line_info: [struct] Internal structure containing many members/info [struct]
* line_info.type: [char] Interaction type of event 't'=Transmit, 'i'=Incoherent, 'c'=Coherent [char]
* line_info.dq: [Angs^-1] Wavevector transfer of last coherent scattering event.
*
* %L
* See also: Single_crystal
* %L
* See <a href="http://icsd.ill.fr">ICSD</a> Inorganic Crystal Structure Database
* %L
* <a href="http://www.ncnr.nist.gov/resources/n-lengths/">Cross sections for single elements</a>
* %L
* <a href="http://www.ncnr.nist.gov/resources/sldcalc.html>Cross sections for compounds</a>
* %L
* <a href="http://www.webelements.com/">Web Elements</a>
* %L
* <a href="http://www.ill.eu/sites/fullprof/index.html">Fullprof</a> powder refinement
* %L
* <a href="http://www.crystallographica.com/">Crystallographica</a> software (free license)
* %L
* <a href="http://www.geomview.org">Geomview and Object File Format (OFF)</a>
* %L
* Java version of Geomview (display only) <a href="http://www.holmes3d.net/graphics/roffview/">jroff.jar</a>
* %L
* <a href="http://qhull.org">qhull</a>
* %L
* <a href="http://www.cs.ucdavis.edu/~amenta/powercrust.html">powercrust</a>
*
* %E
*****************************************************************************/
DEFINE COMPONENT PowderN
DEFINITION PARAMETERS (format=Undefined, mat_format=mat_Undefined)
SETTING PARAMETERS (string reflections=0, string material=0, string geometry=0,
radius=0, yheight=0, xwidth=0, zdepth=0, thickness=0,
pack=1, Vc=0, Delta_d=0, p_inc=0.1, p_transmit=0.1,
DW=0, nb_atoms=1, d_phi=0, d_theta=360, p_interact=0,
concentric=0, density=0, weight=0, barns=1, focus_flip=0)
OUTPUT PARAMETERS (line_info,columns,mat_columns,offdata)
/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */
SHARE
%{
/* used for reading data table from file */
%include "read_table-lib"
%include "interoff-lib"
/* Declare structures and functions only once in each instrument. */
#ifndef POWDERN_DECL
#define POWDERN_DECL
/* format definitions in the order {j d F2 DW Dd inv2d q F} */
#ifndef Crystallographica
#define Crystallographica { 4,5,7,0,0,0,0,0 }
#define Fullprof { 4,0,8,0,0,5,0,0 }
#define Lazy {17,6,0,0,0,0,0,13 }
#define Undefined { 0,0,0,0,0,0,0,0 }
/*format definitions for material data {E my_abs my_cohinc my_coh}*/
#define FFAST_element {1, 4, 0, 5, 0}
#define FFAST_compound {1, 3, 0, 0, 4}
#define XCOM {1, 4, 3, 0, 0}
#define mat_Undefined {0,0,0,0,0}
#endif
struct line_data
{
double F2; /* Value of structure factor */
double q; /* Qvector */
int j; /* Multiplicity */
double DWfactor; /* Debye-Waller factor */
double w; /* Intrinsic line width */
};
struct abs_data
{
double E; /*energy (in keV)*/
double k; /*wavenumber corresponding to E*/
double sigma_a; /*absorption cross section for energy E*/
double mu; /*absoprtion coefficient for the energy E*/
double cohinc; /*coherent +incoherent scattering coefficient*/
};
struct line_info_struct
{
struct line_data *list; /* Reflection array */
int count; /* Number of reflections */
double Dd;
double DWfactor;
double V_0;
double rho;
double at_weight;
double at_nb;
double sigma_a;
double sigma_i;
char compname[256];
double flag_barns;
int shape; /* 0 cylinder, 1 box, 2 sphere, 3 OFF file */
int column_order[8]; /* column signification */
int flag_warning;
char type; /* interaction type of event t=Transmit, i=Incoherent, c=Coherent */
double dq; /* wavevector transfer [Angs-1] */
double Epsilon; /* global strain in ppm */
double XsectionFactor;
double my_s_k2_sum;
double my_a;
double my_inc;
double *w,*q, *my_s_k2;
double radius_i,xwidth_i,yheight_i,zdepth_i;
double k; /* last velocity (cached) */
double Nq;
int nb_reuses, nb_refl, nb_refl_count;
double k_min, k_max;
double xs_Nq[CHAR_BUF_LENGTH];
double xs_sum[CHAR_BUF_LENGTH];
double photon_passed;
long xs_compute, xs_reuse, xs_calls;
t_Table mat_table;
int mat_column_order[5]; /*column signification for the coeff. in material data file*/
};
off_struct offdata;
// PN_list_compare *****************************************************************
int PN_list_compare (void const *a, void const *b)
{
struct line_data const *pa = a;
struct line_data const *pb = b;
double s = pa->q - pb->q;
if (!s) return 0;
else return (s < 0 ? -1 : 1);
} /* PN_list_compare */
int read_abs_data(char *ABS_file, struct abs_data **abs)
{
t_Table table;
char **parsing;
int status,i;
if (!ABS_file || !strlen(ABS_file) || !strcmp(ABS_file, "NULL")) {
fprintf(stderr,"Warning: material file (%s) not found.\n",ABS_file);
*abs=calloc(2,sizeof(struct abs_data));
(*abs)[1].E=-1;
(*abs)[1].k=-1;
(*abs)[1].mu=-1;
}else{
if ( (status=Table_Read(&table,ABS_file,0))==-1){
fprintf(stderr,"Error: %s Could not parse file \"%s\"\n",NAME_CURRENT_COMP,ABS_file);
exit(-1);
}
parsing=Table_ParseHeader(table.header,"Z","A[r]","rho",NULL);
*abs=calloc(table.rows+1,sizeof(struct abs_data));
for (i=0;i<table.rows;i++){
(*abs)[i].E=table.data[i*table.columns];
(*abs)[i].k=(*abs)[i].E*E2K;
// get from cox_76 abs_info.data[i].sigma_a=aba_info.data[i]*
(*abs)[i].mu=100*table.data[i*table.columns + 3];
/* column 4 (but 0-indexed) is photoabsorption in cm^2 / g
* density is in g / cm^3 ->factor 100 to get mu in m^-1.*/
(*abs)[i].cohinc=100*table.data[i*table.columns + 4];
/* column 5 contains the coh+inc mean free paths (mus). in cm^2 / g.*/
}
(*abs)[i].E=-1;
(*abs)[i].k=-1;
(*abs)[i].mu=-1;
Table_Free(&table);
return 1;
}
}
int read_line_data(char *SC_file, struct line_info_struct *info)
{
struct line_data *list = NULL;
int size = 0;
t_Table sTable; /* sample data table structure from SC_file */
int i=0;
int mult_count =0;
char flag=0;
double q_count=0, j_count=0, F2_count=0;
char **parsing;
int list_count=0;
if (!SC_file || !strlen(SC_file) || !strcmp(SC_file, "NULL")) {
MPI_MASTER(
printf("PowderN: %s: Using incoherent elastic scattering only.\n",
info->compname);
);
info->count = 0;
return(0);
}
Table_Read(&sTable, SC_file, 1); /* read 1st block data from SC_file into sTable*/
/* parsing of header */
parsing = Table_ParseHeader(sTable.header,
"Vc","V_0",
"column_j",
"column_d",
"column_F2",
"column_DW",
"column_Dd",
"column_inv2d", "column_1/2d", "column_sintheta/lambda",
"column_q", /* 10 */
"DW", "Debye_Waller",
"Delta_d/d",
"column_F ",
"V_rho",
"density",
"weight",
"nb_atoms","multiplicity",
NULL);
if (parsing) {
if (parsing[0] && !info->V_0) info->V_0 =atof(parsing[0]);
if (parsing[1] && !info->V_0) info->V_0 =atof(parsing[1]);
if (parsing[2]) info->column_order[0]=atoi(parsing[2]);
if (parsing[3]) info->column_order[1]=atoi(parsing[3]);
if (parsing[4]) info->column_order[2]=atoi(parsing[4]);
if (parsing[5]) info->column_order[3]=atoi(parsing[5]);
if (parsing[6]) info->column_order[4]=atoi(parsing[6]);
if (parsing[7]) info->column_order[5]=atoi(parsing[7]);
if (parsing[8]) info->column_order[5]=atoi(parsing[8]);
if (parsing[9]) info->column_order[5]=atoi(parsing[9]);
if (parsing[10]) info->column_order[6]=atoi(parsing[10]);
if (parsing[11] && info->DWfactor<=0) info->DWfactor=atof(parsing[11]);
if (parsing[12] && info->DWfactor<=0) info->DWfactor=atof(parsing[12]);
if (parsing[13] && info->Dd <0) info->Dd =atof(parsing[13]);
if (parsing[14]) info->column_order[7]=atoi(parsing[14]);
if (parsing[15] && !info->V_0) info->V_0 =1/atof(parsing[15]);
if (parsing[16] && !info->rho) info->rho =atof(parsing[16]);
if (parsing[17] && !info->at_weight) info->at_weight =atof(parsing[17]);
if (parsing[18] && info->at_nb <= 1) info->at_nb =atof(parsing[18]);
if (parsing[19] && info->at_nb <= 1) info->at_nb =atof(parsing[19]);
for (i=0; i<=19; i++) if (parsing[i]) free(parsing[i]);
free(parsing);
}
if (!sTable.rows)
exit(fprintf(stderr, "PowderN: %s: Error: The number of rows in %s "
"should be at least %d\n", info->compname, SC_file, 1));
else
size = sTable.rows;
MPI_MASTER(
Table_Info(sTable);
printf("PowderN: %s: Reading %d rows from %s\n",
info->compname, size, SC_file);
);
if (info->column_order[0] == 4 && info->flag_barns !=0)
MPI_MASTER(
printf("PowderN: %s: Powder file probably of type Crystallographica/Fullprof (lau)\n"
"WARNING: but F2 unit is set to barns=1 (barns). Intensity might be 100 times too high.\n",
info->compname);
);
if (info->column_order[0] == 17 && info->flag_barns == 0)
MPI_MASTER(
printf("PowderN: %s: Powder file probably of type Lazy Pulver (laz)\n"
"WARNING: but F2 unit is set to barns=0 (fm^2). Intensity might be 100 times too low.\n",
info->compname);
);
/* allocate line_data array */
list = (struct line_data*)malloc(size*sizeof(struct line_data));
for (i=0; i<size; i++)
{
/* printf("Reading in line %i\n",i);*/
double j=0, d=0, w=0, q=0, DWfactor=0, F2=0;
int index;
if (info->Dd >= 0) w = info->Dd;
if (info->DWfactor > 0) DWfactor = info->DWfactor;
/* get data from table using columns {j d F2 DW Dd inv2d q F} */
/* column indexes start at 1, thus need to substract 1 */
if (info->column_order[0] >0)
j = Table_Index(sTable, i, info->column_order[0]-1);
if (info->column_order[1] >0)
d = Table_Index(sTable, i, info->column_order[1]-1);
if (info->column_order[2] >0)
F2 = Table_Index(sTable, i, info->column_order[2]-1);
if (info->column_order[3] >0)
DWfactor = Table_Index(sTable, i, info->column_order[3]-1);
if (info->column_order[4] >0)
w = Table_Index(sTable, i, info->column_order[4]-1);
if (info->column_order[5] >0)
{ d = Table_Index(sTable, i, info->column_order[5]-1);
d = (d > 0? 1/d/2 : 0); }
if (info->column_order[6] >0)
{ q = Table_Index(sTable, i, info->column_order[6]-1);
d = (q > 0 ? 2*PI/q : 0); }
if (info->column_order[7] >0 && !F2)
{ F2 = Table_Index(sTable, i, info->column_order[7]-1); F2 *= F2; }
/* assign and check values */
j = (j > 0 ? j : 0);
q = (d > 0 ? 2*PI/d : 0); /* this is q */
DWfactor = (DWfactor > 0 ? DWfactor : 1);
w = (w>0 ? w : 0); /* this is q and d relative spreading */
F2 = (F2 >= 0 ? F2 : 0);
if (j == 0 || q == 0) {
MPI_MASTER(
printf("PowderN: %s: line %i has invalid definition\n"
" (mult=0 or q=0 or d=0)\n", info->compname, i);
);
continue;
}
list[list_count].j = j;
list[list_count].q = q;
list[list_count].DWfactor = DWfactor;
list[list_count].w = w;
list[list_count].F2= F2;
/* adjust multiplicity if j-column + multiple d-spacing lines */
/* if d = previous d, increase line duplication index */
if (!q_count) q_count = q;
if (!j_count) j_count = j;
if (!F2_count) F2_count = F2;
if (fabs(q_count-q) < 0.0001*fabs(q)
&& fabs(F2_count-F2) < 0.0001*fabs(F2) && j_count == j) {
mult_count++; flag=0; }
else flag=1;
if (i == size-1) flag=1;
/* else if d != previous d : just passed equivalent lines */
if (flag) {
if (i == size-1) list_count++;
/* if duplication index == previous multiplicity */
/* set back multiplicity of previous lines to 1 */
if ((mult_count && list_count>0)
&& (mult_count == list[list_count-1].j
|| ((list_count < size) && (i == size - 1)
&& (mult_count == list[list_count].j))) ) {
MPI_MASTER(
printf("PowderN: %s: Set multiplicity to 1 for lines [%i:%i]\n"
" (d-spacing %g is duplicated %i times)\n",
info->compname, list_count-mult_count, list_count-1, list[list_count-1].q, mult_count);
);
for (index=list_count-mult_count; index<list_count; list[index++].j = 1);
mult_count = 1;
q_count = q;
j_count = j;
F2_count = F2;
}
if (i == size-1) list_count--;
flag=0;
}
list_count++;
} /* end for */
Table_Free(&sTable);
/* sort the list with increasing q */
qsort(list, list_count, sizeof(struct line_data), PN_list_compare);
MPI_MASTER(
printf("PowderN: %s: Read %i reflections from file '%s'\n",
info->compname, list_count, SC_file);
);
info->list = list;
info->count = list_count;
return(list_count);
} /* read_line_data */
/* computes the number of possible reflections (return value), and the total xsection 'sum' */
/* this routine looks for a pre-computed value in the Nq and sum cache tables */
/* when found, the earch starts from the corresponding lower element in the table */
int calc_xsect(double k, double *q, double *my_s_k2, int count, double tth_max, double *sum,
struct line_info_struct *line_info) {
int Nq = 0, line=0, line0=0;
double sinth=sin(DEG2RAD*tth_max*0.5);
*sum=0;
/* check if a line_info element has been recorded already */
if (k >= line_info->k_min && k <= line_info->k_max && line_info->photon_passed >= CHAR_BUF_LENGTH) {
line = (int)floor(k - line_info->k_min)*CHAR_BUF_LENGTH/(line_info->k_max - line_info->k_min);
Nq = line_info->xs_Nq[line];
*sum = line_info->xs_sum[line];
if (!Nq && *sum == 0) {
/* not yet set: we compute the sum up to the corresponding wavevector in the table cache */
double line_k = line_info->k_min + line*(line_info->k_max - line_info->k_min)/CHAR_BUF_LENGTH;
for(line0=0; line0<count; line0++) {
if (q[line0] <= 2*line_k*sinth) { /* q < 2*kf: restrict structural range */
*sum += my_s_k2[line0];
if (Nq < line0+1) Nq=line0+1; /* determine maximum line index which can scatter */
} else break;
}
line_info->xs_Nq[line] = Nq;
line_info->xs_sum[line]= *sum;
line_info->xs_compute++;
} else line_info->xs_reuse++;
line0 = Nq - 1;
}
line_info->xs_calls++;
for(line=line0; line<count; line++) {
if (q[line] <= 2*k*sinth) { /* q < 2*kf: restrict structural range */
*sum += my_s_k2[line];
if (Nq < line+1) Nq=line+1; /* determine maximum line index which can scatter */
} else break;
}
return(Nq);
} /* calc_xsect */
int calc_abs_xsect(double k, double *abs, struct line_info_struct *line_info){
/* compute the absorption cross section.
* amounts to a table lookup in material data_file by energy*/
double e=k*K2E;
*abs=Table_Value((line_info->mat_table), e, line_info->mat_column_order[1]);
*abs*=100*line_info->rho;
return 0;
}
int calc_inc_xsect(double k, double *inc, struct line_info_struct *line_info){
/*TODO choose one of these bits*/
/* compute the incoherent cross section.
* amounts to a table lookup in material data_file by energy
* In time this should be improved to actually handle Compton scattering.*/
double e=k*K2E;
double correction;
int i=2;
/* To find the inc. scattering column in material datafiles,
* Find the first non-zero column index (btw. 2..4).*/
for (i=2;i<5;i++){
if (line_info->mat_column_order[i]) break;
}
switch (i){
case 2:
/*inc is explicitly in column i - no correction term.*/
correction=0;
break;
case 3:
/*inc is given as inc+ćoh in column i*/
correction=line_info->my_s_k2_sum/(k*k);
break;
case 4:
/*inc is given as total attenuation = inc+ćoh+abs in column i*/
correction=line_info->my_s_k2_sum/(k*k) + line_info->my_a;
break;
case 5:
/*means no non-zero index has been found - set xsect to 0*/
*inc=0;
return 1;
}
*inc=Table_Value((line_info->mat_table), e, line_info->mat_column_order[i]-1)-correction;
if (*inc<0){
*inc=0;
}
*inc*=100*line_info->rho;/*scale by rho and 100 to get in units m^-1*/
return 0;
}
#endif /* !POWDERN_DECL */
%}
DECLARE
%{
struct line_info_struct line_info;
int *columns;
int *mat_columns;
off_struct offdata;
double type;
%}
INITIALIZE
%{
columns = (int[])format;
mat_columns = (int[]) mat_format;
int i=0;
struct line_data *L;
line_info.Dd = Delta_d;
line_info.DWfactor = DW;
line_info.V_0 = Vc;
line_info.rho = density;
line_info.at_weight= weight;
line_info.at_nb = nb_atoms;
line_info.flag_barns=barns;
line_info.shape = 0;
line_info.flag_warning=0;
line_info.radius_i =line_info.xwidth_i=line_info.yheight_i=line_info.zdepth_i=0;
line_info.k = 0;
line_info.Nq = 0;
line_info.k_min = FLT_MAX; line_info.k_max = 0;
line_info.photon_passed=0;
line_info.nb_reuses = line_info.nb_refl = line_info.nb_refl_count = 0;
line_info.xs_compute= line_info.xs_reuse= line_info.xs_calls =0;
for (i=0; i< 8; i++) line_info.column_order[i] = columns[i];
for (i=0; i< 4; i++) line_info.mat_column_order[i] = mat_columns[i];
strncpy(line_info.compname, NAME_CURRENT_COMP, 256);
line_info.shape=-1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape */
if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
if (off_init(geometry, xwidth, yheight, zdepth, 0, &offdata)) {
line_info.shape=3; thickness=0; concentric=0;
}
}
else if (xwidth && yheight && zdepth) line_info.shape=1; /* box */
else if (radius > 0 && yheight) line_info.shape=0; /* cylinder */
else if (radius > 0 && !yheight) line_info.shape=2; /* sphere */
if (line_info.shape < 0)
exit(fprintf(stderr,"PowderN: %s: sample has invalid dimensions.\n"
"ERROR Please check parameter values (xwidth, yheight, zdepth, radius).\n", NAME_CURRENT_COMP));
if (thickness) {
if (radius && (radius < fabs(thickness))) {
MPI_MASTER(
printf("PowderN: %s: hollow sample thickness is larger than its volume (sphere/cylinder).\n"
"WARNING Please check parameter values. Using bulk sample (thickness=0).\n", NAME_CURRENT_COMP);
);
thickness=0;
}
else if (!radius && (xwidth < 2*fabs(thickness) || yheight < 2*fabs(thickness) || zdepth < 2*fabs(thickness))) {
MPI_MASTER(
printf("PowderN: %s: hollow sample thickness is larger than its volume (box).\n"
"WARNING Please check parameter values.\n", NAME_CURRENT_COMP);
);
}
}
if (concentric && thickness==0) {
MPI_MASTER(
printf("PowderN: %s:Can not use concentric mode\n"
"WARNING on non hollow shape. Ignoring.\n",
NAME_CURRENT_COMP);
);
concentric=0;
}
if (thickness>0) {
if (radius>thickness) {
line_info.radius_i=radius-thickness;
} else {
if (xwidth>2*thickness) line_info.xwidth_i =xwidth -2*thickness;
if (yheight>2*thickness) line_info.yheight_i=yheight-2*thickness;
if (zdepth>2*thickness) line_info.zdepth_i =zdepth -2*thickness;
}
} else if (thickness<0) {
thickness = fabs(thickness);
if (radius) {
line_info.radius_i=radius;
radius=line_info.radius_i+thickness;
} else {
line_info.xwidth_i =xwidth;
line_info.yheight_i=yheight;
line_info.zdepth_i =zdepth;
xwidth =xwidth +2*thickness;
yheight =yheight+2*thickness;
zdepth =zdepth +2*thickness;
}
}
if (!line_info.yheight_i) {
line_info.yheight_i = yheight;
}
if (p_interact) {
if (p_interact < p_inc) { double tmp=p_interact; p_interact=p_inc; p_inc=tmp; }
p_transmit = 1-p_interact-p_inc;
}
if (p_inc + p_transmit > 1) {
MPI_MASTER(
printf("PowderN: %s: You have requested an unmeaningful choice of the 'p_inc' and 'p_transmit' parameters (sum is %g, exeeding 1). Fixing.\n",
NAME_CURRENT_COMP, p_inc+p_transmit);
);
if (p_inc > p_transmit) p_transmit=1-2*p_inc;
else p_transmit=1-2*p_inc;
} else if (p_inc + p_transmit == 1) {
MPI_MASTER(
printf("PowderN: %s: You have requested all photons be attenuated\n"
"WARNING or incoherently scattered!\n", NAME_CURRENT_COMP);
);
}
if (concentric) {
MPI_MASTER(
printf("PowderN: %s: Concentric mode - remember to include the 'opposite' copy of this component !\n"
"WARNING The equivalent, 'opposite' comp should have concentric=0\n", NAME_CURRENT_COMP);
);
if (p_transmit == 0) {
MPI_MASTER(
printf("PowderN: %s: Concentric mode and p_transmit==0 !?\n"
"WARNING Don't you want any transmitted xrays?\n", NAME_CURRENT_COMP);
);
}
}
if (reflections && strlen(reflections) && strcmp(reflections, "NULL") && strcmp(reflections, "0")) {
i = read_line_data(reflections, &line_info);
if (i == 0)
exit(fprintf(stderr,"PowderN: %s: reflection file %s is not valid.\n"
"ERROR Please check file format (laz or lau).\n", NAME_CURRENT_COMP, reflections));
}
/* compute the scattering unit density from material weight and density */
/* the weight of the scattering element is the chemical formula molecular weight
* times the nb of chemical formulae in the scattering element (nb_atoms) */
if (!line_info.V_0 && line_info.at_nb > 0
&& line_info.at_weight > 0 && line_info.rho > 0) {
/* molar volume [cm^3/mol] = weight [g/mol] / density [g/cm^3] */
/* atom density per Angs^3 = [mol/cm^3] * N_Avogadro *(1e-8)^3 */
line_info.V_0 = line_info.at_nb
/(line_info.rho/line_info.at_weight/1e24*6.02214199e23);
}
/* the scattering unit cross sections are the chemical formula onces
* times the nb of chemical formulae in the scattering element */
if (line_info.at_nb > 0) {
line_info.sigma_a *= line_info.at_nb; line_info.sigma_i *= line_info.at_nb;
}
if (line_info.sigma_a<0) line_info.sigma_a=0;
if (line_info.sigma_i<0) line_info.sigma_i=0;
if (line_info.V_0 <= 0)
MPI_MASTER(
printf("PowderN: %s: density/unit cell volume is NULL (Vc). Unactivating component.\n", NAME_CURRENT_COMP);
);
if (line_info.V_0 > 0 && p_inc && (!line_info.sigma_i || !material)) {
MPI_MASTER(
printf("PowderN: %s: WARNING: You have requested statistics for incoherent scattering but not set a material datafile (\'%s\') nor defined sigma_inc!\n", NAME_CURRENT_COMP, material);
);
}
if (line_info.flag_barns) { /* Factor 100 to convert from barns to fm^2 */
line_info.XsectionFactor = 100;
} else {
line_info.XsectionFactor = 1;
}
if (line_info.V_0 > 0 && i) {
L = line_info.list;
line_info.q = malloc(line_info.count*sizeof(double));
line_info.w = malloc(line_info.count*sizeof(double));
line_info.my_s_k2 = malloc(line_info.count*sizeof(double));
if (!line_info.q || !line_info.w || !line_info.my_s_k2)
exit(fprintf(stderr,"PowderN: %s: ERROR allocating memory (init)\n", NAME_CURRENT_COMP));
for(i=0; i<line_info.count; i++)
{
line_info.my_s_k2[i] = 4*PI*PI*PI*pack*(L[i].DWfactor ? L[i].DWfactor : 1)
/(line_info.V_0*line_info.V_0)
*(L[i].j * L[i].F2 / L[i].q)*line_info.XsectionFactor;
/* Is not yet divided by k^2 */
/* Squires [3.103] */
line_info.q[i] = L[i].q;
line_info.w[i] = L[i].w;
}
}
if (material || !strlen(material) || !strcmp(material, "NULL")) {
int status;
char **parsing;
if( (status=Table_Read(&(line_info.mat_table),material,0))==-1){
fprintf(stderr,"PowderN: %s Error reading material data from file %s.\n",NAME_CURRENT_COMP,material);
}
parsing=Table_ParseHeader(line_info.mat_table.header,
"column_e","column_abs","column_inc","column_cohinc","column_tot",NULL);
if (parsing){
int i;
for (i=0;i<5;i++){
if (parsing[i]) line_info.mat_column_order[i]=atoi(parsing[i]);
}
}
}
%}
TRACE
%{
double l0, l1, l2, l3, k, k1,l_full, l, l_1, dl, alpha0, alpha, theta, my_s, my_s_n;
double solid_angle;//, type;
double arg, tmp_kx, tmp_ky, tmp_kz, kout_x, kout_y, kout_z, nx, ny, nz, pmul=1;
int line;
char intersect=0;
char intersecti=0;
line_info.type = '\0';
if (line_info.V_0 > 0 && (line_info.count || line_info.my_inc)) {
if (line_info.shape == 1) {
intersect = box_intersect(&l0, &l3, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);
intersecti = box_intersect(&l1, &l2, x, y, z, kx, ky, kz, line_info.xwidth_i, line_info.yheight_i, line_info.zdepth_i);
} else if (line_info.shape == 0) {
intersect = cylinder_intersect(&l0, &l3, x, y, z, kx, ky, kz, radius, yheight);
intersecti = cylinder_intersect(&l1, &l2, x, y, z, kx, ky, kz, line_info.radius_i, line_info.yheight_i);
} else if (line_info.shape == 2) {
intersect = sphere_intersect (&l0, &l3, x,y,z, kx,ky,kz, radius);
intersecti = sphere_intersect (&l1, &l2, x,y,z, kx,ky,kz, line_info.radius_i);
} else if (line_info.shape == 3) {
intersect = off_x_intersect (&l0, &l3, NULL, NULL, x,y,z, kx,ky,kz, offdata);
intersecti = 0;
}
}
if(intersect && l3 >0) {
if (concentric) {
/* Set up for concentric case */
/* 'Remove' the backside of this comp */
if (!intersecti) {
l1 = (l3 + l0) /2;
}
l2 = l1;
l3 = l1;
dl = -1.0*rand01(); /* In case of scattering we will scatter on 'forward' part of sample */
} else {
if (!intersecti) {
l1 = (l3 + l0) /2;
l2 = l1;
}
dl = randpm1(); /* Possibility to scatter at all points in line of sight */
}
/* X-ray enters at t=l0. */
if(l0 < 0) l0=0; /* already in sample */
if(l1 < 0) l1=0; /* already in inner hollow */
if(l2 < 0) l2=0; /* already past inner hollow */
k = sqrt(kx*kx + ky*ky + kz*kz);
l_full =l3 - l2 + l1 - l0;
if (line_info.photon_passed < CHAR_BUF_LENGTH) {
if (k < line_info.k_min) line_info.k_min = k;
if (k > line_info.k_max) line_info.k_max = k;
line_info.photon_passed++;
}
/* Calculate total scattering cross section at relevant wavevector */
if ( fabs(k - line_info.k) < 1e-15) {
line_info.nb_reuses++;
} else {
line_info.Nq = calc_xsect(k, line_info.q, line_info.my_s_k2, line_info.count, d_theta*0.5, &line_info.my_s_k2_sum, &line_info);
calc_abs_xsect(k, &line_info.my_a, &line_info);
calc_inc_xsect(k, &line_info.my_inc, &line_info);
line_info.k = k;
line_info.nb_refl += line_info.Nq;
line_info.nb_refl_count++;
}
if (l3 < 0) {
l3=0; /* Already past sample?! */
if (line_info.flag_warning < 100)
printf("PowderN: %s: Warning: xray has already passed us? (Skipped).\n"
" In concentric geometry, this may be caused by a missing concentric=0 option in 2nd enclosing instance.\n", NAME_CURRENT_COMP);
line_info.flag_warning++;
} else {
if (dl<0) { /* Calculate scattering point position */
dl = fabs(dl)*(l1 - l0); /* 'Forward' part */
} else {
dl = dl * (l3 - l2) + (l2-l0) ; /* Possibly also 'backside' part */
}
my_s=line_info.my_s_k2_sum/(k*k)+line_info.my_inc;
/* Total attenuation from scattering*/
type = rand01();
/* How to handle this one? Transmit (1) / Incoherent (2) / Coherent (3) ? */
if (type < p_transmit) {
type = 1;
l = l_full; /* Passing through, full length */
PROP_DL(l3);
} else if (type >= p_transmit && type < (p_transmit + p_inc)) {
type = 2;
l = dl; /* Penetration in sample */
PROP_DL(dl+l0); /* Point of scattering */
SCATTER;
} else if (type >= p_transmit + p_inc) {
type = 3;
l=dl;
PROP_DL(dl+l0); /* Point of scattering */
SCATTER;
} else {
exit(fprintf(stderr,"PowderN %s: DEAD - this shouldn't happen!\n", NAME_CURRENT_COMP));
}
if (type == 3) { /* Make coherent scattering event */
if (line_info.count > 0) {
/* choose line */
if (line_info.Nq > 1) line=floor(line_info.Nq*rand01()); /* Select between Nq powder lines */
else line = 0;
if (line_info.w[line])
arg = line_info.q[line]*(1+line_info.w[line]*randnorm())/(2.0*k);
else
arg = line_info.q[line]/(2.0*k);
my_s_n = line_info.my_s_k2[line]/(k*k);
if(fabs(arg) > 1)
ABSORB; /* No bragg scattering possible*/
theta = asin(arg); /* Bragg scattering law */
/* Choose point on Debye-Scherrer cone */
if (d_phi)
{ /* relate height of detector to the height on DS cone */
arg = sin(d_phi*DEG2RAD/2)/sin(2*theta);
/* If full Debye-Scherrer cone is within d_phi, don't focus */
if (arg < -1 || arg > 1) d_phi = 0;
/* Otherwise, determine alpha to rotate from scattering plane
into d_phi focusing area*/
else alpha = 2*asin(arg);
}
if (d_phi) {
/* Focusing */
alpha = fabs(alpha);
/* Trick to get scattering for pos/neg theta's */
alpha0= 2*rand01()*alpha;
if (alpha0 > alpha) {
alpha0=M_PI+(alpha0-1.5*alpha);
} else {
alpha0=alpha0-0.5*alpha;
}
if(focus_flip){
alpha0+=M_PI_2;
}
}
else
alpha0 = M_PI*randpm1();
/* now find a nearly vertical rotation axis:
* Either
* (k along Z) x (X axis) -> nearly Y axis
* Or
* (k along X) x (Z axis) -> nearly Y axis
*/
if (fabs(scalar_prod(1,0,0,kx/k,ky/k,kz/k)) < fabs(scalar_prod(0,0,1,kx/k,ky/k,kz/k))) {
nx = 1; ny = 0; nz = 0;
} else {
nx = 0; ny = 0; nz = 1;
}
vec_prod(tmp_kx,tmp_ky,tmp_kz, kx,ky,kz, nx,ny,nz);
/* k_out = rotate 'k' by 2*theta around tmp_k: Bragg angle */
rotate(kout_x,kout_y,kout_z, kx,ky,kz, 2*theta, tmp_kx,tmp_ky,tmp_kz);
/* tmp_k = rotate k_out by alpha0 around 'k' (Debye-Scherrer cone) */
rotate(tmp_kx,tmp_ky,tmp_kz, kout_x,kout_y,kout_z, alpha0, kx, ky, kz);
kx = tmp_kx;
ky = tmp_ky;
kz = tmp_kz;
/*weight the outgoing signal according to polarization*/
if (Ex!=0 || Ey!=0 || Ez!=0){
double EE=sqrt(Ex*Ex+Ey*Ey+Ez*Ez);
double s=scalar_prod(kx,ky,kz,Ex,Ey,Ez)/k/EE;
p*=(1-s)*(1-s);
}else{
/*unpolarized light in - means an effective reduction according to only theta*/
p*=(1+cos(theta)*cos(theta))*0.5;
}
/* Since now scattered and new direction given, calculate path to exit */
if (line_info.shape == 1) {
intersect = box_intersect(&l0, &l3, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);
intersecti = box_intersect(&l1, &l2, x, y, z, kx, ky, kz, line_info.xwidth_i, line_info.yheight_i, line_info.zdepth_i);
} else if (line_info.shape == 0) {
intersect = cylinder_intersect(&l0, &l3, x, y, z, kx, ky, kz, radius, yheight);
intersecti = cylinder_intersect(&l1, &l2, x, y, z, kx, ky, kz, line_info.radius_i, line_info.yheight_i);
} else if (line_info.shape == 2) {
intersect = sphere_intersect (&l0, &l3, x,y,z, kx,ky,kz, radius);
intersecti = sphere_intersect (&l1, &l2, x,y,z, kx,ky,kz, line_info.radius_i);
} else if (line_info.shape == 3) {
intersect = off_x_intersect (&l0, &l3, NULL, NULL, x,y,z, kx,ky,kz, offdata);
intersecti = 0;
}
if (!intersect) {
/* Strange error: did not hit cylinder */
if (line_info.flag_warning < 100)
printf("PowderN: %s: WARNING: Did not hit sample from inside (coh). ABSORB.\n", NAME_CURRENT_COMP);
line_info.flag_warning++;
ABSORB;
}
if (!intersecti) {
l1 = (l3 + l0) /2;
l2 = l1;
}
if (concentric && intersecti) {
/* In case of concentricity, 'remove' backward wall of sample */
l2 = l1;
l3 = l1;
}
if(l0 < 0) l0=0; /* already in sample */
if(l1 < 0) l1=0; /* already in inner hollow */
if(l2 < 0) l2=0; /* already past inner hollow */
l_1 = l3 - l2 + l1 - l0; /* Length to exit */
pmul = line_info.Nq*l_full*my_s_n*exp(-(line_info.my_a+my_s)*(l+l_1))/(1-(p_inc+p_transmit));
/* Correction in case of d_phi focusing - BUT only when d_phi != 0 */
if (d_phi) pmul *= alpha/PI;
line_info.type = 'c';
line_info.dq = line_info.q[line];
} /* else transmit <-- No powder lines in file */
} /* Coherent scattering event */
else if (type == 2) { /* Make incoherent scattering event */
/*should be replaced by Compton scattering and/or TDS*/
if(d_phi) {
randvec_target_rect_angular(&kx, &ky, &kz, &solid_angle,
0, 0, 1,
2*PI, d_phi*DEG2RAD, ROT_A_CURRENT_COMP);
} else {
randvec_target_circle(&kx, &ky, &kz,
&solid_angle, 0, 0, 1, 0);
}
k1 = sqrt(kx*kx+ky*ky+kz*kz);
kx *= k/k1;
ky *= k/k1;
kz *= k/k1;
/* Since now scattered and new direction given, calculate path to exit */
if (line_info.shape == 1) {
intersect = box_intersect(&l0, &l3, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);
intersecti = box_intersect(&l1, &l2, x, y, z, kx, ky, kz, line_info.xwidth_i, line_info.yheight_i, line_info.zdepth_i);
} else if (line_info.shape == 0) {
intersect = cylinder_intersect(&l0, &l3, x, y, z, kx, ky, kz, radius, yheight);
intersecti = cylinder_intersect(&l1, &l2, x, y, z, kx, ky, kz, line_info.radius_i, line_info.yheight_i);
} else if (line_info.shape == 2) {
intersect = sphere_intersect (&l0, &l3, x,y,z, kx,ky,kz, radius);
intersecti = sphere_intersect (&l1, &l2, x,y,z, kx,ky,kz, line_info.radius_i);
} else if (line_info.shape == 3) {
intersect = off_x_intersect (&l0, &l3, NULL, NULL, x,y,z, kx,ky,kz, offdata);
intersecti = 0;
}
if (!intersect) {
/* Strange error: did not hit cylinder */
if (line_info.flag_warning < 100)
printf("PowderN: %s: WARNING: Did not hit sample from inside (inc). ABSORB.\n", NAME_CURRENT_COMP);
line_info.flag_warning++;
ABSORB;
}
if (!intersecti) {
l1 = (l3 + l0) /2;
l2 = l1;
}
if (concentric && intersecti) {
/* In case of concentricity, 'remove' backward wall of sample */
l2 = l1;
l3 = l1;
}
if(l0 < 0) l0=0; /* already in sample */
if(l1 < 0) l1=0; /* already in inner hollow */
if(l2 < 0) l2=0; /* already past inner hollow */
l_1 = (l3 - l2 + l1 - l0); /* Length to exit */
pmul = l_full*line_info.my_inc*exp(-(line_info.my_a+my_s)*(l+l_1))/(p_inc);
pmul *= solid_angle/(4*PI);
line_info.type = 'i';
} /* Incoherent scattering event */
else if (type == 1) {
/* Make transmitted (absorption-corrected) event */
/* No coordinate changes here, simply change xray weight */
pmul = exp(-(line_info.my_a+my_s)*(l))/(p_transmit);
line_info.type = 't';
}
p *= pmul;
} /* Photon leaving since it has passed already */
} /* else transmit non interacting xrays */
%}
FINALLY
%{
free(line_info.list);
free(line_info.q);
free(line_info.w);
free(line_info.my_s_k2);
MPI_MASTER(
if (line_info.flag_warning)
printf("PowderN: %s: Error messages were repeated %i times with absorbed xrays.\n",
NAME_CURRENT_COMP, line_info.flag_warning);
/* in case this instance is used in a SPLIT, we can recommend the
optimal iteration value */
if (line_info.nb_refl_count) {
double split_iterations = (double)line_info.nb_reuses/line_info.nb_refl_count + 1;
double split_optimal = (double)line_info.nb_refl/line_info.nb_refl_count;
if (split_optimal > split_iterations + 5)
printf("PowderN: %s: Info: you may highly improve the computation efficiency by using\n"
" SPLIT %i COMPONENT %s=PowderN(...)\n"
" in the instrument description %s.\n",
NAME_CURRENT_COMP, (int)split_optimal, NAME_CURRENT_COMP, mcinstrument_source);
}
);
%}
MCDISPLAY
%{
if (line_info.V_0) {
if (line_info.shape == 0) { /* cyl */
circle("xz", 0, yheight/2.0, 0, radius);
circle("xz", 0, -yheight/2.0, 0, radius);
line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0);
line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0);
line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius);
line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius);
if (thickness) {
double radius_i=radius-thickness;
circle("xz", 0, yheight/2.0, 0, radius_i);
circle("xz", 0, -yheight/2.0, 0, radius_i);
line(-radius_i, -yheight/2.0, 0, -radius_i, +yheight/2.0, 0);
line(+radius_i, -yheight/2.0, 0, +radius_i, +yheight/2.0, 0);
line(0, -yheight/2.0, -radius_i, 0, +yheight/2.0, -radius_i);
line(0, -yheight/2.0, +radius_i, 0, +yheight/2.0, +radius_i);
}
} else if (line_info.shape == 1) { /* box */
double xmin = -0.5*xwidth;
double xmax = 0.5*xwidth;
double ymin = -0.5*yheight;
double ymax = 0.5*yheight;
double zmin = -0.5*zdepth;
double zmax = 0.5*zdepth;
multiline(5, xmin, ymin, zmin,
xmax, ymin, zmin,
xmax, ymax, zmin,
xmin, ymax, zmin,
xmin, ymin, zmin);
multiline(5, xmin, ymin, zmax,
xmax, ymin, zmax,
xmax, ymax, zmax,
xmin, ymax, zmax,
xmin, ymin, zmax);
line(xmin, ymin, zmin, xmin, ymin, zmax);
line(xmax, ymin, zmin, xmax, ymin, zmax);
line(xmin, ymax, zmin, xmin, ymax, zmax);
line(xmax, ymax, zmin, xmax, ymax, zmax);
if (line_info.zdepth_i) {
xmin = -0.5*line_info.xwidth_i;
xmax = 0.5*line_info.xwidth_i;
ymin = -0.5*line_info.yheight_i;
ymax = 0.5*line_info.yheight_i;
zmin = -0.5*line_info.zdepth_i;
zmax = 0.5*line_info.zdepth_i;
multiline(5, xmin, ymin, zmin,
xmax, ymin, zmin,
xmax, ymax, zmin,
xmin, ymax, zmin,
xmin, ymin, zmin);
multiline(5, xmin, ymin, zmax,
xmax, ymin, zmax,
xmax, ymax, zmax,
xmin, ymax, zmax,
xmin, ymin, zmax);
line(xmin, ymin, zmin, xmin, ymin, zmax);
line(xmax, ymin, zmin, xmax, ymin, zmax);
line(xmin, ymax, zmin, xmin, ymax, zmax);
line(xmax, ymax, zmin, xmax, ymax, zmax);
}
} if (line_info.shape == 2) { /* sphere */
if (line_info.radius_i) {
circle("xy",0,0,0,line_info.radius_i);
circle("xz",0,0,0,line_info.radius_i);
circle("yz",0,0,0,line_info.radius_i);
}
circle("xy",0,0,0,radius);
circle("xz",0,0,0,radius);
circle("yz",0,0,0,radius);
} else if (line_info.shape == 3) { /* OFF file */
off_display(offdata);
}
}
%}
END