/*****************************************************************************
* McStas, neutron ray-tracing package
* Copyright 1997-2000 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: 1.24 $
* Origin: McStas release
* 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)
*
* General powder sample (N lines, single scattering, incoherent scattering)
*
* %D
* General powder sample with
* many scattering vectors
* possibility for intrinsic line broadening
* incoherent backgorund ratio is specified by user.
* No multiple scattering. No secondary extinction.
*
* Based on Powder1/Powder2/Single_crystal.
* Geometry is a powder filled cylinder or a box.
* 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').
* 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.
*
* Example: PowderN(reflections = "c60.lau", d_phi = 15 , radius = 0.01,
* yheight = 0.05, Vc = 1076.89, sigma_abs = 0, Delta_d=0, DW=1,
* format=Crystallographica)
*
* <b>Powder definition and file format</b>
*
* Powder structure is specified with an ascii data file 'reflections'.
* The powder data are free-text column based files.
* 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>
* #sigma_abs <value of Absorption cross section>
* #sigma_inc <value of Incoherent cross section>
* #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.
* Indexes with zero means that the column is 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>
* #column_F2 <index of the squared structure factor '|F|^2' column>
* #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>
*
* %P
* INPUT PARAMETERS
*
* d_phi: Angle corresponding to the vertical angular range
* to focus to, e.g. detector height. 0 for no focusing [deg,0-180]
* radius: Radius of sample in (x,z) plane [m]
* yheight: Height of sample y direction [m]
* Vc: Volume of unit cell=nb atoms per cell/density of atoms [AA^3]
* sigma_abs:Absorption cross section per unit cell at 2200 m/s [fm^2]
* sigma_inc:Incoherent cross section per unit cell [fm^2]
* reflections: Input file for reflections.
* Use only incoherent scattering if NULL or "" [string]
* format: name of the format, or list of column indexes
(see Description). [no quotes]
*
* Optional parameters:
* xwidth: horiz. dimension of sample, as a width [m]
* zthick: thickness of sample [m]
* h: the same as yheight [m]
* Delta_d: global relative Delta_d/d spreading when the 'w' column
* is not available. Use 0 if ideal. [Angs]
* DW: global Debey-Waller factor when the 'DW' column
* is not available. Use 1 if included in F2 [1]
* frac: Fraction of incoherently scattered neutron rays [1]
* pack: Packing factor [1]
* weight: atomic/molecular weight of material [g/mol]
* density: density of material. V_rho=density/weight/1e24*N_A. [g/cm^3]
* nb_atoms: number of atoms per unit cell [1]
*
* %L
* See also: Powder1, Powder2 and PowderN
* %L
* See <a href="http://icsd.ill.fr">ICSD</a> Inorganic Crystal Structure Database
* %L
* Cross sections for single elements: http://www.ncnr.nist.gov/resources/n-lengths/
* %L
* Cross sections for compounds: http://www.ncnr.nist.gov/resources/sldcalc.html
* %L
* Fullprof powder refinement: http://www-llb.cea.fr/fullweb/fp2k/fp2k.htm
* %L
* Web Elements http://www.webelements.com/
*
* %E
*****************************************************************************/
DEFINE COMPONENT PowderN
DEFINITION PARAMETERS (reflections, format=Undefined)
SETTING PARAMETERS (d_phi=0, radius=0.01, yheight=0.05,
pack=1, Vc=0, sigma_abs=0, sigma_inc=0, Delta_d=0, frac=0,
xwidth=0, zthick=0, h=0, DW=0, nb_atoms=1, density=0, weight=0)
OUTPUT PARAMETERS (line_info, Nq, my_s_v2,
my_s_v2_sum, my_a_v, my_inc, q_v, w_v, isrect, columns)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
SHARE
%{
/* used for reading data table from file */
%include "read_table-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,0 }
#define Undefined { 0,0,0,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 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];
int column_order[8]; /* column signification */
};
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)) {
printf("PowderN: %s: Using incoherent elastic scattering only\n",
info->compname);
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",
"sigma_abs","sigma_a ",
"sigma_inc","sigma_i ",
"column_j",
"column_d",
"column_F2",
"column_DW",
"column_Dd",
"column_inv2d", "column_1/2d", "column_sintheta/lambda",
"column_q", /* 14 */
"DW", "Debye_Waller",
"Detla_d/d",
"column_F ",
"V_rho",
"density",
"weight",
"nb_atoms",
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->sigma_a) info->sigma_a=atof(parsing[2]);
if (parsing[3] && !info->sigma_a) info->sigma_a=atof(parsing[3]);
if (parsing[4] && !info->sigma_i) info->sigma_i=atof(parsing[4]);
if (parsing[5] && !info->sigma_i) info->sigma_i=atof(parsing[5]);
if (parsing[6]) info->column_order[0]=atoi(parsing[6]);
if (parsing[7]) info->column_order[1]=atoi(parsing[7]);
if (parsing[8]) info->column_order[2]=atoi(parsing[8]);
if (parsing[9]) info->column_order[3]=atoi(parsing[9]);
if (parsing[10]) info->column_order[4]=atoi(parsing[10]);
if (parsing[11]) info->column_order[5]=atoi(parsing[11]);
if (parsing[12]) info->column_order[5]=atoi(parsing[12]);
if (parsing[13]) info->column_order[5]=atoi(parsing[13]);
if (parsing[14]) info->column_order[6]=atoi(parsing[14]);
if (parsing[15] && info->DWfactor<=0) info->DWfactor=atof(parsing[15]);
if (parsing[16] && info->DWfactor<=0) info->DWfactor=atof(parsing[16]);
if (parsing[17] && info->Dd <0) info->Dd =atof(parsing[17]);
if (parsing[18]) info->column_order[7]=atoi(parsing[18]);
if (parsing[19] && !info->V_0) info->V_0 =1/atof(parsing[19]);
if (parsing[20] && !info->rho) info->rho =atof(parsing[20]);
if (parsing[21] && !info->at_weight) info->at_weight =atof(parsing[21]);
if (parsing[22] && info->at_nb <= 1) info->at_nb =atof(parsing[22]);
for (i=0; i<=22; 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;
Table_Info(sTable);
printf("PowderN: %s: Reading %d rows from %s\n",
info->compname, size, SC_file);
/* 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) {
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[list_count-1].j
|| (mult_count == list[list_count].j && i == size-1)) {
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 */
printf("PowderN: %s: File %s done (%i valid lines).\n", info->compname, SC_file, list_count);
Table_Free(&sTable);
info->list = list;
info->count = list_count;
return(list_count);
}
#endif /* !POWDERN_DECL */
%}
DECLARE
%{
struct line_info_struct line_info;
int Nq=0;
double my_s_v2_sum;
double my_a_v, my_inc;
double *w_v,*q_v, *my_s_v2;
char isrect=0;
int columns[8] = format;
%}
INITIALIZE
%{
int i;
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.sigma_a = sigma_abs;
line_info.sigma_i = sigma_inc;
for (i=0; i< 8; i++) line_info.column_order[i] = columns[i];
strncpy(line_info.compname, NAME_CURRENT_COMP, 256);
if (!radius || !yheight) {
if (!xwidth || !yheight || !zthick) exit(fprintf(stderr,"PowderN: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP));
else isrect=1; }
i = read_line_data(reflections, &line_info);
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);
}
if (line_info.V_0 <= 0)
exit(fprintf(stderr,"PowderN: %s: density/unit cell volume is NULL (Vc)\n", NAME_CURRENT_COMP));
if (i) {
L = line_info.list;
Nq = line_info.count;
q_v = malloc(Nq*sizeof(double));
w_v = malloc(Nq*sizeof(double));
my_s_v2 = malloc(Nq*sizeof(double));
if (!q_v || !w_v || !my_s_v2)
exit(fprintf(stderr,"PowderN: %s: ERROR allocating memory (init)\n", NAME_CURRENT_COMP));
for(i=0; i<Nq; i++)
{
my_s_v2[i] = 4*PI*PI*PI*pack*(L[i].DWfactor ? L[i].DWfactor : 1)
/(line_info.V_0*line_info.V_0*V2K*V2K)
*(L[i].j * L[i].F2 / L[i].q)*100; // Factor 100 to convert from barns to fm^2
/* Is not yet divided by v^2 */
/* Squires [3.103] */
my_s_v2_sum+=my_s_v2[i];
q_v[i] = L[i].q*K2V;
w_v[i] = L[i].w;
}
} else exit(fprintf(stderr,
"PowderN: %s: Data file has no valid powder line definition\n"
" Check format definition\n", NAME_CURRENT_COMP));
/* Is not yet divided by v */
my_a_v = pack*line_info.sigma_a/line_info.V_0*2200*100; // Factor 100 to convert from barns to fm^2
my_inc = pack*line_info.sigma_i/line_info.V_0*100; // Factor 100 to convert from barns to fm^2
my_s_v2_sum=0;
printf("PowderN: %s: Vc=%g [Angs] sigma_abs=%g [barn] sigma_inc=%g [barn]\n",
NAME_CURRENT_COMP, line_info.V_0, line_info.sigma_a, line_info.sigma_i);
%}
TRACE
%{
double t0, t1, v, v1,l_full, l, l_1, dt, alpha0, alpha, theta, my_s, my_s_n;
double solid_angle;
double arg, tmp_vx, tmp_vy, tmp_vz, vout_x, vout_y, vout_z,pmul;
int line,flag;
char intersect=0;
if (isrect)
intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zthick);
else
intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight);
if(intersect && t1 >0)
{
if(t0 < 0) t0=0; /* already in sample */
/* Neutron enters at t=t0. */
v = sqrt(vx*vx + vy*vy + vz*vz);
l_full = v * (t1 - t0); /* Length of full path through sample */
dt = rand01()*(t1 - t0); /* Time of scattering */
PROP_DT(dt+t0); /* Point of scattering */
l = v*dt; /* Penetration in sample */
my_s = my_s_v2_sum/(v*v)+my_inc;
/* Total attenuation from scattering */
if (frac <= 0 || (frac < 1 && rand01() >= frac))
{ /* Make coherent scattering event */
/* choose line */
if (Nq > 1) line=floor(Nq*rand01()); /* Select between Nq powder lines */
else line = 0;
if (w_v[line])
arg = q_v[line]*(1+w_v[line]*randnorm())/(2.0*v);
else
arg = q_v[line]/(2.0*v);
my_s_n = my_s_v2[line]/(v*v);
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=PI+(alpha0-1.5*alpha);
} else {
alpha0=alpha0-0.5*alpha;
}
}
else
alpha0 = PI*randpm1();
/* now find a nearly vertical rotation axis:
* (v along Z) x (X axis) -> nearly Y axis
*/
vec_prod(tmp_vx,tmp_vy,tmp_vz, vx,vy,vz, 1,0,0);
/* handle case where v and aim are parallel */
if (!tmp_vx && !tmp_vy && !tmp_vz) { tmp_vx=tmp_vz=0; tmp_vy=1; }
/* v_out = rotate 'v' by 2*theta around tmp_v: Bragg angle */
rotate(vout_x,vout_y,vout_z, vx,vy,vz, 2*theta, tmp_vx,tmp_vy,tmp_vz);
/* tmp_v = rotate v_out by alpha0 around 'v' (Debye-Scherrer cone) */
rotate(tmp_vx,tmp_vy,tmp_vz, vout_x,vout_y,vout_z, alpha0, vx, vy, vz);
vx = tmp_vx;
vy = tmp_vy;
vz = tmp_vz;
flag=0;
if (isrect && !box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zthick)) flag=1;
else if(!isrect && !cylinder_intersect(&t0, &t1, x, y, z,
vx, vx, vx, radius, yheight)) flag=1;
if (flag) {
/* Strange error: did not hit cylinder */
fprintf(stderr, "PowderN: FATAL ERROR: Did not hit sample from inside.\n");
ABSORB;
}
l_1 = v*t1; /* go to exit */
pmul = Nq*l_full*my_s_n*exp(-(my_a_v/v+my_s)*(l+l_1))/(1-frac);
/* Correction in case of d_phi focusing - BUT only when d_phi != 0 */
if (d_phi) pmul *= alpha/PI;
} /* Coherent scattering event */
else
{ /* Make incoherent scattering event */
v = sqrt(vx*vx+vy*vy+vz*vz);
if(d_phi) {
randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,
0, 0, 1,
2*PI, d_phi*DEG2RAD, ROT_A_CURRENT_COMP);
} else {
randvec_target_circle(&vx, &vy, &vz,
&solid_angle, 0, 0, 1, 0);
}
v1 = sqrt(vx*vx+vy*vy+vz*vz);
vx *= v/v1;
vy *= v/v1;
vz *= v/v1;
flag=0;
if (isrect && !box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zthick)) flag=1;
else if(!isrect && !cylinder_intersect(&t0, &t1, x, y, z,
vx, vx, vx, radius, yheight)) flag=1;
if (flag) {
/* Strange error: did not hit cylinder */
fprintf(stderr, "PowderN: FATAL ERROR: Did not hit sample from inside.\n");
ABSORB;
}
l_1 = v*t1;
pmul = l_full*my_inc*exp(-(my_a_v/v+my_s)*(l+l_1))/(frac);
pmul *= solid_angle/(4*PI);
} /* Incoherent scattering event */
p *= pmul;
SCATTER;
} /* else transmit non interacting neutrons */
%}
MCDISPLAY
%{
double h;
h=yheight;
magnify("xyz");
if (!isrect) {
circle("xz", 0, h/2.0, 0, radius);
circle("xz", 0, -h/2.0, 0, radius);
line(-radius, -h/2.0, 0, -radius, +h/2.0, 0);
line(+radius, -h/2.0, 0, +radius, +h/2.0, 0);
line(0, -h/2.0, -radius, 0, +h/2.0, -radius);
line(0, -h/2.0, +radius, 0, +h/2.0, +radius);
} else {
double xmin = -0.5*xwidth;
double xmax = 0.5*xwidth;
double ymin = -0.5*yheight;
double ymax = 0.5*yheight;
double zmin = -0.5*zthick;
double zmax = 0.5*zthick;
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);
}
%}
END