1/***************************************************************************** 2* McXtrace, x-ray tracing package 3* Copyright (C) Risoe National Laboratory, Roskilde, Denmark 4* 5* Component: PowderN 6* 7* %I 8* Written by: P. Willendrup, L. Chapon, K. Lefmann, A.B.Abrahamsen, N.B.Christensen, E.M.Lauridsen. 9* Date: 4.2.98 10* Version: $Revision$ 11* Origin: McStas release 12* Modified by: KL, KN 20.03.98 (rewrite) 13* Modified by: KL, 28.09.01 (two lines) 14* Modified by: KL, 22.05.03 (background) 15* Modified by: KL, PW 01.05.05 (N lines) 16* Modified by: PW, LC 04.10.05 (Merge with Chapon Powder_multi) 17* Modified by: PW, KL 05.06.07 (Concentricity) 18* Modified by: EF, 17.10.08 (added any shape sample geometry) 19* Modified for X-ray use: EK, 28.10.10 20* 21* General powder sample (N lines, single scattering, incoherent scattering) 22* 23* %D 24* General powder sample with 25* many scattering vectors 26* possibility for intrinsic line broadening 27* incoherent backgorund ratio is specified by user. 28* No multiple scattering. No secondary extinction. 29* 30* Based on Powder1/Powder2/Single_crystal. 31* Geometry is a powder filled cylinder or a box. 32* Incoherent scattering is only provided here to account for a background 33* The efficient is highly improved when restricting the vertical scattering 34* range on the Debye-Scherrer cone (with 'd_phi'). 35* The unit cell volume Vc may also be computed when giving the density, 36* the atomic/molecular weight and the number of atoms per unit cell. 37* 38* <b>Sample shape:</b> 39* Sample shape may be a cylinder, a sphere, a box or any other shape. 40* box/plate: xwidth x yheight x zdepth (thickness=0) 41* hollow box/plate:xwidth x yheight x zdepth and thickness>0 42* cylinder: radius x yheight (thickness=0) 43* hollow cylinder: radius x yheight and thickness>0 44* sphere: radius (yheight=0 thickness=0) 45* hollow sphere: radius and thickness>0 (yheight=0) 46* any shape: geometry=OFF_file 47* 48* The complex geometry option handles any closed non-convex polyhedra. 49* It computes the intersection points of the xray with the object 50* transparently, so that it can be used like a regular sample object. 51* It supports the OFF and NOFF file format but not COFF (colored faces). 52* Such files may be generated from XYZ data using: 53* qhull < coordinates.xyz Qx Qv Tv o > geomview.off 54* and viewed with geomview or java -jar jroff.jar (see below). 55* The default size of the object depends of the OFF file data, but its 56* bounding box may be resized using xwidth,yheight and zdepth. 57* 58* If you use this component and produce valuable scientific results, please 59* cite authors with references bellow (in <a href="#links">Links</a>). 60* 61* Example: PowderN(reflections = "c60.lau", d_phi = 15 , radius = 0.01, 62* yheight = 0.05, Vc = 1076.89, Delta_d=0, DW=1, 63* format=Crystallographica) 64* 65* <b>Powder definition file format</b> 66* Powder structure is specified with an ascii data file 'reflections'. 67* The powder data are free-text column based files. 68* Lines begining by '#' are read as comments (ignored) but they may contain 69* the following keywords (in the header): 70* #Vc <value of unit cell volume Vc [Angs^3]> 71* #Debye_Waller <value of Debye-Waller factor DW> 72* #Delta_d/d <value of Detla_d/d width for all lines> 73* These values are not read if entered as component parameters (Vc=...) 74* 75* The signification of the columns in the numerical block may be 76* set using the 'format' parameter. Built-in formats are: 77* format=Crystallographica 78* format=Fullprof 79* format=Lazy 80* and these specifications it is important NOT to use quotes, as shown. 81* 82* An other possibility to define other formats is to set directly 83* the signification of the columns as a vector of indexes in the order 84* format={j,d,F2,DW,Delta_d/d,1/2d,q,F} 85* Signification of the symbols is given below. Indexes start at 1. 86* Indexes with zero means that the column is not present, so that: 87* Crystallographica={ 4,5,7,0,0,0,0,0 } 88* Fullprof ={ 4,0,8,0,0,5,0,0 } 89* Lazy ={17,6,0,0,0,0,0,13} 90* Here again, NO quotes should be around the 'format' value. 91* 92* At last, the format may be overridden by direct definition of the 93* column indexes in the file itself by using the following keywords 94* in the header (e.g. '#column_j 4'): 95* #column_j <index of the multiplicity 'j' column> 96* #column_d <index of the d-spacing 'd' column [Angs]> 97* #column_F2 <index of the squared str. factor '|F|^2' column [b]> 98* #column_F <index of the structure factor norm '|F|' column> 99* #column_DW <index of the Debye-Waller factor 'DW' column> 100* #column_Dd <index of the relative line width Delta_d/d 'Dd' column> 101* #column_inv2d <index of the 1/2d=sin(theta)/lambda 'inv2d' column> 102* #column_q <index of the scattering wavevector 'q' column [Angs-1]> 103* 104* <b>Concentricity</b> 105* 106* PowderN assumes 'concentric' shape, i.e. can contain other components inside its 107* optional inner hollow. Example, Sample in Al cryostat: 108* 109* 110* COMPONENT Cryo = PowderN(reflections="Al.laz", radius = 0.01, thickness = 0.001, 111* concentric = 1, p_interact=0.1) 112* AT (0,0,0) RELATIVE Somewhere 113* 114* COMPONENT Sample = some_other_component(with geometry FULLY enclosed in the hollow) 115* AT (0,0,0) RELATIVE Somewhere 116* 117* COMPONENT Cryo2 = COPY(Cryo)(concentric = 0) 118* AT (0,0,0) RELATIVE Somewhere 119* 120* 121* (The second instance of the cryostat component can also be written out completely 122* using PowderN(...). In both cases, this second instance needs concentric = 0.) 123* 124* %P 125* INPUT PARAMETERS 126* radius: Outer radius of sample in (x,z) plane [m] 127* xwidth: Horiz. dimension of sample, as a width [m] 128* yheight: Height of sample y direction [m] 129* zdepth: Depth of box sample [m] 130* thickness:Thickness of hollow sample. 131* Negative value extends the hollow volume outside of the box/cylinder. [m] 132* reflections: Input file for reflections. 133* Use only incoherent scattering if NULL or "" [string] 134* 135* Optional parameters: 136* d_phi: Angle corresponding to the vertical angular range 137* to focus to, e.g. detector height. 0 for no focusing [deg,0-180] 138* pack: Packing factor [1] 139* Delta_d: Global relative Delta_d/d spreading when the 'w' column 140* is not available. Use 0 if ideal. [Angs] 141* format: Name of the format, or list of column indexes 142* (see Description). [no quotes] 143* p_inc: Fraction of incoherently scattered xrays [1] 144* p_transmit: Fraction of transmitted (only attenuated) xrays [1] 145* p_interact: Fraction of events interacting with sample, e.g. 1-p_transmit-p_inc [1] 146* concentric: Indicate that this component has a hollow geometry 147* and may contain other components. It should then be duplicated 148* after the inside part (only for box, cylinder, sphere) [1] 149* geometry: Name of an Object File Format (OFF) file for complex geometry. 150* The OFF file may be generated from XYZ coordinates using qhull/powercrust [str] 151* barns: Flag to indicate if |F|^2 from 'reflections' is in barns or fm^2 152* (barns=1 for laz, barns=0 for lau type files).[1] 153* Vc: Volume of unit cell=nb atoms per cell/density of atoms [AA^3] 154* DW: Global Debey-Waller factor when the 'DW' column 155* is not available. Use 1 if included in F2 [1] 156* weight: Atomic/molecular weight of material [g/mol] 157* density: Density of material. rho=density/weight/1e24*N_A. [g/cm^3] 158* nb_atoms: Number of sub-unit per unit cell, that is ratio of sigma for chemical formula to sigma per unit cell [1] 159* 160* %L 161* “Validation of a realistic powder sample using data from DMC at PSI” Willendrup P, Filges U, Keller L, Farhi E, Lefmann K, Physica B-Cond Matt 385 (2006) 1032. 162* %L 163* See also: Powder1, Powder2 and PowderN 164* %L 165* See <a href="http://icsd.ill.fr">ICSD</a> Inorganic Crystal Structure Database 166* %L 167* Cross sections for single elements: http://www.ncnr.nist.gov/resources/n-lengths/ 168* %L 169* Cross sections for compounds: http://www.ncnr.nist.gov/resources/sldcalc.html 170* %L 171* Web Elements http://www.webelements.com/ 172* %L 173* Fullprof powder refinement: http://www.ill.eu/sites/fullprof/index.html 174* %L 175* Crystallographica software: http://www.crystallographica.com/ 176* %L 177* Geomview and Object File Format (OFF) <http://www.geomview.org> 178* %L 179* Java version of Geomview (display only) jroff.jar <http://www.holmes3d.net/graphics/roffview/> 180* %L 181* Powercrust/qhull <http://qhull.org> 182* 183* %E 184*****************************************************************************/ 185DEFINE COMPONENT PowderN 186DEFINITION PARAMETERS (format=Undefined) 187 SETTING PARAMETERS (string reflections=0, string material=0, string geometry=0, 188 radius=0, yheight=0, xwidth=0, zdepth=0, thickness=0, 189 pack=1, Vc=0, Delta_d=0, p_inc=0.1, p_transmit=0.1, 190 DW=0, nb_atoms=1, d_phi=0, p_interact=0, 191 concentric=0, density=0, weight=0, barns=1) 192OUTPUT PARAMETERS (line_info, Nq, my_s_k2, XsectionFactor, 193 my_s_k2_sum, my_a_k, my_inc, my_q, my_w, columns, 194 radius_i, xwidth_i, yheight_i, zdepth_i, offdata) 195/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 196SHARE 197%{ 198 /* used for reading data table from file */ 199 %include "read_table-lib" 200 %include "interoff-lib" 201/* Declare structures and functions only once in each instrument. */ 202#ifndef POWDERN_DECL 203#define POWDERN_DECL 204/* format definitions in the order {j d F2 DW Dd inv2d q F} */ 205#ifndef Crystallographica 206#define Crystallographica { 4,5,7,0,0,0,0,0 } 207#define Fullprof { 4,0,8,0,0,5,0,0 } 208#define Lazy {17,6,0,0,0,0,0,13 } 209#define Undefined { 0,0,0,0,0,0,0,0 } 210#endif 211 212 struct line_data 213 { 214 double F2; /* Value of structure factor */ 215 double q; /* Qvector */ 216 int j; /* Multiplicity */ 217 double DWfactor; /* Debye-Waller factor */ 218 double w; /* Intrinsic line width */ 219 }; 220 221 struct abs_data 222 { 223 double E; /*energy (in keV)*/ 224 double k; /*wavenumber corresponding to E*/ 225 double sigma_a; /*absorption cross section for energy E*/ 226 double mu; /*absoprtion coefficient for the energy E*/ 227 }; 228 struct line_info_struct 229 { 230 struct line_data *list; /* Reflection array */ 231 int count; /* Number of reflections */ 232 double Dd; 233 double DWfactor; 234 double V_0; 235 double rho; 236 double at_weight; 237 double at_nb; 238 double sigma_a; 239 double sigma_i; 240 char compname[256]; 241 double flag_barns; 242 int shape; /* 0 cylinder, 1 box, 2 sphere, 3 OFF file */ 243 int column_order[8]; /* column signification */ 244 int flag_warning; 245 }; 246 247 off_struct offdata; 248 249 int read_abs_data(char *ABS_file, struct abs_data **abs) 250 { 251 t_Table table; 252 char **parsing; 253 int status,i; 254 255 if (!ABS_file || !strlen(ABS_file) || !strcmp(ABS_file, "NULL")) { 256 fprintf(stderr,"Warning: material file (%s) not found.\n",ABS_file); 257 *abs=calloc(2,sizeof(struct abs_data)); 258 (*abs)[1].E=-1; 259 (*abs)[1].k=-1; 260 (*abs)[1].mu=-1; 261 }else{ 262 if ( (status=Table_Read(&table,ABS_file,0))==-1){ 263 fprintf(stderr,"Error: %s Could not parse file \"%s\"\n",NAME_CURRENT_COMP,ABS_file); 264 exit(-1); 265 } 266 parsing=Table_ParseHeader(table.header,"Z","A[r]","rho",NULL); 267 *abs=calloc(table.rows+1,sizeof(struct abs_data)); 268 for (i=0;i<table.rows;i++){ 269 (*abs)[i].E=table.data[i*table.columns]; 270 (*abs)[i].k=(*abs)[i].E*E2K; 271 // get from cox_76 abs_info.data[i].sigma_a=aba_info.data[i]* 272 (*abs)[i].mu=table.data[i*table.columns + 7]; 273 } 274 (*abs)[i].E=-1; 275 (*abs)[i].k=-1; 276 (*abs)[i].mu=-1; 277 Table_Free(&table); 278 return 1; 279 } 280 } 281 282 int read_line_data(char *SC_file, struct line_info_struct *info) 283 { 284 struct line_data *list = NULL; 285 int size = 0; 286 t_Table sTable; /* sample data table structure from SC_file */ 287 int i=0; 288 int mult_count =0; 289 char flag=0; 290 double q_count=0, j_count=0, F2_count=0; 291 char **parsing; 292 int list_count=0; 293 294 if (!SC_file || !strlen(SC_file) || !strcmp(SC_file, "NULL")) { 295 printf("PowderN: %s: Using incoherent elastic scattering only\n", 296 info->compname); 297 info->count = 0; 298 return(0); 299 } 300 Table_Read(&sTable, SC_file, 1); /* read 1st block data from SC_file into sTable*/ 301 302 /* parsing of header */ 303 parsing = Table_ParseHeader(sTable.header, 304 "Vc","V_0", 305 "column_j", 306 "column_d", 307 "column_F2", 308 "column_DW", 309 "column_Dd", 310 "column_inv2d", "column_1/2d", "column_sintheta/lambda", 311 "column_q", /* 10 */ 312 "DW", "Debye_Waller", 313 "Delta_d/d", 314 "column_F ", 315 "V_rho", 316 "density", 317 "weight", 318 "nb_atoms","multiplicity", 319 NULL); 320 321 if (parsing) { 322 if (parsing[0] && !info->V_0) info->V_0 =atof(parsing[0]); 323 if (parsing[1] && !info->V_0) info->V_0 =atof(parsing[1]); 324 if (parsing[2]) info->column_order[0]=atoi(parsing[2]); 325 if (parsing[3]) info->column_order[1]=atoi(parsing[3]); 326 if (parsing[4]) info->column_order[2]=atoi(parsing[4]); 327 if (parsing[5]) info->column_order[3]=atoi(parsing[5]); 328 if (parsing[6]) info->column_order[4]=atoi(parsing[6]); 329 if (parsing[7]) info->column_order[5]=atoi(parsing[7]); 330 if (parsing[8]) info->column_order[5]=atoi(parsing[8]); 331 if (parsing[9]) info->column_order[5]=atoi(parsing[9]); 332 if (parsing[10]) info->column_order[6]=atoi(parsing[10]); 333 if (parsing[11] && info->DWfactor<=0) info->DWfactor=atof(parsing[11]); 334 if (parsing[12] && info->DWfactor<=0) info->DWfactor=atof(parsing[12]); 335 if (parsing[13] && info->Dd <0) info->Dd =atof(parsing[13]); 336 if (parsing[14]) info->column_order[7]=atoi(parsing[14]); 337 if (parsing[15] && !info->V_0) info->V_0 =1/atof(parsing[15]); 338 if (parsing[16] && !info->rho) info->rho =atof(parsing[16]); 339 if (parsing[27] && !info->at_weight) info->at_weight =atof(parsing[17]); 340 if (parsing[18] && info->at_nb <= 1) info->at_nb =atof(parsing[18]); 341 if (parsing[19] && info->at_nb <= 1) info->at_nb =atof(parsing[19]); 342 for (i=0; i<=19; i++) if (parsing[i]) free(parsing[i]); 343 free(parsing); 344 } 345 346 if (!sTable.rows) 347 exit(fprintf(stderr, "PowderN: %s: Error: The number of rows in %s " 348 "should be at least %d\n", info->compname, SC_file, 1)); 349 else size = sTable.rows; 350 Table_Info(sTable); 351 printf("PowderN: %s: Reading %d rows from %s\n", 352 info->compname, size, SC_file); 353 354 if (info->column_order[0] == 4 && info->flag_barns !=0) 355 printf("PowderN: %s: Powder file probably of type Crystallographica/Fullprof (lau)\n" 356 "WARNING: but F2 unit is set to barns=1 (barns). Intensity might be 100 times too high.\n", 357 info->compname); 358 if (info->column_order[0] == 17 && info->flag_barns == 0) 359 printf("PowderN: %s: Powder file probably of type Lazy Pulver (laz)\n" 360 "WARNING: but F2 unit is set to barns=0 (fm^2). Intensity might be 100 times too low.\n", 361 info->compname); 362 /* allocate line_data array */ 363 list = (struct line_data*)malloc(size*sizeof(struct line_data)); 364 365 for (i=0; i<size; i++) 366 { 367 /* printf("Reading in line %i\n",i);*/ 368 double j=0, d=0, w=0, q=0, DWfactor=0, F2=0; 369 int index; 370 371 if (info->Dd >= 0) w = info->Dd; 372 if (info->DWfactor > 0) DWfactor = info->DWfactor; 373 374 /* get data from table using columns {j d F2 DW Dd inv2d q F} */ 375 /* column indexes start at 1, thus need to substract 1 */ 376 if (info->column_order[0] >0) 377 j = Table_Index(sTable, i, info->column_order[0]-1); 378 if (info->column_order[1] >0) 379 d = Table_Index(sTable, i, info->column_order[1]-1); 380 if (info->column_order[2] >0) 381 F2 = Table_Index(sTable, i, info->column_order[2]-1); 382 if (info->column_order[3] >0) 383 DWfactor = Table_Index(sTable, i, info->column_order[3]-1); 384 if (info->column_order[4] >0) 385 w = Table_Index(sTable, i, info->column_order[4]-1); 386 if (info->column_order[5] >0) 387 { d = Table_Index(sTable, i, info->column_order[5]-1); 388 d = (d > 0? 1/d/2 : 0); } 389 if (info->column_order[6] >0) 390 { q = Table_Index(sTable, i, info->column_order[6]-1); 391 d = (q > 0 ? 2*PI/q : 0); } 392 if (info->column_order[7] >0 && !F2) 393 { F2 = Table_Index(sTable, i, info->column_order[7]-1); F2 *= F2; } 394 395 /* assign and check values */ 396 j = (j > 0 ? j : 0); 397 q = (d > 0 ? 2*PI/d : 0); /* this is q */ 398 DWfactor = (DWfactor > 0 ? DWfactor : 1); 399 w = (w>0 ? w : 0); /* this is q and d relative spreading */ 400 F2 = (F2 >= 0 ? F2 : 0); 401 if (j == 0 || q == 0) { 402 printf("PowderN: %s: line %i has invalid definition\n" 403 " (mult=0 or q=0 or d=0)\n", info->compname, i); 404 continue; 405 } 406 list[list_count].j = j; 407 list[list_count].q = q; 408 list[list_count].DWfactor = DWfactor; 409 list[list_count].w = w; 410 list[list_count].F2= F2; 411 412 /* adjust multiplicity if j-column + multiple d-spacing lines */ 413 /* if d = previous d, increase line duplication index */ 414 if (!q_count) q_count = q; 415 if (!j_count) j_count = j; 416 if (!F2_count) F2_count = F2; 417 if (fabs(q_count-q) < 0.0001*fabs(q) 418 && fabs(F2_count-F2) < 0.0001*fabs(F2) && j_count == j) { 419 mult_count++; flag=0; } 420 else flag=1; 421 if (i == size-1) flag=1; 422 /* else if d != previous d : just passed equivalent lines */ 423 if (flag) { 424 if (i == size-1) list_count++; 425 /* if duplication index == previous multiplicity */ 426 /* set back multiplicity of previous lines to 1 */ 427 if (mult_count && list_count>0 && 428 (mult_count == list[list_count-1].j 429 || (mult_count == list[list_count].j && i == size-1)) ) { 430 printf("PowderN: %s: Set multiplicity to 1 for lines [%i:%i]\n" 431 " (d-spacing %g is duplicated %i times)\n", 432 info->compname, list_count-mult_count, list_count-1, list[list_count-1].q, mult_count); 433 for (index=list_count-mult_count; index<list_count; list[index++].j = 1); 434 mult_count = 1; 435 q_count = q; 436 j_count = j; 437 F2_count = F2; 438 } 439 if (i == size-1) list_count--; 440 flag=0; 441 } 442 list_count++; 443 } /* end for */ 444 printf("PowderN: %s: File %s done (%i valid lines).\n", info->compname, SC_file, list_count); 445 Table_Free(&sTable); 446 info->list = list; 447 info->count = list_count; 448 449 return(list_count); 450 } 451#endif /* !POWDERN_DECL */ 452 453%} 454 455DECLARE 456%{ 457 struct line_info_struct line_info; 458 struct abs_data *a_info; 459 int Nq=0; 460 double my_s_k2_sum=0; 461 double my_a_k=0, my_inc=0; 462 double *my_w,*my_q, *my_s_k2; 463 int columns[8] = format; 464 double XsectionFactor; 465 double radius_i=0, xwidth_i=0, yheight_i=0, zdepth_i=0; 466 467 off_struct offdata; 468%} 469INITIALIZE 470%{ 471 int i=0; 472 struct line_data *L; 473 line_info.Dd = Delta_d; 474 line_info.DWfactor = DW; 475 line_info.V_0 = Vc; 476 line_info.rho = density; 477 line_info.at_weight= weight; 478 line_info.at_nb = nb_atoms; 479 line_info.flag_barns=barns; 480 line_info.shape = 0; 481 line_info.flag_warning=0; 482 for (i=0; i< 8; i++) line_info.column_order[i] = columns[i]; 483 strncpy(line_info.compname, NAME_CURRENT_COMP, 256); 484 485 line_info.shape=-1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape */ 486 if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) { 487 if (off_init(geometry, xwidth, yheight, zdepth, 0, &offdata)) { 488 line_info.shape=3; thickness=0; concentric=0; 489 } 490 } 491 else if (xwidth && yheight && zdepth) line_info.shape=1; /* box */ 492 else if (radius > 0 && yheight) line_info.shape=0; /* cylinder */ 493 else if (radius > 0 && !yheight) line_info.shape=2; /* sphere */ 494 495 if (line_info.shape < 0) 496 exit(fprintf(stderr,"PowderN: %s: sample has invalid dimensions.\n" 497 "ERROR Please check parameter values (xwidth, yheight, zdepth, radius).\n", NAME_CURRENT_COMP)); 498 if (thickness) { 499 if (radius && (radius < fabs(thickness))) { 500 fprintf(stderr,"PowderN: %s: hollow sample thickness is larger than its volume (sphere/cylinder).\n" 501 "WARNING Please check parameter values. Using bulk sample (thickness=0).\n", NAME_CURRENT_COMP); 502 thickness=0; 503 } 504 else if (!radius && (xwidth < 2*fabs(thickness) || yheight < 2*fabs(thickness) || zdepth < 2*fabs(thickness))) { 505 fprintf(stderr,"PowderN: %s: hollow sample thickness is larger than its volume (box).\n" 506 "WARNING Please check parameter values.\n", NAME_CURRENT_COMP); 507 } 508 } 509 510 if (concentric && thickness==0) { 511 printf("PowderN: %s:Can not use concentric mode\n" 512 "WARNING on non hollow shape. Ignoring.\n", 513 NAME_CURRENT_COMP); 514 concentric=0; 515 } 516 517 if (thickness>0) { 518 if (radius>thickness) { 519 radius_i=radius-thickness; 520 } else { 521 if (xwidth>2*thickness) xwidth_i =xwidth -2*thickness; 522 if (yheight>2*thickness) yheight_i=yheight-2*thickness; 523 if (zdepth>2*thickness) zdepth_i =zdepth -2*thickness; 524 } 525 } else if (thickness<0) { 526 thickness = fabs(thickness); 527 if (radius) { 528 radius_i=radius; 529 radius=radius_i+thickness; 530 } else { 531 xwidth_i =xwidth; 532 yheight_i=yheight; 533 zdepth_i =zdepth; 534 xwidth =xwidth +2*thickness; 535 yheight =yheight+2*thickness; 536 zdepth =zdepth +2*thickness; 537 } 538 } 539 540 if (!yheight_i) { 541 yheight_i = yheight; 542 } 543 if (p_interact) { 544 if (p_interact < p_inc) { double tmp=p_interact; p_interact=p_inc; p_inc=tmp; } 545 p_transmit = 1-p_interact-p_inc; 546 } 547 548 if (p_inc + p_transmit > 1) { 549 fprintf(stderr,"PowderN: %s: You have requested an unmeaningful choice of the 'p_inc' and 'p_transmit' parameters (sum is %g, exeeding 1). Fixing.\n", 550 NAME_CURRENT_COMP, p_inc+p_transmit); 551 if (p_inc > p_transmit) p_transmit=1-2*p_inc; 552 else p_transmit=1-2*p_inc; 553 } else if (p_inc + p_transmit == 1) { 554 printf("PowderN: %s: You have requested all xrays be attenuated\n" 555 "WARNING or incoherently scattered!\n", NAME_CURRENT_COMP); 556 } 557 558 if (concentric) { 559 printf("PowderN: %s WARNING: Concentric mode - remember to include the 'opposite' copy of this component !\n (The equivalent, 'opposite' comp should have concentric=0)\n", NAME_CURRENT_COMP); 560 if (p_transmit == 0) { 561 printf("PowderN: %s: Concentric mode and p_transmit==0 !?\n" 562 "WARNING Don't you want any transmitted xrays?\n", NAME_CURRENT_COMP); 563 } 564 } 565 566 if (reflections && strlen(reflections) && strcmp(reflections, "NULL") && strcmp(reflections, "0")) { 567 i = read_line_data(reflections, &line_info); 568 if (i == 0) 569 exit(fprintf(stderr,"PowderN: %s: reflection file %s is not valid.\n" 570 "ERROR Please check file format (laz or lau).\n", NAME_CURRENT_COMP, reflections)); 571 } 572 573 /* compute the scattering unit density from material weight and density */ 574 /* the weight of the scattering element is the chemical formula molecular weight 575 * times the nb of chemical formulae in the scattering element (nb_atoms) */ 576 if (!line_info.V_0 && line_info.at_nb > 0 577 && line_info.at_weight > 0 && line_info.rho > 0) { 578 /* molar volume [cm^3/mol] = weight [g/mol] / density [g/cm^3] */ 579 /* atom density per Angs^3 = [mol/cm^3] * N_Avogadro *(1e-8)^3 */ 580 line_info.V_0 = line_info.at_nb 581 /(line_info.rho/line_info.at_weight/1e24*6.02214199e23); 582 } 583 584 /* the scattering unit cross sections are the chemical formula onces 585 * times the nb of chemical formulae in the scattering element */ 586 if (line_info.at_nb > 0) { 587 line_info.sigma_a *= line_info.at_nb; line_info.sigma_i *= line_info.at_nb; 588 } 589 590 if (line_info.V_0 <= 0) 591 fprintf(stderr,"PowderN: %s: density/unit cell volume is NULL (Vc). Unactivating component.\n", NAME_CURRENT_COMP); 592 593 if (line_info.V_0 && p_inc && !line_info.sigma_i) { 594 printf("PowderN %s WARNING: You have requested statistics for incoherent scattering but not defined sigma_inc!\n", NAME_CURRENT_COMP); 595 } 596 597 if (line_info.flag_barns) { /* Factor 100 to convert from barns to fm^2 */ 598 XsectionFactor = 100; 599 } else { 600 XsectionFactor = 1; 601 } 602 603 if (line_info.V_0 && i) { 604 L = line_info.list; 605 606 Nq = line_info.count; 607 my_q = malloc(Nq*sizeof(double)); 608 my_w = malloc(Nq*sizeof(double)); 609 my_s_k2 = malloc(Nq*sizeof(double)); 610 if (!my_q || !my_w || !my_s_k2) 611 exit(fprintf(stderr,"PowderN: %s: ERROR allocating memory (init)\n", NAME_CURRENT_COMP)); 612 for(i=0; i<Nq; i++) 613 { 614 my_s_k2[i] = 4*PI*PI*PI*pack*(L[i].DWfactor ? L[i].DWfactor : 1) 615 /(line_info.V_0*line_info.V_0) 616 *(L[i].j * L[i].F2 / L[i].q)*XsectionFactor; 617 /* Is not yet divided by k^2 */ 618 /* Squires [3.103] */ 619 my_q[i] = L[i].q; 620 my_w[i] = L[i].w; 621 } 622 } 623 if (!read_abs_data(material,&a_info)){ 624 fprintf(stderr,"PowderN: %s Error reading absorption data from file %s - will proceed without absorption.\n",NAME_CURRENT_COMP,material); 625 } 626// if (line_info.V_0) { 627// /* Is not yet divided by v */ 628// my_abs = pack*line_info.sigma_a/line_info.V_0*2200*100; // Factor 100 to convert from barns to fm^2 629// my_inc = pack*line_info.sigma_i/line_info.V_0*100; // Factor 100 to convert from barns to fm^2 630// 631// printf("PowderN: %s: Vc=%g [Angs] sigma_abs=%g [barn] sigma_inc=%g [barn] reflections=%s\n", 632// NAME_CURRENT_COMP, line_info.V_0, line_info.sigma_a, line_info.sigma_i, reflections && strlen(reflections) ? reflections : "NULL"); 633// } 634 635%} 636TRACE 637%{ 638 double l0, l1, l2, l3, k, k1,l_full, l, l_1, dl, alpha0, alpha, theta, my_s, my_s_n; 639 double solid_angle, type; 640 double arg, tmp_kx, tmp_ky, tmp_kz, kout_x, kout_y, kout_z, nx, ny, nz, my_abs, pmul=1; 641 int line; 642 char intersect=0; 643 char intersecti=0; 644 645 if (line_info.V_0 && (Nq || my_inc)) { 646 if (line_info.shape == 1) { 647 intersect = box_intersect(&l0, &l3, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); 648 intersecti = box_intersect(&l1, &l2, x, y, z, kx, ky, kz, xwidth_i, yheight_i, zdepth_i); 649 } else if (line_info.shape == 0) { 650 intersect = cylinder_intersect(&l0, &l3, x, y, z, kx, ky, kz, radius, yheight); 651 intersecti = cylinder_intersect(&l1, &l2, x, y, z, kx, ky, kz, radius_i, yheight_i); 652 } else if (line_info.shape == 2) { 653 intersect = sphere_intersect (&l0, &l3, x,y,z, kx,ky,kz, radius); 654 intersecti = sphere_intersect (&l1, &l2, x,y,z, kx,ky,kz, radius_i); 655 } else if (line_info.shape == 3) { 656 intersect = off_intersect (&l0, &l3, NULL, NULL, x,y,z, kx,ky,kz, offdata); 657 intersecti = 0; 658 } 659 } 660 661 if(intersect && l3 >0) { 662 663 if (concentric) { 664 /* Set up for concentric case */ 665 /* 'Remove' the backside of this comp */ 666 if (!intersecti) { 667 l1 = (l3 + l0) /2; 668 } 669 l2 = l1; 670 l3 = l1; 671 dl = -1.0*rand01(); /* In case of scattering we will scatter on 'forward' part of sample */ 672 } else { 673 if (!intersecti) { 674 l1 = (l3 + l0) /2; 675 l2 = l1; 676 } 677 dl = randpm1(); /* Possibility to scatter at all points in line of sight */ 678 } 679 680 /* X-ray enters at t=l0. */ 681 if(l0 < 0) l0=0; /* already in sample */ 682 if(l1 < 0) l1=0; /* already in inner hollow */ 683 if(l2 < 0) l2=0; /* already past inner hollow */ 684 k = sqrt(kx*kx + ky*ky + kz*kz); 685 l_full =l3 - l2 + l1 - l0; 686 687 /* Calculate total scattering cross section at relevant velocity */ 688 my_s_k2_sum=0; 689 for(line=0; line<Nq; line++) { 690 if (my_q[line] <= 2*k) { 691 my_s_k2_sum+=my_s_k2[line]; 692 } 693 } 694 695 696 697 698 if (l3 < 0) { 699 l3=0; /* Already past sample?! */ 700 if (line_info.flag_warning < 100) 701 printf("PowderN: %s: Warning: xray has already passed us? (Skipped).\n" 702 " In concentric geometry, this may be caused by a missing concentric=0 option in 2nd enclosing instance.\n", NAME_CURRENT_COMP); 703 line_info.flag_warning++; 704 } else { 705 if (dl<0) { /* Calculate scattering point position */ 706 dl = fabs(dl)*(l1 - l0); /* 'Forward' part */ 707 } else { 708 dl = dl * (l3 - l2) + (l2-l0) ; /* Possibly also 'backside' part */ 709 } 710 711 my_s = my_s_k2_sum/(k*k); 712 /* Total attenuation from scattering */ 713 714 /* Total attenuation from absorption */ 715 int ii=0; 716 double delta; 717 while(a_info[ii].k!=-1 && k>a_info[ii].k){ 718 ii++; 719 } 720 if (a_info[ii].k==-1) { 721 my_abs=density*a_info[ii-1].mu; 722 }else{ 723 delta=(k-a_info[ii-1].k)/(a_info[ii].k-a_info[ii-1].k); 724 my_abs=density*( (1-delta)*a_info[ii-1].mu + delta*a_info[ii].mu ); 725 } 726 //printf("myabs=%g, k=%g\n",my_abs,k); 727 728 type = rand01(); 729 p_inc=0; /*is turned off till we fix the Compton scattering code*/ 730 /* How to handle this one? Transmit (1) / Incoherent (2) / Coherent (3) ? */ 731 if (type < p_transmit) { 732 type = 1; 733 l = l_full; /* Passing through, full length */ 734 PROP_DL(l3); 735// } else if (type >= p_transmit && type < (p_transmit + p_inc)) { 736// type = 2; 737// l = v*dt; /* Penetration in sample */ 738// PROP_DL(dl+l0); /* Point of scattering */ 739// SCATTER; 740 } else if (type >= p_transmit + p_inc) { 741 type = 3; 742 PROP_DL(dl+l0); /* Point of scattering */ 743 SCATTER; 744 } else { 745 exit(fprintf(stderr,"PowderN %s: DEAD - this shouldn't happen!\n", NAME_CURRENT_COMP)); 746 } 747 748 if (type == 3) { /* Make coherent scattering event */ 749 if (line_info.count > 0) { 750 /* choose line */ 751 if (Nq > 1) line=floor(Nq*rand01()); /* Select between Nq powder lines */ 752 else line = 0; 753 if (my_w[line]) 754 arg = my_q[line]*(1+my_w[line]*randnorm())/(2.0*k); 755 else 756 arg = my_q[line]/(2.0*k); 757 my_s_n = my_s_k2[line]/(k*k); 758 if(fabs(arg) > 1) 759 ABSORB; /* No bragg scattering possible*/ 760 theta = asin(arg); /* Bragg scattering law */ 761 762 /* Choose point on Debye-Scherrer cone */ 763 if (d_phi) 764 { /* relate height of detector to the height on DS cone */ 765 arg = sin(d_phi*DEG2RAD/2)/sin(2*theta); 766 /* If full Debye-Scherrer cone is within d_phi, don't focus */ 767 if (arg < -1 || arg > 1) d_phi = 0; 768 /* Otherwise, determine alpha to rotate from scattering plane 769 into d_phi focusing area*/ 770 else alpha = 2*asin(arg); 771 } 772 if (d_phi) { 773 /* Focusing */ 774 alpha = fabs(alpha); 775 /* Trick to get scattering for pos/neg theta's */ 776 alpha0= 2*rand01()*alpha; 777 if (alpha0 > alpha) { 778 alpha0=PI+(alpha0-1.5*alpha); 779 } else { 780 alpha0=alpha0-0.5*alpha; 781 } 782 } 783 else 784 alpha0 = PI*randpm1(); 785 786 /* now find a nearly vertical rotation axis: 787 * Either 788 * (k along Z) x (X axis) -> nearly Y axis 789 * Or 790 * (k along X) x (Z axis) -> nearly Y axis 791 */ 792 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))) { 793 nx = 1; ny = 0; nz = 0; 794 } else { 795 nx = 0; ny = 0; nz = 1; 796 } 797 vec_prod(tmp_kx,tmp_ky,tmp_kz, kx,ky,kz, nx,ny,nz); 798 799 /* k_out = rotate 'k' by 2*theta around tmp_k: Bragg angle */ 800 rotate(kout_x,kout_y,kout_z, kx,ky,kz, 2*theta, tmp_kx,tmp_ky,tmp_kz); 801 802 /* tmp_k = rotate k_out by alpha0 around 'k' (Debye-Scherrer cone) */ 803 rotate(tmp_kx,tmp_ky,tmp_kz, kout_x,kout_y,kout_z, alpha0, kx, ky, kz); 804 kx = tmp_kx; 805 ky = tmp_ky; 806 kz = tmp_kz; 807 /* Since now scattered and new direction given, calculate path to exit */ 808 if (line_info.shape == 1) { 809 intersect = box_intersect(&l0, &l3, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); 810 intersecti = box_intersect(&l1, &l2, x, y, z, kx, ky, kz, xwidth_i, yheight_i, zdepth_i); 811 } else if (line_info.shape == 0) { 812 intersect = cylinder_intersect(&l0, &l3, x, y, z, kx, ky, kz, radius, yheight); 813 intersecti = cylinder_intersect(&l1, &l2, x, y, z, kx, ky, kz, radius_i, yheight_i); 814 } else if (line_info.shape == 2) { 815 intersect = sphere_intersect (&l0, &l3, x,y,z, kx,ky,kz, radius); 816 intersecti = sphere_intersect (&l1, &l2, x,y,z, kx,ky,kz, radius_i); 817 } else if (line_info.shape == 3) { 818 intersect = off_intersect (&l0, &l3, NULL, NULL, x,y,z, kx,ky,kz, offdata); 819 intersecti = 0; 820 } 821 822 if (!intersect) { 823 /* Strange error: did not hit cylinder */ 824 if (line_info.flag_warning < 100) 825 fprintf(stderr, "PowderN: %s: FATAL ERROR: Did not hit sample from inside (coh).\n", NAME_CURRENT_COMP); 826 line_info.flag_warning++; 827 ABSORB; 828 } 829 830 if (!intersecti) { 831 l1 = (l3 + l0) /2; 832 l2 = l1; 833 } 834 835 if (concentric && intersecti) { 836 /* In case of concentricity, 'remove' backward wall of sample */ 837 l2 = l1; 838 l3 = l1; 839 } 840 841 if(l0 < 0) l0=0; /* already in sample */ 842 if(l1 < 0) l1=0; /* already in inner hollow */ 843 if(l2 < 0) l2=0; /* already past inner hollow */ 844 845 846 l_1 = l3 - l2 + l1 - l0; /* Length to exit */ 847 848 pmul = Nq*l_full*my_s_n*exp(-(my_abs+my_s)*(l+l_1))/(1-(p_inc+p_transmit)); 849 /* Correction in case of d_phi focusing - BUT only when d_phi != 0 */ 850 if (d_phi) pmul *= alpha/PI; 851 } /* else transmit <-- No powder lines in file */ 852 } /* Coherent scattering event */ 853// else if (type == 2) { /* Make incoherent scattering event */ 854 /*should be replaced by Compton scattering and/or TDS*/ 855// if(d_phi) { 856// randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle, 857// 0, 0, 1, 858// 2*PI, d_phi*DEG2RAD, ROT_A_CURRENT_COMP); 859// } else { 860// randvec_target_circle(&vx, &vy, &vz, 861// &solid_angle, 0, 0, 1, 0); 862// } 863// v1 = sqrt(vx*vx+vy*vy+vz*vz); 864// vx *= v/v1; 865// vy *= v/v1; 866// vz *= v/v1; 867// 868// /* Since now scattered and new direction given, calculate path to exit */ 869// if (line_info.shape == 1) { 870// intersect = box_intersect(&l0, &l3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); 871// intersecti = box_intersect(&l1, &l2, x, y, z, vx, vy, vz, xwidth_i, yheight_i, zdepth_i); 872// } else if (line_info.shape == 0) { 873// intersect = cylinder_intersect(&l0, &l3, x, y, z, vx, vy, vz, radius, yheight); 874// intersecti = cylinder_intersect(&l1, &l2, x, y, z, vx, vy, vz, radius_i, yheight_i); 875// } else if (line_info.shape == 2) { 876// intersect = sphere_intersect (&l0, &l3, x,y,z, vx,vy,vz, radius); 877// intersecti = sphere_intersect (&l1, &l2, x,y,z, vx,vy,vz, radius_i); 878// } else if (line_info.shape == 3) { 879// intersect = off_intersect (&l0, &l3, NULL, NULL, x,y,z, vx,vy,vz, offdata); 880// intersecti = 0; 881// } 882// 883// if (!intersect) { 884// /* Strange error: did not hit cylinder */ 885// if (line_info.flag_warning < 100) 886// fprintf(stderr, "PowderN: %s: FATAL ERROR: Did not hit sample from inside (inc).\n", NAME_CURRENT_COMP); 887// line_info.flag_warning++; 888// ABSORB; 889// } 890// 891// if (!intersecti) { 892// l1 = (l3 + l0) /2; 893// l2 = l1; 894// } 895// 896// if (concentric && intersecti) { 897// /* In case of concentricity, 'remove' backward wall of sample */ 898// l2 = l1; 899// l3 = l1; 900// } 901// 902// if(l0 < 0) l0=0; /* already in sample */ 903// if(l1 < 0) l1=0; /* already in inner hollow */ 904// if(l2 < 0) l2=0; /* already past inner hollow */ 905// 906// 907// l_1 = v*(l3 - l2 + l1 - l0); /* Length to exit */ 908// 909// pmul = l_full*my_inc*exp(-(my_a_v/v+my_s)*(l+l_1))/(p_inc); 910// pmul *= solid_angle/(4*PI); 911// 912// } /* Incoherent scattering event */ 913 else if (type == 1) { 914 /* Make transmitted (absorption-corrected) event */ 915 /* No coordinate changes here, simply change xray weight */ 916 pmul = exp(-(my_abs+my_s)*(l))/(p_transmit); 917 } 918 p *= pmul; 919 } /* Neutron leaving since it has passed already */ 920 } /* else transmit non interacting xrays */ 921 922%} 923 924FINALLY 925%{ 926 free(line_info.list); 927 free(my_q); 928 free(my_w); 929 free(my_s_k2); 930 if (line_info.flag_warning) 931 fprintf(stderr, "PowderN: %s: Error message was repeated %i times with absorbed xrays.\n", 932 NAME_CURRENT_COMP, line_info.flag_warning); 933 934%} 935 936MCDISPLAY 937%{ 938 if (line_info.V_0) { 939 magnify("xyz"); 940 if (line_info.shape == 0) { /* cyl */ 941 circle("xz", 0, yheight/2.0, 0, radius); 942 circle("xz", 0, -yheight/2.0, 0, radius); 943 line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0); 944 line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0); 945 line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius); 946 line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius); 947 if (thickness) { 948 double radius_i=radius-thickness; 949 circle("xz", 0, yheight/2.0, 0, radius_i); 950 circle("xz", 0, -yheight/2.0, 0, radius_i); 951 line(-radius_i, -yheight/2.0, 0, -radius_i, +yheight/2.0, 0); 952 line(+radius_i, -yheight/2.0, 0, +radius_i, +yheight/2.0, 0); 953 line(0, -yheight/2.0, -radius_i, 0, +yheight/2.0, -radius_i); 954 line(0, -yheight/2.0, +radius_i, 0, +yheight/2.0, +radius_i); 955 } 956 } else if (line_info.shape == 1) { /* box */ 957 double xmin = -0.5*xwidth; 958 double xmax = 0.5*xwidth; 959 double ymin = -0.5*yheight; 960 double ymax = 0.5*yheight; 961 double zmin = -0.5*zdepth; 962 double zmax = 0.5*zdepth; 963 multiline(5, xmin, ymin, zmin, 964 xmax, ymin, zmin, 965 xmax, ymax, zmin, 966 xmin, ymax, zmin, 967 xmin, ymin, zmin); 968 multiline(5, xmin, ymin, zmax, 969 xmax, ymin, zmax, 970 xmax, ymax, zmax, 971 xmin, ymax, zmax, 972 xmin, ymin, zmax); 973 line(xmin, ymin, zmin, xmin, ymin, zmax); 974 line(xmax, ymin, zmin, xmax, ymin, zmax); 975 line(xmin, ymax, zmin, xmin, ymax, zmax); 976 line(xmax, ymax, zmin, xmax, ymax, zmax); 977 if (zdepth_i) { 978 xmin = -0.5*xwidth_i; 979 xmax = 0.5*xwidth_i; 980 ymin = -0.5*yheight_i; 981 ymax = 0.5*yheight_i; 982 zmin = -0.5*zdepth_i; 983 zmax = 0.5*zdepth_i; 984 multiline(5, xmin, ymin, zmin, 985 xmax, ymin, zmin, 986 xmax, ymax, zmin, 987 xmin, ymax, zmin, 988 xmin, ymin, zmin); 989 multiline(5, xmin, ymin, zmax, 990 xmax, ymin, zmax, 991 xmax, ymax, zmax, 992 xmin, ymax, zmax, 993 xmin, ymin, zmax); 994 line(xmin, ymin, zmin, xmin, ymin, zmax); 995 line(xmax, ymin, zmin, xmax, ymin, zmax); 996 line(xmin, ymax, zmin, xmin, ymax, zmax); 997 line(xmax, ymax, zmin, xmax, ymax, zmax); 998 } 999 } if (line_info.shape == 2) { /* sphere */ 1000 if (radius_i) { 1001 circle("xy",0,0,0,radius_i); 1002 circle("xz",0,0,0,radius_i); 1003 circle("yz",0,0,0,radius_i); 1004 } 1005 circle("xy",0,0,0,radius); 1006 circle("xz",0,0,0,radius); 1007 circle("yz",0,0,0,radius); 1008 } else if (line_info.shape == 3) { /* OFF file */ 1009 off_display(offdata); 1010 } 1011 } 1012%} 1013END 1014