1/***************************************************************************** 2* 3* McStas, neutron ray-tracing package 4* Copyright(C) 2007 Risoe National Laboratory. 5* 6* Component: Isotropic_Sqw 7* 8* %I 9* Written by: E. Farhi, V. Hugouvieux 10* Date: August 2003 11* Origin: ILL 12* Modified by: E. Farhi, Jul 2005: made it work, concentric mode, multiple use 13* Modified by: E. Farhi, Mar 2007: improved implementation, correct small bugs 14* Modified by: E. Farhi, Oct 2008: added any shape sample geometry 15* Modified by: E. Farhi, Oct 2012: improved sampling scheme, correct bug in powder S(q) 16* 17* Isotropic sample handling multiple scattering and absorption for a general 18* S(q,w) (coherent and/or incoherent/self) 19* 20* %D 21* An isotropic sample handling multiple scattering and including as input the 22* dynamic structure factor of the chosen sample (e.g. from Molecular 23* Dynamics). Handles elastic/inelastic, coherent and incoherent scattering - 24* depending on the input S(q,w) - with multiple scattering and absorption. 25* Only the norm of q is handled (not the vector), and thus suitable for 26* liquids, gazes, amorphous and powder samples. 27* 28* If incoherent/self S(q,w) file is specified as empty (0 or "") then the 29* scattering is constant isotropic (Vanadium like). 30* In case you only have one S(q,w) data containing both coherent and 31* incoherent contributions you should e.g. use 'Sqw_coh' and set 'sigma_coh' 32* to the total scattering cross section. 33* The implementation assumes that the S(q,w) data is normalized, i.e. S(q) 34* goes to 1 for large q. If this is not the case, the component can do that 35* when 'norm=-1'. Alternatively, the S(q,w) data will be multiplied by 36* 'norm' for positive values. Both non symmetric and classical S(q,w) 37* data sets are handled by mean of the 'classical' parameter (see below). 38* 39* Additionally, for single order scattering (order=1), you may restrict the 40* vertical spreading of the scattering area using d_phi parameter. 41* 42* An important option to enhance statistics is to set 'p_interact' to, say, 43* 30 percent (0.3) in order to force a fraction of the beam to scatter. This 44* will result on a larger number of scattered events, retaining intensity. 45* 46* If you use this component and produce valuable scientific results, please 47* cite authors with references bellow (in <a href="#links">Links</a>). 48* E. Farhi et al, J Comp Phys 228 (2009) 5251 49* 50* <b>Sample shape:</b> 51* Sample shape may be a cylinder, a sphere, a box or any other shape 52* box/plate: xwidth x yheight x zdepth (thickness=0) 53* hollow box/plate:xwidth x yheight x zdepth and thickness>0 54* cylinder: radius x yheight (thickness=0) 55* hollow cylinder: radius x yheight and thickness>0 56* sphere: radius (yheight=0 thickness=0) 57* hollow sphere: radius and thickness>0 (yheight=0) 58* any shape: geometry=OFF file 59* 60* The complex geometry option handles any closed non-convex polyhedra. 61* It computes the intersection points of the neutron ray with the object 62* transparently, so that it can be used like a regular sample object. 63* It supports the OFF, PLY and NOFF file format but not COFF (colored faces). 64* Such files may be generated from XYZ data using: 65* qhull < coordinates.xyz Qx Qv Tv o > geomview.off 66* or 67* powercrust coordinates.xyz 68* and viewed with geomview or java -jar jroff.jar (see below). 69* The default size of the object depends of the OFF file data, but its 70* bounding box may be resized using xwidth,yheight and zdepth. 71* 72* <b>Concentric components:</b> 73* This component has the ability to contain other components when used in 74* hollow cylinder geometry (namely sample environment, e.g. cryostat and 75* furnace structure). Such component 'shells' should be split into input and 76* output side surrounding the 'inside' components. First part must then use 77* 'concentric=1' flag to enter the inside part. The component itself must be 78* repeated to mark the end of the concentric zone. The number of concentric 79* shells and number of components inside is not limited. 80* 81* COMPONENT S_in = Isotropic_Sqw(Sqw_coh="Al.laz", concentric=1, ...) 82* AT (0,0,0) RELATIVE sample_position 83* 84* COMPONENT something_inside ... // e.g. the sample itself or other materials 85* 86* COMPONENT S_out = COPY(S_in)(concentric=0) 87* AT (0,0,0) RELATIVE sample_position 88* 89* <b>Sqw file format:</b> 90* File format for S(Q,w) (coherent and incoherent) should contain 3 numerical 91* blocks, defining q axis values (vector), then energy axis values (vector), 92* then a matrix with one line per q axis value, containing Sqw values for 93* each energy axis value. Comments (starting with '#') and non numerical lines 94* are ignored and used to separate blocks. Sampling must be regular. 95* Some parameters can be specified in comment lines, namely (00 is a numerical value): 96* 97* # sigma_abs 00 absorption scattering cross section in [barn] 98* # sigma_inc 00 coherent scattering cross section in [barn] 99* # sigma_coh 00 incoherent scattering cross section in [barn] 100* # Temperature 00 in [K] 101* # V_rho 00 atom density per Angs^3 102* # density 00 in [g/cm^3] 103* # weight 00 in [g/mol] 104* # classical 00 [0=contains Bose factor (measurement) ; 1=classical symmetric] 105* 106* Example: 107* # q axis values 108* # vector of m values in Angstroem-1 109* 0.001000 .... 3.591000 110* # w axis values 111* # vector of n values in meV 112* 0.001391 ... 1.681391 113* # sqw values (one line per q axis value) 114* # matrix of S(q,w) values (m rows x n values), one line per q value, 115* 9.721422 10.599145 ... 0.000000 116* 10.054191 11.025244 ... 0.000000 117* ... 118* 0.000000 ... 3.860253 119* 120* See for instance file He4_liq_coh.sqw. Such files may be obtained from e.g. INX, 121* Nathan, Lamp and IDA softwares, as well as Molecular Dynamics (nMoldyn). 122* When the provided S(q,w) data is obtained from the classical correlation function 123* G(r,t), which is real and symmetric in time, the 'classical=1' parameter 124* should be set in order to multiply the file data with exp(hw/2kT). Otherwise, 125* the S(q,w) is NOT symmetrised (classical). If the S(q,w) data set includes both 126* negative and positive energy values, setting 'classical=-1' will attempt to 127* guess what type of S(q,w) it is. The temperature can also be determined this way. 128* In case you do not know if the data is classical or quantum, assume it is classical 129* at high temperatures, and quantum otherwise (T < typical mode excitations). 130* The positive energy values correspond to Stokes processes, i.e. material gains 131* energy, and neutrons loose energy. The energy range is symmetrized to allow up 132* and down scattering, taking into account detailed balance exp(-hw/2kT). 133* 134* You may also generate such S(q,w) 2D files using <a href="http://ifit.mccode.org/McStas.html#mozTocId297488">iFit </a> 135* 136* <b>Powder file format:</b> 137* Files for coherent elastic powder scattering may also be used. 138* Format specification follows the same principle as in the PowderN 139* component, with parameters: 140* 141* powder_format=Crystallographica 142* or powder_format=Fullprof 143* or powder_format=Lazy 144* or powder_format={j,d,F2,DW,Delta_d/d,1/2d,q,F,strain} (column indexes 1:n) 145* 146* or column indexes (starting from 1) given as comments in the file header 147* (e.g. '#column_j 4'). Refer to the PowderN component for more details. 148* Delta_d/d and Debye-Waller factor may be specified for all lines with the 149* 'powder_Dd' and 'powder_DW' parameters. 150* The reflection list should be ordered by decreasing d-spacing values. 151* 152* Additionally a special [q,Sq] format is also defined with: 153* powder_format=qSq 154* for which column 1 is 'q' and column 2 is 'S(q)'. 155* 156* <b>Examples:</b> 157* 1- Vanadium-like incoherent elastic scattering 158* Isotropic_Sqw(radius=0.005, yheight=0.01, V_rho=1/13.827, 159* sigma_abs=5.08, sigma_inc=4.935, sigma_coh=0) 160* 161* 2- liq-4He parameters 162* Isotropic_Sqw(..., Sqw_coh="He4_liq_coh.sqw", T=10, p_interact=0.3) 163* 164* 3- powder sample 165* Isotropic_Sqw(..., Sqw_coh="Al.laz") 166* 167* %BUGS: 168* When used in concentric mode, multiple bouncing scattering 169* (traversing the hollow part) is not taken into account. 170* 171* %VALIDATION 172* For Vanadium incoherent scattering mode, V_sample, PowderN, Single_crystal 173* and Isotropic_Sqw produce equivalent results, eventhough the two later are 174* more accurate (geometry, multiple scattering). Isotropic_Sqw gives same 175* powder patterns as PowderN, with an intensity within 20 %. 176* 177* %P 178* INPUT PARAMETERS: 179* Sqw_coh: [str] Name of the file containing the values of Q, w and S(Q,w) Coherent part; Q in Angs-1, E in meV, S(q,w) in meV-1. Use 0, NULL or "" to disable. 180* Sqw_inc: [str] Name of the file containing the values of Q, w and S(Q,w). Incoherent (self) part. Use 0, NULL or "" to scatter isotropically (V-like). 181* sigma_coh: [barns] Coherent Scattering cross-section. Use -1 to unactivate. 182* sigma_inc: [barns] Incoherent Scattering cross-section. Use -1 to unactivate. 183* sigma_abs: [barns] Absorption cross-section at 2200 m/s. Use -1 to unactivate. 184* rho: [AA-3] Density of scattering elements (nb atoms/unit cell V_0). 185* T: [K] Temperature of sample, detailed balance. Use T=0 to disable it. and T=-1 to automatically set it from non-classical S(q,w). 186* 187* Geometry parameters: 188* radius: [m] Outer radius of sample in (x,z) plane. cylinder/sphere. 189* xwidth: [m] width for a box sample shape 190* yheight: [m] Height of sample in vertical direction for box/cylinder shapes 191* zdepth: [m] depth for a box sample shape 192* thickness: [m] Thickness of hollow sample Negative value extends the hollow volume outside of the box/cylinder. 193* 194* OPTIONAL PARAMETERS: 195* 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) [1] 196* 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 197* order: [1] Limit multiple scattering up to given order 0:all (default), 1:single, 2:double, ... 198* verbose: [1] Verbosity level (0:silent, 1:normal, 2:verbose, 3:debug). A verbosity>1 also computes dispersions and S(q,w) analysis. 199* d_phi: [deg] scattering vertical angular spreading (usually the height of the next component/detector). Use 0 for full space. This is only relevant for single scattering (order=1). 200* weight: [g/mol] atomic/molecular weight of material 201* density: [g/cm^3] density of material. V_rho=density/weight/1e24*N_A 202* threshold: [1] Value under which S(Q,w) is not accounted for. to set according to the S(Q,w) values, i.e. not too low. 203* p_interact: [1] Force a given fraction of the beam to scatter, keeping intensity right, to enhance small signals (-1 unactivate). 204* norm: [1] Normalize S(q,w) when -1 (default). Use raw data when 0, multiply S(q,w) when norm>0. 205* classical: [1] Assumes the S(q,w) data from the files is a classical S(q,w), and multiply that data by exp(hw/2kT) on up/down energy sides. Use 0 when obtained from raw experiments, 1 from molecular dynamics. Use -1 to guess from a data set including both energy sides. 206* quantum_correction: [str] Specify the type of quantum correction to use "Boltzmann"/"Schofield","harmonic"/"Bader" or "standard"/"Frommhold" (default) 207* 208* POWDER ELASTIC SCATTERING PARAMETERS (see PowderN for more details): 209* powder_Dd: [1] global Delta_d/d spreading, or 0 if ideal. 210* powder_DW: [1] global Debey-Waller factor, if not in |F2| or 1. 211* powder_format: [no quotes] name or definition of column indexes in file 212* powder_Vc: [AA^3] volume of the unit cell 213* powder_barns: [1] 0 when |F2| data in powder file are fm^2, 1 when in barns (barns=1 for laz, barns=0 for lau type files). 214* 215* OUTPUT PARAMETERS: 216* VarSqw: [] internal structure including dq=wavevector transfer Kf-Ki in Angs-1; dw=energy transfer Ef-Ei in meV; type=interaction type of event as 'c' (coherent), 't' (transmitted) 'i' (incoherent) 'v' (isotropic incoherent, Vanadium-like). 217* SCATTERED: [] order of multiple scattering 218* 219* %Links 220* E. Farhi, V. Hugouvieux, M.R. Johnson, and W. Kob, Journal of Computational Physics 228 (2009) 5251-5261 "Virtual experiments: Combining realistic neutron scattering instrument and sample simulations" 221* %L 222* Hugouvieux V, Farhi E, Johnson MR, Physica B, 350 (2004) 151 "Virtual neutron scattering experiments" 223* %L 224* Hugouvieux V, PhD, University of Montpellier II, France (2004). 225* %L 226* H. Schober, Collection SFN 10 (2010) 159-336 227* %L 228* H.E. Fischer, A.C. Barnes, and P.S. Salmon. Rep. Prog. Phys., 69 (2006) 233 229* %L 230* P.A. Egelstaff, An Introduction to the Liquid State, 2nd Ed., Oxford Science Pub., Clarendon Press (1992). 231* %L 232* S. W. Lovesey, Theory of Neutron Scattering from Condensed Matter, Vol1, Oxford Science Pub., Clarendon Press (1984). 233* %L 234* <a href="http://www.ncnr.nist.gov/resources/n-lengths/">Cross sections for single elements</a> 235* %L 236* <a href="http://www.ncnr.nist.gov/resources/sldcalc.html>Cross sections for compounds</a> 237* %L 238* <a href="http://www.webelements.com/">Web Elements</a> 239* %L 240* <a href="http://www.ill.eu/sites/fullprof/index.html">Fullprof</a> powder refinement 241* %L 242* <a href="http://www.crystallographica.com/">Crystallographica</a> software 243* %L 244* Example data file <a href="../data/He4_liq_coh.sqw">He4_liq_coh.sqw</a> 245* %L 246* The <a href="PowderN.html">PowderN</a> component. 247* %L 248* The test/example instrument <a href="../examples/Test_Isotropic_Sqw.instr">Test_Isotropic_Sqw.instr</a>. 249* %L 250* <a href="http://www.geomview.org">Geomview and Object File Format (OFF)</a> 251* %L 252* Java version of Geomview (display only) <a href="http://www.holmes3d.net/graphics/roffview/">jroff.jar</a> 253* %L 254* <a href="http://qhull.org">qhull</a> 255* %L 256* <a href="http://www.cs.ucdavis.edu/~amenta/powercrust.html">powercrust</a> 257* %E 258***********************************************************/ 259 260DEFINE COMPONENT Isotropic_Sqw 261DEFINITION PARAMETERS (powder_format=Undefined) 262SETTING PARAMETERS(string Sqw_coh=0, string Sqw_inc=0, string geometry=0, 263 radius=0,thickness=0, 264 xwidth=0, yheight=0, zdepth=0, 265 threshold=1e-20, int order=0, T=0, verbose=1, d_phi=0, int concentric=0, 266 rho=0, sigma_abs=0, sigma_coh=0, sigma_inc=0, classical=-1, 267 powder_Dd=0, powder_DW=0, powder_Vc=0, density=0, weight=0, 268 p_interact=-1, norm=-1, powder_barns=1, string quantum_correction="Frommhold") 269OUTPUT PARAMETERS (VarSqw, columns, offdata) 270/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ 271/*****************************************************************************/ 272/*****************************************************************************/ 273 274/* SHARE functions: 275* void Sqw_Data_init (struct Sqw_Data_struct *Sqw_Data) 276* t_Table *Sqw_read_PowderN(struct Sqw_sample_struct *Sqw, t_Table sqwTable) 277* int Sqw_search_SW(struct Sqw_Data_struct Sqw, double randnum) 278* int Sqw_search_Q_proba_per_w(struct Sqw_Data_struct Sqw, double randnum, int index) 279* double Sqw_init(struct Sqw_sample_struct *Sqw, char *file_coh, char *file_inc) 280* double Sqw_integrate_iqSq(struct Sqw_Data_struct *Sqw_Data, double Ei) 281* void Sqw_diagnosis(struct Sqw_sample_struct *Sqw, struct Sqw_Data_struct *Sqw_Data) 282* struct Sqw_Data_struct *Sqw_readfile( 283struct Sqw_sample_struct *Sqw, char *file, struct Sqw_Data_struct *Sqw_Data) 284*/ 285SHARE 286%{ 287 288#ifndef ISOTROPIC_SQW 289#define ISOTROPIC_SQW $Revision$ 290 291/* {j d F2 DW Dd inv2d q F} + { Sq if j == -1}*/ 292#ifndef Crystallographica 293#define Crystallographica { 4,5,7,0,0,0,0, 0,0 } 294#define Fullprof { 4,0,8,0,0,5,0, 0,0 } 295#define Undefined { 0,0,0,0,0,0,0, 0,0 } 296#define Lazy {17,6,0,0,0,0,0,13,0 } 297#endif 298/* special case for [q,Sq] table */ 299#define qSq {-1,0,0,0,0,0,1, 0,0 } 300 301%include "read_table-lib" 302%include "interoff-lib" 303 304/* For the density of states S(w) */ 305struct Sqw_W_struct 306{ 307 double omega; /* omega value for the data block */ 308 double cumul_proba; /* cumulated intensity (between 0 and 1) */ 309}; 310 311/* For the S(q|w) probabilities */ 312struct Sqw_Q_struct 313{ 314 double Q; /* omega value for the data block */ 315 double cumul_proba; /* normalized cumulated probability */ 316}; 317 318struct Sqw_Data_struct /* contains normalized Sqw data for probabilities, coh and inc */ 319{ 320 struct Sqw_W_struct *SW; /* P(w) ~ density of states */ 321 struct Sqw_Q_struct **SQW; /* P(Q|w)= probability of each Q with w */ 322 323 long *SW_lookup; 324 long **QW_lookup; 325 t_Table Sqw; /* S(q,w) rebin from file in range -w_max:w_max and 0:q_max, with exp(-hw/kT) weight */ 326 t_Table iqSq; /* sigma(Ei) = sigma/2/Ki^2 * \int q S(q,w) dq dw up to 2*Ki_max */ 327 long q_bins; 328 long w_bins; /* length of q and w vectors/axes from file */ 329 double q_max, q_step; /* min=0 */ 330 double w_max, w_step; /* min=-w_max */ 331 long lookup_length; 332 char filename[80]; 333 double intensity; 334 double Ei_max; /* max neutron incoming energy for Sigma=iqSq table */ 335 long iqSq_length; 336 char type; 337 double q_min_file; 338}; 339 340struct Sqw_sample_struct { /* global parameters gathered as a structure */ 341 char compname[256]; 342 343 struct Sqw_Data_struct Data_inc; 344 struct Sqw_Data_struct Data_coh; 345 346 double s_abs, s_coh, s_inc; /* material constants */ 347 double my_s; 348 double my_a_v; 349 double mat_rho; 350 double mat_weight; 351 double mat_density; 352 double Temperature; /* temperature from the data file */ 353 int shape; /* 0:cylinder, 1:box, 2:sphere 3:any shape*/ 354 355 double sqw_threshold; /* options to tune S(q,w) */ 356 double sqw_classical; 357 double sqw_norm; 358 359 double barns; /* for powders */ 360 double Dd, DWfactor; 361 362 double T2E; /* constants */ 363 char Q_correction[256]; 364 double sqSE2K; 365 366 int maxloop; /* flags to monitor caught warnings */ 367 int minevents; 368 long neutron_removed; 369 long neutron_enter; 370 long neutron_pmult; 371 long neutron_exit; 372 char verbose_output; 373 int column_order[9]; /* column signification */ 374 long lookup_length; 375 376 double dq, dw; /* q/w transfer */ 377 char type; /* interaction type: c(coherent), i(incoherent), 378 V(isotropic incoherent), t(transmitted) */ 379 /* store information from the last event */ 380 double ki_x,ki_y,ki_z,kf_x,kf_y,kf_z; 381 double ti, tf; 382 double vi, vf; 383 double ki, kf; 384 double theta; 385 386 double mean_scatt; /* stat to show at the end */ 387 double mean_abs; 388 double psum_scatt; 389 double single_coh; 390 double single_inc; 391 double multi; 392 393 double rw, rq; 394}; 395 396#include <stdio.h> 397#include <math.h> 398 399/* sets a Data S(q,w) to 'NULL' */ 400void Sqw_Data_init(struct Sqw_Data_struct *Sqw_Data) 401{ 402 Sqw_Data->q_bins =0; 403 Sqw_Data->w_bins =0; 404 Sqw_Data->q_max =0; 405 Sqw_Data->q_step =1; 406 Sqw_Data->w_max =0; 407 Sqw_Data->w_step =1; 408 Sqw_Data->Ei_max = 0; 409 Sqw_Data->lookup_length=100; /* length of lookup tables */ 410 Sqw_Data->intensity =0; 411 strcpy(Sqw_Data->filename, ""); 412 Sqw_Data->SW =NULL; 413 Sqw_Data->SQW =NULL; 414 Sqw_Data->SW_lookup =NULL; 415 Sqw_Data->QW_lookup =NULL; 416 Sqw_Data->iqSq_length =100; 417 Sqw_Data->type = ' '; 418 Sqw_Data->q_min_file = 0; 419} 420 421off_struct offdata; 422 423/* gaussian distribution to appply around Bragg peaks in a powder */ 424double Sqw_powder_gauss(double x, double mean, double rms) { 425 return (exp(-(x-mean)*(x-mean)/(2*rms*rms))/(sqrt(2*PI)*rms)); 426} 427 428/* Sqw_quantum_correction 429* 430* Return the 'quantum correction factor Q so that: 431* 432* S(q, w) = Q(w) S*(q,w) 433* S(q,-w) = exp(-hw/kT) S(q,w) 434* S(q, w) = exp( hw/kT) S(q,-w) 435* 436* with S*=classical limit and Q(w) defined below. For omega > 0, S(q,w) > S(q,-w) 437* 438* input: 439* w: energy [meV] 440* T: temperature [K] 441* type: 'Schofield' or 'Boltzmann' Q = exp(hw/kT/2) 442* 'harmonic' or 'Bader' Q = hw/kT./(1-exp(-hw/kT)) 443* 'standard' or 'Frommhold' Q = 2./(1+exp(-hw/kT)) [recommended] 444* 445* References: 446* B. Hehr, http://www.lib.ncsu.edu/resolver/1840.16/7422 PhD manuscript (2010). 447* S. A. Egorov, K. F. Everitt and J. L. Skinner. J. Phys. Chem., 103, 9494 (1999). 448* P. Schofield. Phys. Rev. Lett., 4, 239 (1960). 449* J. S. Bader and B. J. Berne. J. Chem. Phys., 100, 8359 (1994). 450* T. D. Hone and G. A. Voth. J. Chem. Phys., 121, 6412 (2004). 451* L. Frommhold. Collision-induced absorption in gases, 1 st ed., Cambridge 452* Monographs on Atomic, Molecular, and Chemical Physics, Vol. 2, 453* Cambridge Univ. Press: London (1993). 454 455 */ 456double Sqw_quantum_correction(double hw, double T, char *type) { 457 double Q = 1; 458 double kT = T/11.605; /* [K] -> [meV = 1000*KB/e] */ 459 if (!hw || !T) return 1; 460 if (type == NULL || !strcmp(type, "standard") 461 || !strcmp(type, "Frommhold") || !strcmp(type, "default")) 462 Q = 2/(1+exp(-hw/kT)); 463 if (!strcmp(type, "Schofield") || !strcmp(type, "Boltzmann")) 464 Q = exp(hw/kT/2); 465 if (!strcmp(type, "harmonic") || !strcmp(type, "Bader")) 466 Q = hw/kT/(1-exp(-hw/kT)); 467 468 return Q; 469} 470 471/***************************************************************************** 472* Sqw_read_PowderN: Read PowderN data files 473* Returns t_Table array or NULL in case of error 474* Used in : Sqw_readfile (1) 475*****************************************************************************/ 476t_Table *Sqw_read_PowderN(struct Sqw_sample_struct *Sqw, t_Table sqwTable) 477{ 478 struct line_data 479 { 480 double F2; /* Value of structure factor */ 481 double q; /* Q vector */ 482 int j; /* Multiplicity */ 483 double DWfactor; /* Debye-Waller factor */ 484 double w; /* Intrinsic line width */ 485 }; 486 struct line_data *list = NULL; 487 double q_count=0, j_count=0, F2_count=0; 488 int mult_count =0; 489 double q_step =FLT_MAX; 490 long size =0; 491 int i, index; 492 double q_min=0, q_max=0; 493 char flag=0; 494 int list_count=0; 495 double q_step_cur; 496 char flag_qSq = 0; 497 498 t_Table *retTable; 499 500 flag_qSq = (Sqw->column_order[8]>0 && Sqw->column_order[6]>0); 501 502 if (Sqw->column_order[0] == 4 && Sqw->barns !=0) 503 printf("Isotropic_sqw: %s: Powder file probably of type Crystallographica/Fullprof (lau)\n" 504 "WARNING: but F2 unit is set to powder_barns=1 (barns). Intensity might be 100 times too high.\n", 505 Sqw->compname); 506 if (Sqw->column_order[0] == 17 && Sqw->barns == 0) 507 printf("Isotropic_sqw: %s: Powder file probably of type Lazy Pulver (laz)\n" 508 "WARNING: but F2 unit is set to powder_barns=0 (fm^2). Intensity might be 100 times too low.\n", 509 Sqw->compname); 510 511 size = sqwTable.rows; 512 MPI_MASTER( 513 if (Sqw->verbose_output > 0) { 514 printf("Isotropic_sqw: Converting %ld %s from %s into S(q,w) data\n", 515 size, flag_qSq ? "S(q)" : "powder lines", sqwTable.filename); 516 } 517 ); 518 /* allocate line_data array */ 519 list = (struct line_data*)malloc(size*sizeof(struct line_data)); 520 521 for (i=0; i<size; i++) 522 { 523 /* printf("Reading in line %i\n",i);*/ 524 double j=0, d=0, w=0, DWfactor=0, F2=0, Sq=-1, q=0; 525 int index; 526 527 if (Sqw->Dd >= 0) w = Sqw->Dd; 528 if (Sqw->DWfactor > 0) DWfactor = Sqw->DWfactor; 529 530 /* get data from table using columns {j d F2 DW Dd inv2d q} + { Sq }*/ 531 /* column indexes start at 1, thus need to substract 1 */ 532 if (Sqw->column_order[0]>0) 533 j = Table_Index(sqwTable, i, Sqw->column_order[0]-1); 534 if (Sqw->column_order[1]>0) 535 d = Table_Index(sqwTable, i, Sqw->column_order[1]-1); 536 if (Sqw->column_order[2]>0) 537 F2 = Table_Index(sqwTable, i, Sqw->column_order[2]-1); 538 if (Sqw->column_order[3]>0) 539 DWfactor = Table_Index(sqwTable, i, Sqw->column_order[3]-1); 540 if (Sqw->column_order[4]>0) 541 w = Table_Index(sqwTable, i, Sqw->column_order[4]-1); 542 if (Sqw->column_order[5]>0) { 543 d = Table_Index(sqwTable, i, Sqw->column_order[5]-1); if (d) d = 1/d/2; } 544 if (Sqw->column_order[6]>0) 545 q = Table_Index(sqwTable, i, Sqw->column_order[6]-1); 546 if (Sqw->column_order[7]>0 && !F2) 547 {F2= Table_Index(sqwTable, i, Sqw->column_order[7]-1); F2 *= F2;} 548 549 if (Sqw->column_order[8]>0) 550 Sq= Table_Index(sqwTable, i, Sqw->column_order[8]-1); 551 552 if (q > 0 && Sq >= 0) F2 = Sq; 553 if (d > 0 && q <= 0) q = 2*PI/d; 554 555 /* assign and check values */ 556 j = (j > 0 ? j : 0); 557 if (flag_qSq) j=1; 558 DWfactor = (DWfactor > 0 ? DWfactor : 1); 559 w = (w>0 ? w : 0); 560 F2 = (F2 >= 0 ? F2 : 0); 561 d = (q > 0 ? 2*PI/d : 0); 562 if (j == 0 || d == 0 || q == 0) { 563 printf("Isotropic_sqw: %s: Warning: line %i has invalid definition\n" 564 " (mult=0 or q=0 or d=0)\n", Sqw->compname, i); 565 continue; 566 } 567 list[list_count].j = j; 568 list[list_count].q = q; 569 list[list_count].DWfactor = DWfactor; 570 list[list_count].w = w; 571 list[list_count].F2= F2; /* or S(q) if flag_qSq */ 572 573 if (q_max < d) q_max = q; 574 if (q_min > d) q_min = q; 575 if (list_count > 1) { 576 q_step_cur = fabs(list[list_count].q - list[list_count-1].q); 577 if (q_step_cur > 1e-5 && (!q_step || q_step_cur < q_step)) 578 q_step = q_step_cur; 579 } 580 581 /* adjust multiplicity if j-column + multiple d-spacing lines */ 582 /* if d = previous d, increase line duplication index */ 583 if (!q_count) q_count = q; 584 if (!j_count) j_count = j; 585 if (!F2_count) F2_count= F2; 586 if (fabs(q_count-q) < 0.0001*fabs(q) 587 && fabs(F2_count-F2) < 0.0001*fabs(F2) && j_count == j) { 588 mult_count++; flag=0; } 589 else flag=1; 590 if (i == size-1) flag=1; 591 /* else if d != previous d : just passed equivalent lines */ 592 if (flag) { 593 if (i == size-1) list_count++; 594 /* if duplication index == previous multiplicity */ 595 /* set back multiplicity of previous lines to 1 */ 596 if (Sqw->verbose_output > 2 && (mult_count == list[list_count-1].j 597 || (mult_count == list[list_count].j && i == size-1))) { 598 printf("Isotropic_Sqw: %s: Setting multiplicity to 1 for lines [%i:%i]\n" 599 " (d-spacing %g is duplicated %i times)\n", 600 Sqw->compname, list_count-mult_count, list_count-1, list[list_count-1].q, mult_count); 601 for (index=list_count-mult_count; index<list_count; list[index++].j = 1); 602 mult_count = 1; 603 q_count = q; 604 j_count = j; 605 F2_count= F2; 606 } 607 if (i == size-1) list_count--; 608 flag=0; 609 } 610 list_count++; 611 } /* end for */ 612 613 /* now builds new Table_Array to continue with Sqw_readfile */ 614 if (q_max == q_min || !q_step) return(NULL); 615 if (!flag_qSq) 616 size = 3*q_max/q_step; /* set a default of 3 q values per line */ 617 else size = list_count; 618 /* update the value of q_step */ 619 q_step = q_max/size; 620 MPI_MASTER( 621 if (Sqw->verbose_output > 0) 622 printf("Isotropic_sqw: q range [%g:%g], creating %li elements vector\n", 623 q_min, q_max, size); 624 ); 625 626 retTable = (t_Table*)calloc(4, sizeof(t_Table)); 627 if (!retTable) printf("Isotropic_Sqw: ERROR: Cannot allocate PowderN->Sqw table.\n"); 628 else { 629 char *header; 630 if (!Table_Init(&retTable[0], size, 1)) 631 { printf("Isotropic_Sqw: ERROR Cannot allocate q-axis [%li] from Powder lines.\n", size); return(NULL); } 632 if (!Table_Init(&retTable[1], 1, 1)) 633 { printf("Isotropic_Sqw: ERROR Cannot allocate w-axis from Powder lines.\n"); return(NULL); } 634 if (!Table_Init(&retTable[2], size, 1)) 635 { printf("Isotropic_Sqw: ERROR Cannot allocate Sqw [%li] from Powder lines.\n", size); return(NULL); } 636 Table_Init(&retTable[3], 0,0); 637 638 header = malloc(64); if (header) 639 { retTable[0].header = header; strcpy(retTable[0].header, "q"); } 640 header = malloc(64); if (header) 641 { retTable[1].header = header; strcpy(retTable[1].header, "w"); } 642 header = malloc(64); if (header) 643 { retTable[2].header = header; strcpy(retTable[2].header, "Sqw"); } 644 for (i=0; i < 4; i++) { 645 retTable[i].array_length = 3; 646 retTable[i].block_number = i+1; 647 } 648 if (!flag_qSq) 649 for (i=0; i<size; i++) 650 retTable[0].data[i] = i*q_max/size; 651 for (i=0; i<list_count; i++) { /* loop on each Bragg peak */ 652 double peak_qmin, peak_qmax,factor,q; 653 if (list[i].w > 0 && !flag_qSq) { 654 peak_qmin = list[i].q*(1 - list[i].w*3); 655 peak_qmax = list[i].q*(1 + list[i].w*3); 656 } else { /* Dirac peak, no width */ 657 peak_qmin = peak_qmax = list[i].q; 658 } 659 /* S(q) intensity is here */ 660 factor = list[i].j*(list[i].DWfactor ? list[i].DWfactor : 1) 661 *Sqw->mat_rho*PI/2 662 /(Sqw->type == 'c' ? Sqw->s_coh : Sqw->s_inc)*list[i].F2/list[i].q/list[i].q; 663 if (Sqw->barns) factor *= 100; 664 for (q=peak_qmin; q <= peak_qmax; q += q_step) { 665 index = (long)floor(size*q/q_max); 666 if (index < 0) index=0; 667 else if (index >= size) index = size-1; 668 if (flag_qSq) { 669 retTable[2].data[index] += list[i].F2; 670 retTable[0].data[index] = list[i].q; 671 } else { 672 if (list[i].w <=0 || list[i].w*q < q_step) /* step function */ 673 retTable[2].data[index] += factor/q_step; 674 else /* gaussian */ 675 retTable[2].data[index] += factor 676 * Sqw_powder_gauss(q, list[i].q, list[i].w*list[i].q); 677 } 678 } 679 } /* end for i */ 680 Table_Stat(&retTable[0]); Table_Stat(&retTable[1]); Table_Stat(&retTable[2]); 681 Sqw->sqw_norm = 0; /* F2 are normalized already */ 682 } 683 684 return(retTable); 685} /* Sqw_read_PowderN */ 686 687/***************************************************************************** 688* Sqw_search_SW: For a given random number 'randnum', search for the bin 689* containing the corresponding Sqw->SW 690* Choose an energy in the projected S(w) distribution 691* Used in : TRACE (1) 692*****************************************************************************/ 693int Sqw_search_SW(struct Sqw_Data_struct Sqw, double randnum) 694{ 695 int index_w=0; 696 697 if (randnum <0) randnum=0; 698 if (randnum >1) randnum=1; 699 700 if (Sqw.w_bins == 1) return(0); 701 /* benefit from fast lookup table if exists */ 702 if (Sqw.SW_lookup) { 703 index_w = Sqw.SW_lookup[(long)floor(randnum*Sqw.lookup_length)]-1; 704 if (index_w<0) index_w=0; 705 } 706 707 while (index_w < Sqw.w_bins && (&(Sqw.SW[index_w]) != NULL) && (randnum > Sqw.SW[index_w].cumul_proba)) 708 index_w++; 709 if (index_w >= Sqw.w_bins) index_w = Sqw.w_bins-1; 710 711 if (&(Sqw.SW[index_w]) == NULL) 712 { 713 printf("Isotropic_Sqw: Warning: No corresponding value in the SW. randnum too big.\n"); 714 printf(" index_w=%i ; randnum=%f ; Sqw.SW[index_w-1].cumul_proba=%f (Sqw_search_SW)\n", 715 index_w, randnum, Sqw.SW[index_w-1].cumul_proba); 716 return index_w-1; 717 } 718 else 719 return (index_w); 720} 721 722/***************************************************************************** 723* Sqw_search_Q_proba_per_w: For a given random number randnum, search for 724* the bin containing the corresponding Sqw.SW in the Q probablility grid 725* Choose a momentum in the S(q|w) distribution 726* index is given by Sqw_search_SW 727* Used in : TRACE (1) 728*****************************************************************************/ 729int Sqw_search_Q_proba_per_w(struct Sqw_Data_struct Sqw, double randnum, int index_w) 730{ 731 int index_q=0; 732 733 if (randnum <0) randnum=0; 734 if (randnum >1) randnum=1; 735 736 /* benefit from fast lookup table if exists */ 737 if (Sqw.QW_lookup && Sqw.QW_lookup[index_w]) { 738 index_q = Sqw.QW_lookup[index_w][(long)floor(randnum*Sqw.lookup_length)]-1; 739 if (index_q<0) index_q=0; 740 } 741 742 while (index_q < Sqw.q_bins && (&(Sqw.SQW[index_w][index_q]) != NULL) 743 && (randnum > Sqw.SQW[index_w][index_q].cumul_proba)) { 744 index_q++; 745 } 746 if (index_q >= Sqw.q_bins) index_q = Sqw.q_bins-1; 747 748 if (&(Sqw.SQW[index_w][index_q]) == NULL) 749 return -1; 750 else 751 return (index_q); 752} 753 754/***************************************************************************** 755* compute the effective total cross section \int q S(q,w) dw dq 756* for incoming neutron energy 0 < Ei < 2*w_max, and 757* integration range w=-w_max:Ei and q=Q0:Q1 with 758* Q0 = SE2Q*(sqrt(Ei)-sqrt(Ei-w))=|Ki-Kf| 759* Q1 = SE2Q*(sqrt(Ei)+sqrt(Ei-w))=|Ki+Kf| 760* The data to use is Sqw_Data->Sqw, and the limits are Sqw_Data->w_max Sqw_Data->q_max 761* Returns the integral value 762* Used in: Sqw_readfile (1) 763*****************************************************************************/ 764double Sqw_integrate_iqSq(struct Sqw_Data_struct *Sqw_Data, double Ei) 765{ 766 long index_w; 767 double iqSq = 0; 768 /* w=Ei-Ef q=ki-kf w>0 neutron looses energy, Stokes, Ef = Ei-w > 0, Kf =|Ki-q| > 0 */ 769 for (index_w=0; index_w < Sqw_Data->w_bins; index_w++) { 770 long index_q; 771 double w = -Sqw_Data->w_max + index_w * Sqw_Data->w_step; /* in the Sqw table */ 772 if (w <= Ei) { /* integration range w=-w_max:Ei, Ef = Ei-w > 0 */ 773 double sq=0, Q0=0, Q1=0; 774 sq = sqrt(Ei-w); /* always real as test was true before */ 775 Q0 = SE2V*V2K*(sqrt(Ei)-sq); 776 Q1 = SE2V*V2K*(sqrt(Ei)+sq); 777 778 for (index_q=0; index_q < Sqw_Data->q_bins; index_q++) { 779 double q=(double)index_q * Sqw_Data->q_step; 780 /* add 'pixel' = q S(q,w) */ 781 if (Q0 <= q && q <= Q1) iqSq += q*Table_Index(Sqw_Data->Sqw, index_q, index_w); 782 } 783 } 784 } 785 /* multiply by 'pixel' size = dq dw */ 786 return(iqSq * Sqw_Data->q_step * Sqw_Data->w_step); 787} /* Sqw_integrate_iqSq */ 788 789 790/***************************************************************************** 791* Sqw_diagnosis: Computes Sqw_classical, moments and physical quantities 792* make consistency checks, and output some data files 793* Return: output files and information displayed 794* Used in: Sqw_init (2) 795*****************************************************************************/ 796void Sqw_diagnosis(struct Sqw_sample_struct *Sqw, struct Sqw_Data_struct *Sqw_Data) 797{ 798 799 t_Table Sqw_cl; /* the Sqw symmetric/classical version (T-> Inf) */ 800 t_Table Gqw; /* the generalized density of states as of Carpenter and Price, J Non Cryst Sol 92 (1987) 153 */ 801 t_Table Sqw_moments[7]; /* M0=S(q) M1=E_r M3 w_c w_l M0_cl=S_cl(q) G(w) */ 802 t_Table w_c, w_l; 803 long index_q, index_w; 804 char c[CHAR_BUF_LENGTH]; /* temporary variable */ 805 long q_min_index = 0; 806 807 char do_coh=0, do_inc=0; 808 double q_min =0; 809 double u2 =0, S0=1; 810 long u2_count=0; 811 812 if (!Sqw_Data || !Sqw_Data->intensity) return; /* nothing to do with empty S(q,w) */ 813 814 if (Sqw_Data->type=='c') do_coh = 1; 815 if (Sqw_Data->type=='i') do_inc = 1; 816 817 q_min = Sqw_Data->q_min_file; 818 if (q_min <= 0) q_min = Sqw_Data->q_step; 819 820 /* test if there is only one S(q,w) available */ 821 if (!((Sqw->Data_inc).intensity) || !((Sqw->Data_coh).intensity)) 822 do_coh = do_inc = 1; /* do both if only one file given */ 823 824 if (Sqw->Temperature > 0) { 825 if (!Table_Init(&Sqw_cl, Sqw_Data->q_bins, Sqw_Data->w_bins)) { 826 printf("Isotropic_Sqw: %s: Cannot allocate S_cl(q,w) Table (%lix%i).\n" 827 "WARNING Skipping S(q,w) diagnosis.\n", 828 Sqw->compname, Sqw_Data->q_bins, 1); 829 return; 830 } 831 sprintf(Sqw_cl.filename, 832 "S(q,w)_cl from %s (dynamic structure factor, classical)", 833 Sqw_Data->filename); 834 Sqw_cl.block_number = 1; 835 Sqw_cl.min_x = 0; 836 Sqw_cl.max_x = Sqw_Data->q_max; 837 Sqw_cl.step_x = Sqw_Data->q_step; 838 } 839 840 /* initialize moments and 1D stuff */ 841 for (index_q=0; index_q < 6; index_q++) { 842 if (!Table_Init(&Sqw_moments[index_q], Sqw_Data->q_bins, 1)) { 843 printf("Isotropic_Sqw: %s: Cannot allocate S(q,w) moment %ld Table (%lix%i).\n" 844 "WARNING Skipping S(q,w) diagnosis.\n", 845 Sqw->compname, index_q, Sqw_Data->q_bins, 1); 846 Table_Free(&Sqw_cl); 847 return; 848 } 849 Sqw_moments[index_q].block_number = 1; 850 Sqw_moments[index_q].min_x = 0; 851 Sqw_moments[index_q].max_x = Sqw_Data->q_max; 852 Sqw_moments[index_q].step_x = Sqw_Data->q_step; 853 } 854 index_q=6; 855 Table_Init(&Sqw_moments[index_q], Sqw_Data->w_bins, 1); 856 Sqw_moments[index_q].block_number = 1; 857 Sqw_moments[index_q].min_x = -Sqw_Data->w_max; 858 Sqw_moments[index_q].max_x = Sqw_Data->w_max; 859 Sqw_moments[index_q].step_x = Sqw_Data->w_step; 860 861 /* set Table titles */ 862 sprintf(Sqw_moments[0].filename, 863 "S(q)=M0(q) from %s [int S(q,w) dw]", 864 Sqw_Data->filename); 865 sprintf(Sqw_moments[1].filename, 866 "M1(q) 1-st moment from %s [int w S(q,w) dw] = HBAR^2*q^2/2/m (f-sum rule, recoil, Lovesey T1 Eq 3.63 p72, Egelstaff p196)", 867 Sqw_Data->filename); 868 sprintf(Sqw_moments[2].filename, 869 "M3(q) 3-rd moment from %s [int w^3 S(q,w) dw] = M1(q)*w_l^2(q)", 870 Sqw_Data->filename); 871 sprintf(Sqw_moments[3].filename, 872 "w_c(q) = sqrt(M1(q)/M0(q)*2kT) collective excitation from %s (Lovesey T1 Eq 5.38 p180, p211 Eq 5.204). Gaussian half-width of the S(q,w) classical", 873 Sqw_Data->filename); 874 sprintf(Sqw_moments[4].filename, 875 "w_l(q) = sqrt(M3(q)/M1(q)) harmonic frequency from %s (Lovesey T1 5.39 p 180)", 876 Sqw_Data->filename); 877 sprintf(Sqw_moments[5].filename, 878 "S_cl(q)=M0_cl(q) from %s [int S_cl(q,w) dw]", 879 Sqw_Data->filename); 880 sprintf(Sqw_moments[6].filename, 881 "G(w) generalized effective density of states from %s (Carpenter J Non Cryst Sol 92 (1987) 153)", 882 Sqw_Data->filename); 883 884 for (index_q=0; index_q < Sqw_Data->q_bins; index_q++) { 885 double q = index_q*Sqw_Data->q_step; /* q value in Sqw_full ; q_min = 0 */ 886 double sq = 0; /* S(q) = w0 = 0-th moment */ 887 double w1 = 0; /* first moment \int w Sqw dw */ 888 double w3 = 0; /* third moment \int w^3 Sqw dw */ 889 double sq_cl = 0; /* S(q) = M0 = 0-th moment classical */ 890 double w_c = 0; 891 double w_l = 0; 892 893 for (index_w=0; index_w < Sqw_Data->w_bins; index_w++) { 894 895 double w = -Sqw_Data->w_max + index_w*Sqw_Data->w_step; /* w value in Sqw_full */ 896 double sqw_cl =0; 897 double sqw_full =0; 898 899 sqw_full = Table_Index(Sqw_Data->Sqw, index_q, index_w); 900 901 /* Sqw moments */ 902 if (w && Sqw_Data->w_bins) { 903 double tmp; 904 tmp = sqw_full*Sqw_Data->w_step; 905 tmp *= w; w1 += tmp; 906 tmp *= w*w; w3 += tmp; 907 } 908 909 /* compute classical Sqw and S(q)_cl */ 910 if (Sqw->Temperature > 0) { 911 double n; 912 sqw_cl = sqw_full * Sqw_quantum_correction(-w,Sqw->Temperature,Sqw->Q_correction); 913 if (!Table_SetElement(&Sqw_cl, index_q, index_w, sqw_cl)) 914 printf("Isotropic_Sqw: %s: " 915 "Error when setting Sqw_cl[%li q=%g,%li w=%g]=%g from file %s\n", 916 Sqw->compname, index_q, q, index_w, w, sqw_cl, Sqw_Data->filename); 917 sq_cl += sqw_cl; 918 } 919 sq += sqw_full; 920 } /* for index_w */ 921 922 sq *= Sqw_Data->w_step; /* S(q) = \int S(q,w) dw = structure factor */ 923 sq_cl *= Sqw_Data->w_step; 924 /* find minimal reliable q value (not interpolated) */ 925 if (q >= q_min && !q_min_index && sq) { 926 q_min_index = index_q; 927 q_min = q; 928 if (0.9 < sq) 929 S0 = sq; /* minimum reliable S(q) */ 930 else S0 = 1; 931 } 932 /* compute <u^2> = <3 * ln(S(q)) / q^2> */ 933 if (q_min_index && q && S0 && sq) { 934 u2 += 3 * log(sq/S0) /q/q; 935 u2_count++; 936 } 937 938 /* store moment values (q) as M0=S(q) M1=E_r M3 w_c w_l M0_cl=S_cl(q) */ 939 Table_SetElement(&Sqw_moments[0], index_q, 0, sq); 940 Table_SetElement(&Sqw_moments[1], index_q, 0, w1); 941 Table_SetElement(&Sqw_moments[2], index_q, 0, w3); 942 if (w1 > 0 && sq && Sqw->Temperature > 0) { 943 double w_c = sqrt(w1/sq*2*Sqw->Temperature*Sqw->T2E); /* HBAR^2 q^2 kT /m/ S(q) */ 944 Table_SetElement(&Sqw_moments[3], index_q, 0, w_c); /* collective dispersion */ 945 } 946 if (w1 && w3*w1 > 0) { 947 double w_l = sqrt(w3/w1); 948 Table_SetElement(&Sqw_moments[4], index_q, 0, w_l); /* harmonic dispersion */ 949 } 950 if (Sqw->Temperature > 0) 951 Table_SetElement(&Sqw_moments[5], index_q, 0, sq_cl); 952 953 } /* for index_q */ 954 955 956 957 /* display some usefull information, only once in MPI mode (MASTER) */ 958 if (Sqw->Temperature > 0) { 959 double Da = 1.660538921e-27; /* [kg] unified atomic mass unit = Dalton = 1 g/mol */ 960 #ifndef KB 961 double KB = 1.3806503e-23; /* [J/K] */ 962 /* HBAR = 1.05457168e-34 */ 963 #endif 964 double CELE = 1.602176487e-19; /* [C] Elementary charge CODATA 2006 'e' */ 965 double meV2Hz = CELE/HBAR/1000/2/PI; /* 1 meV = 241.80e9 GHz */ 966 double gqw_sum = 0; 967 968 /* classical Sqw */ 969 sprintf(c, "%s_%s_cl.sqw", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 970 Table_Write(Sqw_cl, c, "Momentum [Angs-1]", "'S(q,w)*exp(hw/2kT) classical limit' Energy [meV]", 971 0,Sqw_Data->q_max,-Sqw_Data->w_max,Sqw_Data->w_max); 972 Table_Free(&Sqw_cl); 973 974 if (u2_count) u2 /= u2_count; 975 976 MPI_MASTER( 977 if (do_coh || do_inc) 978 printf("Isotropic_Sqw: %s: " 979 "Physical constants from the S(q,w) %s for T=%g [K]. Values are estimates.\n", 980 Sqw->compname, Sqw_Data->filename, Sqw->Temperature); 981 if (do_coh) { 982 if (Sqw->mat_weight) { 983 double LAMBDA = HBAR*2*PI/sqrt(2*PI*Sqw->mat_weight*Da*KB*Sqw->Temperature)*1e10; /* in [Angs] */ 984 double z = Sqw->mat_rho * LAMBDA*LAMBDA*LAMBDA; /* fugacity , rho=N/V in [Angs-3]*/ 985 double mu = KB*Sqw->Temperature*log(z); /* perfect gas chemical potential */ 986 printf("# De Broglie wavelength LAMBDA=%g [Angs]\n", LAMBDA); 987 printf("# Fugacity z=%g (from Egelstaff p32 Eq 2.31)\n", z); 988 printf("# Chemical potential mu=%g [eV] (eq. perfect gas)\n", mu/CELE); 989 } 990 991 /* compute isothermal sound velocity and compressibility */ 992 /* get the S(q_min) value and the corresponding w_c */ 993 994 if (q_min_index > 0 && q_min && q_min < 0.6) { 995 double w_c = Table_Index(Sqw_moments[3], q_min_index, 0); /* meV */ 996 /* HBAR = [J*s] */ 997 double c_T = 2*PI*w_c*meV2Hz/q_min/1e10; /* meV*Angs -> m/s */ 998 double ChiT= S0/(KB*Sqw->Temperature*Sqw->mat_rho*1e30); 999 printf("# Isothermal compressibility Chi_T=%g [Pa-1] (Egelstaff p201 Eq 10.21) at q=%g [Angs-1]\n", 1000 ChiT, q_min); 1001 printf("# Isothermal sound velocity c_T=%g [m/s] (Lovesey T1 p210 Eq 5.197) at q=%g [Angs-1]\n", 1002 c_T, q_min); 1003 1004 /* Computation if C11 is rather tricky as it is obtained from w_l, which is usually quite noisy 1005 * This means that the obtained values are not reliable from C = rho c_l^2 (Egelstaff Eq 14.10b p284) 1006 * C44 = rho c_c^2 ~ C11/3 1007 */ 1008 double w_l = Table_Index(Sqw_moments[4], q_min_index, 0); /* meV */ 1009 double c_l = 2*PI*w_l*meV2Hz/q_min/1e10; /* meV*Angs -> m/s */ 1010 double C11 = (Sqw->mat_weight*Da)*(Sqw->mat_rho*1e30)*c_l*c_l; 1011 printf("# Elastic modulus C11=%g [GPa] (Egelstaff Eq 14.10b p284) [rough estimate] at q=%g [Angs-1]\n", 1012 C11/1e9, q_min); 1013 } 1014 } 1015 if (do_inc) { 1016 /* display the mean square displacement from S(q) = exp(-<u^2>q^2/3) 1017 <u^2>= <3 * ln(S(q)) / q^2> 1018 */ 1019 if (u2_count && u2) { 1020 printf("# Mean square displacement <u^2>=%g [Angs^2] (<3 * ln(S(q)) / q^2>)\n", u2); 1021 } 1022 1023 /* compute the mean diffusion coefficient D=w_c/q^2 */ 1024 /* FWHM of gaussian is Gamma*RMS2FWHM, only in diffusive regime (Q < 0.2 Angs-1) */ 1025 if (q_min_index > 0 && q_min && q_min < 0.6) { 1026 double w_c = Table_Index(Sqw_moments[3], q_min_index, 0); 1027 double D = 2*PI*w_c*meV2Hz/q_min/q_min/1e14*RMS2FWHM/2; /* meV*Angs^2 -> mm^2/s */ 1028 printf("# Diffusion coefficient D=%g [mm^2/s] (Egelstaff p220)\n", D); 1029 if (u2_count && u2 && D) 1030 printf("# Jump relaxation time tau=%g [ns] (Egelstaff Eq 11.8 p220)\n", u2*1e-2/6/D); 1031 } 1032 } 1033 ); /* end MASTER in MPI mode */ 1034 1035 /* density of states (generalized) */ 1036 if (!Table_Init(&Gqw, Sqw_Data->q_bins, Sqw_Data->w_bins)) { 1037 printf("Isotropic_Sqw: %s: Cannot allocate G(q,w) Table (%lix%i).\n" 1038 "WARNING Skipping S(q,w) diagnosis.\n", 1039 Sqw->compname, Sqw_Data->q_bins, 1); 1040 return; 1041 } 1042 sprintf(Gqw.filename, 1043 "G(q,w) from %s (generalized density of states, Carpenter J Non Cryst Sol 92 (1987) 153)", 1044 Sqw_Data->filename); 1045 Gqw.block_number = 1; 1046 Gqw.min_x = 0; 1047 Gqw.max_x = Sqw_Data->q_max; 1048 Gqw.step_x = Sqw_Data->q_step; 1049 1050 for (index_w=0; index_w < Sqw_Data->w_bins; index_w++) { 1051 double w = -Sqw_Data->w_max + index_w*Sqw_Data->w_step; /* w value in Sqw_full */ 1052 double gw = 0; 1053 for (index_q=0; index_q < Sqw_Data->q_bins; index_q++) { 1054 double q = index_q*Sqw_Data->q_step; /* q value in Sqw_full ; q_min = 0 */ 1055 double sqw_full = Table_Index(Sqw_Data->Sqw, index_q, index_w); 1056 double n = 1/(exp(w/(Sqw->Temperature*Sqw->T2E))-1); /* Bose factor */ 1057 double DW = q && u2 ? exp(2*u2*q*q/6) : 1; /* Debye-Weller factor */ 1058 double gqw = q && n+1 ? sqw_full*DW*2*(Sqw->mat_weight*Da)*w/(n+1)/q/q : 0; 1059 if (!Table_SetElement(&Gqw, index_q, index_w, gqw)) 1060 printf("Isotropic_Sqw: %s: " 1061 "Error when setting Gqw[%li q=%g,%li w=%g]=%g from file %s\n", 1062 Sqw->compname, index_q, q, index_w, w, gqw, Sqw_Data->filename); 1063 gw += gqw; 1064 gqw_sum += gqw; 1065 } 1066 Table_SetElement(&Sqw_moments[6], index_w, 0, gw); 1067 } 1068 1069 /* normalize the density of states */ 1070 for (index_w=0; index_w < Sqw_Data->w_bins; index_w++) { 1071 double gw = Table_Index(Sqw_moments[6], index_w, 0); 1072 Table_SetElement(&Sqw_moments[6], index_w, 0, gw / gqw_sum); 1073 for (index_q=0; index_q < Sqw_Data->q_bins; index_q++) { 1074 double gqw = Table_Index(Gqw, index_q, index_w); 1075 Table_SetElement(&Gqw, index_q, index_w, gqw / gqw_sum); 1076 } 1077 } 1078 1079 /* write Gqw and free memory */ 1080 if (Sqw_Data->w_bins > 1) { 1081 sprintf(c, "%s_%s.gqw", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 1082 Table_Write(Gqw, c, "Momentum [Angs-1]", "'Generalized density of states' Energy [meV]", 1083 0,Sqw_Data->q_max,-Sqw_Data->w_max,Sqw_Data->w_max); 1084 Table_Free(&Gqw); 1085 } 1086 } /* if T>0 */ 1087 1088 /* write all tables to disk M0=S(q) M1=E_r M3 w_c w_l M0_cl=S_cl(q) */ 1089 if (Sqw_Data->w_bins > 1) { 1090 sprintf(c, "%s_%s.m1", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 1091 Table_Write(Sqw_moments[1], c, "Momentum [Angs-1]", "int w S(q,w) dw (recoil) q^2/2m [meV]", 1092 0,Sqw_Data->q_max,0,0); 1093 sprintf(c, "%s_%s.w_l", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 1094 Table_Write(Sqw_moments[4], c, "Momentum [Angs-1]", "w_l(q) harmonic frequency [meV]", 1095 0,Sqw_Data->q_max,0,0); 1096 sprintf(c, "%s_%s.sqw", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 1097 Table_Write(Sqw_Data->Sqw, c, "Momentum [Angs-1]", "'S(q,w) dynamical structure factor [meV-1]' Energy [meV]", 1098 0,Sqw_Data->q_max,-Sqw_Data->w_max,Sqw_Data->w_max); 1099 1100 if (Sqw->Temperature > 0) { 1101 sprintf(c, "%s_%s.w_c", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 1102 Table_Write(Sqw_moments[3], c, "Momentum [Angs-1]", "w_c(q) collective excitation [meV]", 0,Sqw_Data->q_max,0,0); 1103 sprintf(c, "%s_%s_cl.sq", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 1104 Table_Write(Sqw_moments[5], c, "Momentum [Angs-1]", "int S_cl(q,w) dw", 1105 0,Sqw_Data->q_max,0,0); 1106 sprintf(c, "%s_%s.gw", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 1107 Table_Write(Sqw_moments[6], c, "Energy [meV]", "'Generalized effective density of states' Energy [meV]", 1108 -Sqw_Data->w_max,Sqw_Data->w_max,0,0); 1109 1110 } 1111 } 1112 sprintf(c, "%s_%s.sq", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 1113 Table_Write(Sqw_moments[0], c, "Momentum [Angs-1]","S(q) = int S(q,w) dw", 0,Sqw_Data->q_max,0,0); 1114 sprintf(c, "%s_%s.sigma", Sqw->compname, Sqw_Data->type=='c' ? "coh" : "inc"); 1115 Table_Write(Sqw_Data->iqSq, c, "Energy [meV]", "sigma kf/ki int q S(q,w) dw scattering cross section [barns]", 0,0,0,0); 1116 1117 /* free Tables */ 1118 for (index_q=0; index_q < 7; Table_Free(&Sqw_moments[index_q++])); 1119 1120} /* Sqw_diagnosis */ 1121 1122/***************************************************************************** 1123* Sqw_readfile: Read Sqw data files 1124* Returns Sqw_Data_struct or NULL in case of error 1125* Used in : Sqw_init (2) 1126*****************************************************************************/ 1127struct Sqw_Data_struct *Sqw_readfile( 1128 struct Sqw_sample_struct *Sqw, char *file, struct Sqw_Data_struct *Sqw_Data) 1129{ 1130 1131 t_Table *Table_Array= NULL; 1132 long nblocks = 0; 1133 char flag = 0; 1134 1135 t_Table Sqw_full, iqSq; /* the Sqw (non symmetric) and total scattering X section */ 1136 1137 double sum=0; 1138 double mat_at_nb=1; 1139 double iq2Sq=0; 1140 long *SW_lookup=NULL; 1141 long **QW_lookup=NULL; 1142 char **parsing =NULL; 1143 1144 long index_q, index_w; 1145 double q_min_file, q_max_file, q_step_file; 1146 long q_bins_file; 1147 double w_min_file, w_max_file, w_step_file; 1148 long w_bins_file; 1149 double q_max, q_step; 1150 long q_bins; 1151 double w_max, w_step; 1152 long w_bins; 1153 1154 double alpha=0; 1155 1156 double M1 = 0; 1157 double M1_cl = 0; 1158 double T = 0; 1159 double T_file = 0; 1160 long T_count = 0; 1161 long M1_count = 0; 1162 long M1_cl_count = 0; 1163 1164 /* setup default */ 1165 Sqw_Data_init(Sqw_Data); 1166 1167 if (!file || !strlen(file) || !strcmp(file, "NULL") || !strcmp(file, "0")) return(Sqw_Data); 1168 /* read the Sqw file */ 1169 Table_Array = Table_Read_Array(file, &nblocks); 1170 strncpy(Sqw_Data->filename, file, 80); 1171 if (!Table_Array) return(NULL); 1172 1173 /* (1) parsing of header ================================================== */ 1174 parsing = Table_ParseHeader(Table_Array[0].header, 1175 "Vc","V_0", 1176 "sigma_abs","sigma_a ", 1177 "sigma_inc","sigma_i ", 1178 "column_j", /* 6 */ 1179 "column_d", 1180 "column_F2", 1181 "column_DW", 1182 "column_Dd", 1183 "column_inv2d", "column_1/2d", "column_sintheta_lambda", 1184 "column_q", /* 14 */ 1185 "sigma_coh","sigma_c ", 1186 "Temperature", 1187 "column_Sq", 1188 "column_F ", /* 19 */ 1189 "V_rho", 1190 "density", 1191 "weight", 1192 "nb_atoms","multiplicity", 1193 "classical", 1194 NULL); 1195 if (parsing) { 1196 int i; 1197 if (parsing[0] && !Sqw->mat_rho) Sqw->mat_rho =1/atof(parsing[0]); 1198 if (parsing[1] && !Sqw->mat_rho) Sqw->mat_rho =1/atof(parsing[1]); 1199 if (parsing[2] && !Sqw->s_abs) Sqw->s_abs = atof(parsing[2]); 1200 if (parsing[3] && !Sqw->s_abs) Sqw->s_abs = atof(parsing[3]); 1201 if (parsing[4] && !Sqw->s_inc) Sqw->s_inc = atof(parsing[4]); 1202 if (parsing[5] && !Sqw->s_inc) Sqw->s_inc = atof(parsing[5]); 1203 if (parsing[6]) Sqw->column_order[0]=atoi(parsing[6]); 1204 if (parsing[7]) Sqw->column_order[1]=atoi(parsing[7]); 1205 if (parsing[8]) Sqw->column_order[2]=atoi(parsing[8]); 1206 if (parsing[9]) Sqw->column_order[3]=atoi(parsing[9]); 1207 if (parsing[10]) Sqw->column_order[4]=atoi(parsing[10]); 1208 if (parsing[11]) Sqw->column_order[5]=atoi(parsing[11]); 1209 if (parsing[12]) Sqw->column_order[5]=atoi(parsing[12]); 1210 if (parsing[13]) Sqw->column_order[5]=atoi(parsing[13]); 1211 if (parsing[14]) Sqw->column_order[6]=atoi(parsing[14]); 1212 if (parsing[15] && !Sqw->s_coh) Sqw->s_coh=atof(parsing[15]); 1213 if (parsing[16] && !Sqw->s_coh) Sqw->s_coh=atof(parsing[16]); 1214 if (parsing[17] && !Sqw->Temperature) Sqw->Temperature=atof(parsing[17]); /* from user or file */ 1215 if (parsing[17] && !T_file) T_file=atof(parsing[17]); /* from file */ 1216 if (parsing[18]) Sqw->column_order[8]=atoi(parsing[18]); 1217 if (parsing[19]) Sqw->column_order[7]=atoi(parsing[19]); 1218 if (parsing[20] && !Sqw->mat_rho) Sqw->mat_rho =atof(parsing[20]); 1219 if (parsing[21] && !Sqw->mat_density) Sqw->mat_density=atof(parsing[21]); 1220 if (parsing[22] && !Sqw->mat_weight) Sqw->mat_weight =atof(parsing[22]); 1221 if (parsing[23] ) mat_at_nb =atof(parsing[23]); 1222 if (parsing[24] ) mat_at_nb =atof(parsing[24]); 1223 if (parsing[25] ) { /* classical is found in the header */ 1224 char *endptr; 1225 double value = strtod(parsing[25], &endptr); 1226 if (*endptr == *parsing[25]) { 1227 if (Sqw->sqw_classical < 0) Sqw->sqw_classical = 1; 1228 } else Sqw->sqw_classical = value; 1229 } 1230 for (i=0; i<=25; i++) if (parsing[i]) free(parsing[i]); 1231 free(parsing); 1232 } 1233 1234 /* compute the scattering unit density from material weight and density */ 1235 /* the weight of the scattering element is the chemical formula molecular weight 1236 * times the nb of chemical formulae in the scattering element (nb_atoms) */ 1237 if (!Sqw->mat_rho && Sqw->mat_density > 0 && Sqw->mat_weight > 0 && mat_at_nb > 0) { 1238 /* molar volume [cm^3/mol] = weight [g/mol] / density [g/cm^3] */ 1239 /* atom density per Angs^3 = [mol/cm^3] * N_Avogadro *(1e-8)^3 */ 1240 Sqw->mat_rho = Sqw->mat_density/(Sqw->mat_weight*mat_at_nb)/1e24*NA; 1241 MPI_MASTER( 1242 if (Sqw->verbose_output > 0) 1243 printf("Isotropic_Sqw: %s: Computing scattering unit density V_rho=%g [AA^-3] from density=%g [g/cm^3] weight=%g [g/mol].\n", 1244 Sqw->compname, Sqw->mat_rho, Sqw->mat_density, Sqw->mat_weight); 1245 ); 1246 } 1247 1248 /* the scattering unit cross sections are the chemical formula ones 1249 * times the nb of chemical formulae in the scattering element */ 1250 if (mat_at_nb > 0) { 1251 Sqw->s_abs *= mat_at_nb; Sqw->s_inc *= mat_at_nb; Sqw->s_coh *= mat_at_nb; 1252 } 1253 1254 if (nblocks) { 1255 if (nblocks == 1) { 1256 /* import Powder file */ 1257 t_Table *newTable = NULL; 1258 newTable = Sqw_read_PowderN(Sqw, Table_Array[0]); 1259 if (!newTable) { 1260 printf("Isotropic_Sqw: %s: ERROR importing powder line file %s.\n" 1261 " Check format definition.\n", 1262 Sqw->compname, file); 1263 exit(-1); 1264 } else flag=0; 1265 Table_Free_Array(Table_Array); 1266 Table_Array = newTable; 1267 } else if (nblocks != 3) { 1268 printf("Isotropic_Sqw: %s: ERROR " 1269 "File %s contains %li block%s instead of 3.\n", 1270 Sqw->compname, file, nblocks, (nblocks == 1 ? "" : "s")); 1271 } else { flag=0; Sqw->barns=0; /* Sqw files do not use powder_barns */ } 1272 } 1273 1274 /* print some info about Sqw files */ 1275 if (flag) Sqw->verbose_output = 2; 1276 1277 if (flag) { 1278 if (nblocks) printf("ERROR Wrong file format.\n" 1279 " Disabling contribution.\n" 1280 " File must contain 3 blocks for [q,w,sqw] or Powder file (1 block, laz,lau).\n"); 1281 return(Sqw_Data); 1282 } 1283 1284 sprintf(Table_Array[0].filename, "%s#q", file); 1285 sprintf(Table_Array[1].filename, "%s#w", file); 1286 sprintf(Table_Array[2].filename, "%s#sqw", file); 1287 1288 MPI_MASTER( 1289 if (nblocks && Sqw->verbose_output > 2) { 1290 printf("Isotropic_Sqw: %s file read, analysing...\n", file); 1291 Table_Info_Array(Table_Array); 1292 } 1293 ); 1294 1295 /* (2) compute range for full +/- w and allocate S(q,w) =================== */ 1296 1297 /* get the q,w extend of the table from the file */ 1298 q_bins_file = Table_Array[0].rows*Table_Array[0].columns; 1299 w_bins_file = Table_Array[1].rows*Table_Array[1].columns; 1300 1301 /* is there enough qw data in file to proceed ? */ 1302 if (q_bins_file <= 1 || w_bins_file <= 0) { 1303 printf("Isotropic_Sqw: %s: Data file %s has incomplete q or omega information (%lix%li).\n" 1304 "ERROR Exiting.\n", 1305 Sqw->compname, file, q_bins_file, w_bins_file); 1306 return(Sqw_Data); 1307 } 1308 1309 q_min_file = Table_Array[0].min_x; q_max_file = Table_Array[0].max_x; 1310 q_step_file = Table_Array[0].step_x ? Table_Array[0].step_x : (q_max_file - q_min_file)/(Table_Array[0].rows*Table_Array[0].columns); 1311 w_min_file = Table_Array[1].min_x; w_max_file = Table_Array[1].max_x; 1312 w_step_file = Table_Array[1].step_x; 1313 1314 /* create a regular extended q,w and Sqw tables applying the exp(-hw/kT) factor */ 1315 q_max = q_max_file; 1316 q_bins = (q_step_file ? q_max/q_step_file : q_bins_file)+1; 1317 q_step = q_bins-1 > 0 ? q_max/(q_bins-1) : 1; 1318 w_max = fabs(w_max_file); 1319 if (fabs(w_min_file) > fabs(w_max_file)) w_max = fabs(w_min_file); 1320 /* w_min =-w_max */ 1321 w_bins = (w_step_file ? (long)(2*w_max/w_step_file) : 0)+1; /* twice the initial w range */ 1322 w_step = w_bins-1 > 0 ? 2*w_max/(w_bins-1) : 1; /* that is +/- w_max */ 1323 1324 /* create the Sqw table in full range */ 1325 if (!Table_Init(&Sqw_full, q_bins, w_bins)) { 1326 printf("Isotropic_Sqw: %s: Cannot allocate Sqw_full Table (%lix%li).\n" 1327 "ERROR Exiting.\n", 1328 Sqw->compname, q_bins, w_bins); 1329 return(NULL); 1330 } 1331 sprintf(Sqw_full.filename, "S(q,w) from %s (dynamic structure factor)", file); 1332 Sqw_full.block_number = 1; 1333 1334 Sqw_Data->q_bins = q_bins; Sqw_Data->q_max = q_max; Sqw_Data->q_step= q_step; 1335 Sqw_Data->w_bins = w_bins; Sqw_Data->w_max = w_max; Sqw_Data->w_step= w_step; 1336 Sqw_Data->q_min_file = q_min_file; 1337 1338 /* build an energy symmetric Sqw data set with detailed balance there-in, so 1339 * that we can both compute effective scattering Xsection, probability distributions 1340 * that is S(q) and \int q S(q). 1341 * We scan the new Sqw table elements with regular qw binning and search for their 1342 * equivalent element in the Sqw file data set. This is slower than doing the opposite. 1343 * We could be scanning all file elements, and fill the new table, but in the 1344 * process some empty spaces may appear when the initial file binning is not regular 1345 * in qw, leading to gaps in the new table. 1346 */ 1347 1348 /* (3) we build q and w lookup table for conversion file -> sqw_full ====== */ 1349 MPI_MASTER( 1350 if (Sqw->verbose_output > 2) 1351 printf("Isotropic_Sqw: %s: Creating Sqw_full... (%s, %s)\n", 1352 Sqw->compname, file, Sqw->type=='c' ? "coh" : "inc"); 1353 ); 1354 1355 double w_file2full[w_bins]; /* lookup table for fast file -> Sqw_full allocation */ 1356 1357 for (index_w=0; index_w < w_bins; w_file2full[index_w++]=0); 1358 1359 for (index_w=0; index_w < w_bins; index_w++) { 1360 1361 double w = -w_max + index_w*w_step; /* w value in Sqw_full */ 1362 double index_w_file=0; /* w index in Sqw file */ 1363 char found=0; 1364 for (index_w_file=0; index_w_file < w_bins_file; index_w_file++) { 1365 double w0=Table_Index(Table_Array[1], (long)index_w_file, 0); 1366 double w1=Table_Index(Table_Array[1], (long)index_w_file+1,0); 1367 /* test if we are in Stokes */ 1368 if (w0 > w1) { 1369 double tmp=w0; w0=w1; w1=tmp; 1370 } 1371 if (w0 <= w && w < w1) { 1372 /* w ~ w_file exists in file, usually on w > 0 side Stokes, neutron looses energy */ 1373 index_w_file += w1-w0 ? (w-w0)/(w1-w0) : 0; /* may correspond with a position in-betwwen two w elements */ 1374 found=1; 1375 break; 1376 } 1377 } 1378 /* test if we are in anti-Stokes */ 1379 if (!found) 1380 for (index_w_file=0; index_w_file < w_bins_file; index_w_file++) { 1381 double w0=Table_Index(Table_Array[1], (long)index_w_file, 0); 1382 double w1=Table_Index(Table_Array[1], (long)index_w_file+1,0); 1383 /* test if we are in anti-Stokes */ 1384 if (w0 > w1) { 1385 double tmp=w0; w0=w1; w1=tmp; 1386 } 1387 if (w0 <= -w && -w < w1) { /* w value is mirrored from the opposite side in file */ 1388 index_w_file += w1-w0 ? (-w-w0)/(w1-w0) : 0; 1389 index_w_file = -index_w_file; /* in this case, index value is set to negative */ 1390 break; 1391 } 1392 } 1393 w_file2full[index_w] = index_w_file; 1394 } 1395 1396 double q_file2full[q_bins]; 1397 for (index_q=0; index_q < q_bins; q_file2full[index_q++]=0); 1398 1399 for (index_q=0; index_q < q_bins; index_q++) { 1400 1401 double q = index_q*q_step; /* q value in Sqw_full ; q_min = 0 */ 1402 double index_q_file= 0; /* q index in Sqw file */ 1403 1404 /* search for q value in the initial file data set */ 1405 if (q <= q_min_file) index_q_file=0; 1406 else if (q >= q_max_file) index_q_file=q_bins_file-1; 1407 else 1408 for (index_q_file=0; index_q_file < q_bins_file; index_q_file++) { 1409 double q0=Table_Index(Table_Array[0], (long)index_q_file, 0); 1410 double q1=Table_Index(Table_Array[0], (long)index_q_file+1,0); 1411 if (q0 <= q && q <= q1) { 1412 index_q_file += q1-q0 ? (q-q0)/(q1-q0) : 0; /* may correspond with a position in-betwwen two q elements */ 1413 break; 1414 } 1415 } 1416 q_file2full[index_q] = index_q_file; 1417 } 1418 1419 /* (4) now we build Sqw on full Q,W ranges, using the Q,W table lookup above */ 1420 for (index_q=0; index_q < q_bins; index_q++) { 1421 1422 double q = index_q*q_step; /* q value in Sqw_full ; q_min = 0 */ 1423 double index_q_file= 0; /* q index in Sqw file */ 1424 1425 /* get q value in the initial file data set */ 1426 index_q_file = q_file2full[index_q]; 1427 1428 /* now scan energy elements in Sqw full, and search these in file data */ 1429 for (index_w=0; index_w < w_bins; index_w++) { 1430 double w = -w_max + index_w*w_step; /* w value in Sqw_full */ 1431 double index_w_file=0; /* w index in Sqw file */ 1432 double sqw_file =0; /* Sqw(index_q, index_w) value interpolated from file */ 1433 1434 /* search for w value in the file data set, negative when mirrored */ 1435 index_w_file = w_file2full[index_w]; 1436 /* get Sqw_file element, with bi-linear interpolation from file */ 1437 /* when the initial file does not contain the energy, the opposite element (-w) is used */ 1438 sqw_file = Table_Value2d(Table_Array[2], index_q_file, abs(index_w_file)); 1439 /* apply the minimum threshold to remove noisy background in S(q,w) */ 1440 if (sqw_file < Sqw->sqw_threshold) sqw_file = 0; 1441 else if (index_w_file < 0) sqw_file = -sqw_file; /* negative == mirrored from other side */ 1442 1443 if (!Table_SetElement(&Sqw_full, index_q, index_w, sqw_file)) 1444 printf("Isotropic_Sqw: %s: " 1445 "Error when setting Sqw[%li q=%g,%li w=%g]=%g from file %s\n", 1446 Sqw->compname, index_q, q, index_w, w, fabs(sqw_file), file); 1447 } /* for index_w */ 1448 } /* for index_q */ 1449 1450 /* free memory and store limits for new full Sqw table */ 1451 Table_Free_Array(Table_Array); 1452 1453 /* if only one S(q,w) side is given, it is symmetrised by mirroring, then M1=0 */ 1454 1455 /* (5) test if the Sqw_full is classical or not by computing the 1st moment (=0 for classical) */ 1456 /* also compute temperature (quantum case) if not set */ 1457 for (index_q=0; index_q < q_bins; index_q++) { 1458 1459 double q = index_q*q_step; /* q value in Sqw_full ; q_min = 0 */ 1460 1461 for (index_w=0; index_w < w_bins; index_w++) { 1462 double w = -w_max + index_w*w_step; /* w value in Sqw_full */ 1463 double sqw_full = Table_Index(Sqw_full, index_q, index_w); 1464 long index_mw = w_bins-1-index_w; /* opposite w index in S(q,w) */ 1465 double sqw_opp = Table_Index(Sqw_full, index_q, index_mw); 1466 1467 /* the analysis must be done only on values which exist on both sides */ 1468 /* as integrals must be symmetric, and Bose factor requires both sides as well */ 1469 if (sqw_full > 0 && sqw_opp > 0) { 1470 /* compute temperature from Bose factor */ 1471 if (sqw_opp != sqw_full) { 1472 T += fabs(w/log(sqw_opp/sqw_full)/Sqw->T2E); 1473 T_count++; 1474 } 1475 /* we first assume Sqw is quantum. M1_cl should be 0, M1 should be recoil */ 1476 M1 += w*sqw_full*w_step; 1477 M1_count++; 1478 /* we assume it is quantum (non symmetric) and check that its symmetrized version has M1_cl=0 */ 1479 if (Sqw->Temperature > 0) { 1480 sqw_opp = sqw_full * Sqw_quantum_correction(-w,Sqw->Temperature,Sqw->Q_correction); /* Sqw_cl */ 1481 M1_cl += w*sqw_opp*w_step; 1482 M1_cl_count++; 1483 } else if (Sqw->mat_weight) { 1484 /* T=0 ? would compute the M1_cl = M1 - recoil energy, assuming we have a quantum S(q,w) in file */ 1485 /* the M1(quantum) = (MNEUTRON/m)*2.0725*q^2 recoil energy */ 1486 double Da = 1.660538921e-27; /* atomic mass unit */ 1487 double Er = (MNEUTRON/Sqw->mat_weight/Da)*2.0725*q*q; /* recoil for one scattering unit in the cell [meV] Schober JDN16 p239 */ 1488 M1_cl += M1 - Er; 1489 M1_cl_count++; 1490 } 1491 } /* both side from file */ 1492 } /*index_w */ 1493 } /*index_q */ 1494 1495 if (T_count) T /= T_count; /* mean temperature from Bose ratio */ 1496 if (M1_count) M1 /= M1_count; 1497 if (M1_cl_count) M1_cl /= M1_cl_count; /* mean energy value along q range */ 1498 1499 /* determine if we use a classical or quantum S(q,w) */ 1500 if (Sqw->sqw_classical < 0) { 1501 if (fabs(M1) < 2*w_step) { 1502 Sqw->sqw_classical = 1; /* the initial Sqw from file seems to be centered, thus classical */ 1503 } else if (fabs(M1_cl) < fabs(M1)) { 1504 /* M1 for classical is closer to 0 than for quantum one */ 1505 Sqw->sqw_classical = 0; /* initial data from file seems to be quantum (non classical) */ 1506 } else { /* M1_cl > M1 > 2*w_step */ 1507 printf("Isotropic_Sqw: %s: I do not know if S(q,w) data is classical or quantum.\n" 1508 "WARNING First moment M1=%g M1_cl=%g for file %s. Defaulting to classical case.\n", 1509 Sqw->compname, M1, M1_cl, file); 1510 } 1511 } 1512 if (Sqw->sqw_classical < 0) Sqw->sqw_classical=1; /* default when quantum/classical choice is not set */ 1513 1514 if (!T_file && T) { 1515 if (Sqw->Temperature < 0) Sqw->Temperature = fabs(T); 1516 MPI_MASTER( 1517 if (Sqw->verbose_output > 0 && T < 3000) { 1518 printf( "Isotropic_Sqw: %s: Temperature computed from S(q,w) data from %s is T=%g [K] (not set here).\n", 1519 Sqw->compname, file, T); 1520 if (Sqw->Temperature == 0) 1521 printf("Warning: %s: Use T=-1 to set it. Currently using T=%g, i.e. no detailed balance.\n", 1522 Sqw->compname, Sqw->Temperature); 1523 } 1524 if (!Sqw->sqw_classical && Sqw->Temperature>0 && fabs(Sqw->Temperature - T) > 0.1*Sqw->Temperature) 1525 printf("WARNING: %s: The temperature %g [K] guessed from the non-classical\n" 1526 " S(q,w) %s does not match the requested T=%g [K].\n", 1527 Sqw->compname, T, file, Sqw->Temperature); 1528 ); 1529 } else T=T_file; 1530 1531 MPI_MASTER( 1532 if (Sqw->verbose_output > 0 && w_bins > 1) 1533 printf("Isotropic_Sqw: %s: S(q,w) data from %s (%s) assumed to be %s.\n", 1534 Sqw->compname, file, Sqw->type=='c' ? "coh" : "inc", 1535 Sqw->sqw_classical ? "classical (symmetrised in energy)" : "non-classical (includes Bose factor, non symmetric in energy)"); 1536 ); 1537 1538 /* (6) apply detailed balance on Sqw_full for classical or quantum S(q,w) */ 1539 /* compute the \int q^2 S(q) for normalisation */ 1540 1541 /* temperatures defined 1542 * T_file: temperature read from the .sqw file 1543 * T: temperature computed from the quantum Sqw using detailed balance 1544 * Sqw->Temperature: temperature defined by user, or read from file when not set 1545 */ 1546 1547 MPI_MASTER( 1548 if (Sqw->sqw_classical && Sqw->verbose_output > 0 && Sqw->Temperature > 0) 1549 printf("Isotropic_Sqw: %s: Applying exp(hw/2kT) factor with T=%g [K] on %s file (classical/symmetric) using '%s' quantum correction\n", 1550 Sqw->compname, Sqw->Temperature, file, Sqw->Q_correction); 1551 ); 1552 for (index_q=0; index_q < q_bins; index_q++) { 1553 double sq = 0; 1554 double q = index_q*q_step; /* q value in Sqw_full ; q_min = 0 */ 1555 for (index_w=0; index_w < w_bins; index_w++) { 1556 double w = -w_max + index_w*w_step; /* w value in Sqw_full */ 1557 double balance = 1; /* detailed balance factor, default is 1 */ 1558 double sqw_full = Table_Index(Sqw_full, index_q, index_w); 1559 1560 /* do we use a symmetric S(q,w) from real G(r,t) from e.g. MD ? */ 1561 1562 if (Sqw->sqw_classical && Sqw->Temperature > 0) /* data is symmetric, we apply Bose factor */ 1563 /* un-symmetrize Sqw(file) * exp(hw/kT/2) on both sides */ 1564 balance = Sqw_quantum_correction(w, Sqw->Temperature, Sqw->Q_correction); 1565 else if (!Sqw->sqw_classical) { /* data is quantum (contains Bose) */ 1566 if (sqw_full < 0) { /* quantum but mirrored/symmetric value (was missing in file) */ 1567 if (T) 1568 /* prefer to use T computed from file for mirroring */ 1569 balance *= exp(w/(T*Sqw->T2E)); /* apply Bose where missing */ 1570 else if (Sqw->Temperature) 1571 balance *= exp(w/(Sqw->Temperature*Sqw->T2E)); /* apply Bose where missing */ 1572 } 1573 /* test if T computed from file matches requested T, else apply correction */ 1574 if (T && Sqw->Temperature > 0 && Sqw->Temperature != T) { 1575 balance *= exp(-w/(T*Sqw->T2E)/2); /* make it a classical data set: remove computed/read T from quantum data file */ 1576 balance *= exp( w/(Sqw->Temperature*Sqw->T2E)/2); /* then apply Bose to requested temperature */ 1577 } 1578 } 1579 1580 /* update Sqw value */ 1581 sqw_full = fabs(sqw_full)*balance; 1582 Table_SetElement(&Sqw_full, index_q, index_w, sqw_full); 1583 /* sum up the S(q) (0-th moment) = w0 */ 1584 sq += sqw_full; 1585 } /* index_w */ 1586 sq *= w_step; /* S(q) = \int S(q,w) dw = structure factor */ 1587 iq2Sq += q*q*sq*q_step; /* iq2Sq = \int q^2 S(q) dq = used for g-sum rule (normalisation) */ 1588 sum += sq*q_step; /* |S| = \int S(q,w) dq dw = total integral value in file */ 1589 } /* index_q */ 1590 1591 if (!sum) { 1592 printf("Isotropic_Sqw: %s: No valid data in the selected (Q,w) range: sum(S)=0\n" 1593 "ERROR Available Sqw data is\n", 1594 Sqw->compname); 1595 printf(" q=[%g:%g] w=[%g:%g]\n", 1596 q_min_file, q_max_file, 1597 w_min_file, w_max_file); 1598 Table_Free(&Sqw_full); 1599 return(NULL); 1600 } 1601 1602 /* norm S(q ,w) = sum(S)*q_range/q_bins on total q,w range from file */ 1603 sum *= (q_max_file - q_min_file)/q_bins_file; 1604 1605 /* (7) renormalization of S(q,w) ========================================== */ 1606 1607 if ( Sqw->sqw_norm >0) alpha=Sqw->sqw_norm; 1608 else if (!Sqw->sqw_norm) alpha=1; 1609 1610 if (!alpha && iq2Sq) { /* compute theoretical |S| norm */ 1611 /* Eq (2.44) from H.E. Fischer et al, Rep. Prog. Phys., 69 (2006) 233 */ 1612 alpha = 1613 (q_max*q_max*q_max/3 - (Sqw->type == 'c' ? 2*PI*PI*Sqw->mat_rho : 0)) 1614 /iq2Sq; 1615 } 1616 1617 if (alpha < 0) { 1618 MPI_MASTER( 1619 printf("Isotropic_Sqw: %s: normalisation factor is negative. rho=%g [Angs^-3] may be too high.\n" 1620 "WARNING Disabling renormalization i.e. keeping initial S(q,w).\n", 1621 Sqw->compname, Sqw->mat_rho); 1622 ); 1623 alpha=0; 1624 } 1625 1626 /* apply normalization on S(q,w) */ 1627 if (alpha && alpha != 1) { 1628 sum *= alpha; 1629 for (index_q=0; index_q < q_bins ; index_q++) { 1630 for (index_w=0; index_w < w_bins; index_w++) 1631 Table_SetElement(&Sqw_full, index_q, index_w, 1632 Table_Index(Sqw_full, index_q, index_w) * alpha); 1633 } 1634 } 1635 1636 Sqw_Data->intensity = sum; 1637 1638 Table_Stat(&Sqw_full); 1639 Sqw_full.min_x = 0; 1640 Sqw_full.max_x = q_max; 1641 Sqw_full.step_x = q_step; 1642 1643 MPI_MASTER( 1644 if (Sqw->verbose_output > 0) { 1645 printf("Isotropic_Sqw: %s: Generated %s %scoherent Sqw\n" 1646 " q=[%g:%g Angs-1] w=[%g:%g meV] |S|=%g size=[%lix%li] sigma=%g [barns]\n", 1647 Sqw->compname, file, (Sqw->type == 'i' ? "in" : ""), 1648 q_min_file, q_max_file, 1649 w_min_file, w_max_file, Sqw_Data->intensity, 1650 q_bins, Sqw_Data->w_bins, 1651 (Sqw->type == 'i' ? Sqw->s_inc : Sqw->s_coh)); 1652 if (w_max < 1e-2) 1653 printf(" Mainly elastic scattering.\n"); 1654 if (Sqw->sqw_norm >0 && Sqw->sqw_norm != 1) 1655 printf(" normalization factor S(q,w)*%g (user)\n", alpha); 1656 else if (Sqw->sqw_norm<0) 1657 printf(" normalization factor S(q,w)*%g (auto) \\int q^2 S(q) dq=%g\n", alpha, iq2Sq); 1658 } 1659 ); 1660 1661 /* (8) Compute total cross section ======================================== */ 1662 1663 /* now compute the effective total cross section (Sqw_integrate_iqSq) 1664 sigma(Ei) = sigma/2/Ki^2 * \int q S(q,w) dw dq 1665 * for each incoming neutron energy 0 < Ei < 2*w_max, and 1666 * integration range w=-Ei:w_max and q=Q0:Q1 with 1667 * Q0 = SE2Q*(sqrt(E)-sqrt(E+w)) 1668 * Q1 = SE2Q*(sqrt(E)+sqrt(E+w)) 1669 */ 1670 1671 Sqw_Data->lookup_length = Sqw->lookup_length; 1672 Sqw_Data->iqSq_length = Sqw->lookup_length; 1673 /* this length should be greater when w_max=0 for e.g. elastic only */ 1674 if (w_bins <= 1) Sqw_Data->iqSq_length = q_bins; 1675 1676 if (!Table_Init(&iqSq, Sqw_Data->iqSq_length, 1)) { 1677 /* from 0 to 2*Ki_max */ 1678 printf("Isotropic_Sqw: %s: Cannot allocate [int q S(q,w) dq dw] array (%li bytes).\n" 1679 "ERROR Exiting.\n", 1680 Sqw->compname, Sqw->lookup_length*sizeof(double)); 1681 Table_Free(&Sqw_full); 1682 return(NULL); 1683 } 1684 1685 /* compute the maximum incoming energy that can be handled */ 1686 Sqw_Data->Ei_max = 2*w_max; 1687 1688 /* Checked in different ways in Powder and "proper inelastic" case */ 1689 if (w_step==1) { 1690 /* Powder */ 1691 double Ei_max_q = (q_max*K2V)*(q_max*K2V)*VS2E/2; 1692 if (Ei_max_q > Sqw_Data->Ei_max) Sqw_Data->Ei_max = Ei_max_q; 1693 } else { 1694 /* Proper inelastic */ 1695 /* check if the q-range will limit the integration */ 1696 if ((q_max*K2V)*(q_max*K2V)*VS2E/2 > Sqw_Data->Ei_max) { 1697 /* then scan Ei until we pass q_max */ 1698 for (index_w=0; index_w < Sqw_Data->iqSq_length; index_w++) { 1699 double Ei = index_w*2*w_max/Sqw_Data->iqSq_length; 1700 if ( (Ei > w_max && sqrt(Ei)+sqrt(Ei-w_max) >= q_max/(SE2V*V2K)) 1701 || sqrt(Ei)+sqrt(Ei+w_max) >= q_max/(SE2V*V2K)) 1702 if (Ei < Sqw_Data->Ei_max) { 1703 Sqw_Data->Ei_max = Ei; 1704 break; 1705 } 1706 } 1707 } 1708 } 1709 1710 MPI_MASTER( 1711 if (Sqw->verbose_output >= 2) 1712 printf("Isotropic_Sqw: %s: Creating Sigma(Ei=0:%g [meV]) with %li entries...(%s %s)\n", 1713 Sqw->compname, Sqw_Data->Ei_max, Sqw_Data->iqSq_length, file, Sqw->type=='c' ? "coh" : "inc"); 1714 ); 1715 Sqw_Data->Sqw = Sqw_full; /* store the S(q,w) Table (matrix) for Sqw_integrate_iqSq */ 1716 1717 /* this loop takes time: use MPI when available */ 1718 1719 for (index_w=0; index_w < Sqw_Data->iqSq_length; index_w++) { 1720 1721 /* Ei = energy of incoming neutron, typically 0:w_max or 0:2*q_max */ 1722 double Ei; /* neutron energy value in Sqw_full, up to 2*w_max */ 1723 double Ki, Vi; 1724 double Sigma=0; 1725 Ei = index_w*Sqw_Data->Ei_max/Sqw_Data->iqSq_length; 1726 Vi = sqrt(Ei/VS2E); 1727 Ki = V2K*Vi; 1728 /* sigma(Ei) = sigma/2/Ki^2 * \int q S(q,w) dq dw */ 1729 /* Eq (6) from E. Farhi et al. J. Comp. Phys. 228 (2009) 5251 */ 1730 Sigma = Ki <= 0 ? 0 : (Sqw->type=='c' ? Sqw->s_coh : Sqw->s_inc) 1731 /2/Ki/Ki * Sqw_integrate_iqSq(Sqw_Data, Ei); 1732 Table_SetElement(&iqSq, index_w, 0, Sigma ); 1733 } 1734 1735 sprintf(iqSq.filename, "[sigma/2Ki^2 int q S(q,w) dq dw] from %s", file); 1736 iqSq.min_x = 0; 1737 iqSq.max_x = Sqw_Data->Ei_max; 1738 iqSq.step_x = Sqw_Data->Ei_max/Sqw_Data->iqSq_length; 1739 iqSq.block_number = 1; 1740 1741 Sqw_Data->iqSq = iqSq; /* store the sigma(Ei) = \int q S(q,w) dq dw Table (vector) */ 1742 1743 /* (9) Compute P(w) probability =========================================== */ 1744 1745 /* set up 'density of states' */ 1746 /* uses: Sqw_full and w axes */ 1747 Sqw_Data->SW = 1748 (struct Sqw_W_struct*)calloc(w_bins, sizeof(struct Sqw_W_struct)); 1749 1750 if (!Sqw_Data->SW) { 1751 printf("Isotropic_Sqw: %s: Cannot allocate SW (%li bytes).\n" 1752 "ERROR Exiting.\n", 1753 Sqw->compname, w_bins*sizeof(struct Sqw_W_struct)); 1754 Table_Free(&Sqw_full); 1755 Table_Free(&iqSq); 1756 return(NULL); 1757 } 1758 sum = 0; 1759 for (index_w=0; index_w < w_bins ; index_w++) { 1760 double local_val = 0; 1761 double w = -w_max + index_w * w_step; 1762 for (index_q=0; index_q < q_bins ; index_q++) { /* integrate on all q values */ 1763 local_val += Table_Index(Sqw_full, index_q, index_w)*q_step*index_q*q_step; /* q*S(q,w) */ 1764 } 1765 Sqw_Data->SW[index_w].omega = w; 1766 sum += local_val; /* S(w)=\int S(q,w) dq */ 1767 /* compute cumulated probability */ 1768 Sqw_Data->SW[index_w].cumul_proba = local_val; 1769 if (index_w) Sqw_Data->SW[index_w].cumul_proba += Sqw_Data->SW[index_w-1].cumul_proba; 1770 else Sqw_Data->SW[index_w].cumul_proba = 0; 1771 } 1772 1773 if (!sum) { 1774 printf("Isotropic_Sqw: %s: Total S(q,w) intensity is NULL.\n" 1775 "ERROR Exiting.\n", Sqw->compname); 1776 Table_Free(&Sqw_full); 1777 Table_Free(&iqSq); 1778 return(NULL); 1779 } 1780 1781 /* normalize Pw distribution to 0:1 range */ 1782 for (index_w=0; index_w < w_bins ; index_w++) { 1783 Sqw_Data->SW[index_w].cumul_proba /= Sqw_Data->SW[w_bins-1].cumul_proba; 1784 } 1785 1786 if (Sqw->verbose_output > 2) 1787 printf("Isotropic_Sqw: %s: Generated normalized SW[%li] in range [0:%g]\n", 1788 Sqw->compname, w_bins, Sqw_Data->SW[w_bins-1].cumul_proba); 1789 1790 /* (10) Compute P(Q|w) probability ======================================== */ 1791 1792 /* set up Q probability table per w bin */ 1793 /* uses: Sqw_full */ 1794 Sqw_Data->SQW = 1795 (struct Sqw_Q_struct**)calloc(w_bins, sizeof(struct Sqw_Q_struct*)); 1796 1797 if (!Sqw_Data->SQW) { 1798 printf("Isotropic_Sqw: %s: Cannot allocate P(Q|w) array (%li bytes).\n" 1799 "ERROR Exiting.\n", 1800 Sqw->compname, w_bins*sizeof(struct Sqw_Q_struct*)); 1801 Table_Free(&Sqw_full); 1802 Table_Free(&iqSq); 1803 return(NULL); 1804 } 1805 for (index_w=0; index_w < w_bins ; index_w++) { 1806 Sqw_Data->SQW[index_w]= 1807 (struct Sqw_Q_struct*)calloc(q_bins, sizeof(struct Sqw_Q_struct)); 1808 1809 if (!Sqw_Data->SQW[index_w]) { 1810 printf("Isotropic_Sqw: %s: Cannot allocate P(Q|w)[%li] (%li bytes).\n" 1811 "ERROR Exiting.\n", 1812 Sqw->compname, index_w, q_bins*sizeof(struct Sqw_Q_struct)); 1813 Table_Free(&Sqw_full); 1814 Table_Free(&iqSq); 1815 return(NULL); 1816 } 1817 /* set P(Q|W) and compute total intensity */ 1818 for (index_q=0; index_q < q_bins ; index_q++) { 1819 double q = index_q * q_step; 1820 Sqw_Data->SQW[index_w][index_q].Q = q; 1821 1822 /* compute cumulated probability and take into account Jacobian with additional 'q' factor */ 1823 Sqw_Data->SQW[index_w][index_q].cumul_proba = q*Table_Index(Sqw_full, index_q, index_w); /* q*S(q,w) */ 1824 if (index_q) Sqw_Data->SQW[index_w][index_q].cumul_proba += Sqw_Data->SQW[index_w][index_q-1].cumul_proba; 1825 else Sqw_Data->SQW[index_w][index_q].cumul_proba = 0; 1826 } 1827 /* normalize P(q|w) distribution to 0:1 range */ 1828 for (index_q=0; index_q < q_bins ; 1829 Sqw_Data->SQW[index_w][index_q++].cumul_proba /= Sqw_Data->SQW[index_w][q_bins-1].cumul_proba 1830 ); 1831 1832 } 1833 if (Sqw->verbose_output > 2) 1834 printf("Isotropic_Sqw: %s: Generated P(Q|w)\n", 1835 Sqw->compname); 1836 1837 /* (11) generate quick lookup tables for SW and SQW ======================= */ 1838 1839 SW_lookup = (long*)calloc(Sqw->lookup_length, sizeof(long)); 1840 1841 if (!SW_lookup) { 1842 printf("Isotropic_Sqw: %s: Cannot allocate SW_lookup (%li bytes).\n" 1843 "Warning Will be slower.\n", 1844 Sqw->compname, Sqw->lookup_length*sizeof(long)); 1845 } else { 1846 int i; 1847 for (i=0; i < Sqw->lookup_length; i++) { 1848 double w = (double)i/(double)Sqw->lookup_length; /* a random number tabulated value */ 1849 SW_lookup[i] = Sqw_search_SW(*Sqw_Data, w); 1850 } 1851 Sqw_Data->SW_lookup = SW_lookup; 1852 } 1853 QW_lookup = (long**)calloc(w_bins, sizeof(long*)); 1854 1855 if (!QW_lookup) { 1856 printf("Isotropic_Sqw: %s: Cannot allocate QW_lookup (%li bytes).\n" 1857 "Warning Will be slower.\n", 1858 Sqw->compname, w_bins*sizeof(long*)); 1859 } else { 1860 for (index_w=0; index_w < w_bins ; index_w++) { 1861 QW_lookup[index_w] = 1862 (long*)calloc(Sqw->lookup_length, sizeof(long)); 1863 if (!QW_lookup[index_w]) { 1864 printf("Isotropic_Sqw: %s: Cannot allocate QW_lookup[%li] (%li bytes).\n" 1865 "Warning Will be slower.\n", 1866 Sqw->compname, index_w, Sqw->lookup_length*sizeof(long)); 1867 free(QW_lookup); Sqw_Data->QW_lookup = QW_lookup = NULL; break; 1868 } else { 1869 int i; 1870 for (i=0; i < Sqw->lookup_length; i++) { 1871 double w = (double)i/(double)Sqw->lookup_length; /* a random number tabulated value */ 1872 QW_lookup[index_w][i] = Sqw_search_Q_proba_per_w(*Sqw_Data, w, index_w); 1873 } 1874 } 1875 } 1876 Sqw_Data->QW_lookup = QW_lookup; 1877 } 1878 if ((Sqw_Data->QW_lookup || Sqw_Data->SW_lookup) && Sqw->verbose_output > 2) 1879 printf("Isotropic_Sqw: %s: Generated lookup tables with %li entries\n", 1880 Sqw->compname, Sqw->lookup_length); 1881 1882 return(Sqw_Data); 1883} /* end Sqw_readfile */ 1884 1885/***************************************************************************** 1886* Sqw_init_read: Read coherent/incoherent Sqw data files 1887* Returns Sqw total intensity, or 0 (error) 1888* Used in : INITIALIZE (1) 1889*****************************************************************************/ 1890double Sqw_init(struct Sqw_sample_struct *Sqw, char *file_coh, char *file_inc) 1891{ 1892 double ret=0; 1893 1894 /* read files */ 1895 struct Sqw_Data_struct *d_inc, *d_coh; 1896 Sqw->type = 'i'; 1897 d_inc = Sqw_readfile(Sqw, file_inc, &(Sqw->Data_inc)); 1898 Sqw->type = 'c'; 1899 d_coh = Sqw_readfile(Sqw, file_coh, &(Sqw->Data_coh)); 1900 1901 if (d_inc && !d_inc->intensity && Sqw->s_inc>0) { 1902 MPI_MASTER( 1903 if (Sqw->verbose_output > 0) 1904 printf("Isotropic_Sqw: %s: Using Isotropic elastic incoherent scattering (sigma=%g [barns])\n", Sqw->compname, Sqw->s_inc); 1905 ); 1906 ret=1; 1907 } 1908 1909 if (!d_inc || !d_coh) return(0); 1910 1911 d_coh->type = 'c'; 1912 Sqw->Data_inc.type = 'i'; 1913 if (d_coh && !d_coh->intensity && Sqw->s_coh) 1914 printf("Isotropic_Sqw: %s: Coherent scattering Sqw intensity is null.\n" 1915 "Warning Disabling coherent scattering.\n", Sqw->compname); 1916 if (d_inc && d_coh && d_inc->intensity && d_coh->intensity) { 1917 char msg[80]; 1918 strcpy(msg, ""); 1919 /* check dimensions/limits for Q, Energy in coh and inc Tables */ 1920 if (d_inc->q_bins != d_coh->q_bins) 1921 strcpy(msg, "Q axis size"); 1922 if (d_inc->w_bins != d_coh->w_bins) 1923 strcpy(msg, "Energy axis size"); 1924 if (d_inc->q_max != d_coh->q_max) 1925 strcpy(msg, "Q axis limits"); 1926 if (d_inc->w_max != d_coh->w_max) 1927 strcpy(msg, "Energy axis limits"); 1928 if (strlen(msg)) { 1929 printf("Isotropic_Sqw: %s: Sqw data from files %s and %s do not match\n" 1930 "WARNING wrong %s\n", 1931 Sqw->compname, file_coh, file_inc, msg); 1932 } 1933 } 1934 1935 if (!ret) ret=d_inc->intensity+d_coh->intensity; 1936 return(ret); 1937} /* Sqw_init */ 1938 1939#endif /* definied ISOTROPIC_SQW */ 1940%} 1941 1942 1943/*****************************************************************************/ 1944/*****************************************************************************/ 1945 1946DECLARE 1947%{ 1948 struct Sqw_sample_struct VarSqw; 1949 int *columns; 1950 off_struct offdata; 1951%} 1952 1953/*****************************************************************************/ 1954/*****************************************************************************/ 1955 1956INITIALIZE 1957%{ 1958 int i; 1959 /* check for parameters */ 1960 columns = (int[])powder_format; 1961 1962 VarSqw.verbose_output= verbose; 1963 VarSqw.shape = -1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape */ 1964 if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) { 1965 if (off_init(geometry, xwidth, yheight, zdepth, 0, &offdata)) { 1966 VarSqw.shape=3; thickness=0; concentric=0; 1967 } 1968 } 1969 else if (xwidth && yheight && zdepth) VarSqw.shape=1; /* box */ 1970 else if (radius > 0 && yheight) VarSqw.shape=0; /* cylinder */ 1971 else if (radius > 0 && !yheight) VarSqw.shape=2; /* sphere */ 1972 1973 if (VarSqw.shape < 0) 1974 exit(fprintf(stderr,"Isotropic_Sqw: %s: sample has invalid dimensions.\n" 1975 "ERROR Please check parameter values (xwidth, yheight, zdepth, radius).\n", NAME_CURRENT_COMP)); 1976 1977 1978 1979 if (thickness) { 1980 if (radius && (radius < fabs(thickness) )) { 1981 fprintf(stderr,"Isotropic_Sqw: %s: hollow sample thickness is larger than its volume (sphere/cylinder).\n" 1982 "WARNING Please check parameter values. Using bulk sample (thickness=0).\n", NAME_CURRENT_COMP); 1983 thickness=0; 1984 } 1985 else if (!radius && (xwidth < 2*fabs(thickness) || yheight < 2*fabs(thickness) || zdepth < 2*fabs(thickness))) { 1986 fprintf(stderr,"Isotropic_Sqw: %s: hollow sample thickness is larger than its volume (box).\n" 1987 "WARNING Please check parameter values.\n", NAME_CURRENT_COMP); 1988 } 1989 } 1990 MPI_MASTER( 1991 if (VarSqw.verbose_output) { 1992 switch (VarSqw.shape) { 1993 case 0: printf("Isotropic_Sqw: %s: is a %scylinder: radius=%f thickness=%f height=%f [J Comp Phys 228 (2009) 5251]\n", 1994 NAME_CURRENT_COMP, (thickness ? "hollow " : ""), 1995 radius,fabs(thickness),yheight); 1996 break; 1997 case 1: printf("Isotropic_Sqw: %s: is a %sbox: width=%f height=%f depth=%f \n", 1998 NAME_CURRENT_COMP, (thickness ? "hollow " : ""), xwidth,yheight,zdepth); 1999 break; 2000 case 2: printf("Isotropic_Sqw: %s: is a %ssphere: radius=%f thickness=%f\n", 2001 NAME_CURRENT_COMP, (thickness ? "hollow " : ""), 2002 radius,fabs(thickness)); 2003 break; 2004 case 3: printf("Isotropic_Sqw: %s: is a volume defined from file %s\n", 2005 NAME_CURRENT_COMP, geometry); 2006 } 2007 } 2008 ); 2009 2010 if (concentric && !thickness) { 2011 printf("Isotropic_Sqw: %s:Can not use concentric mode\n" 2012 "WARNING on non hollow shape. Ignoring.\n", 2013 NAME_CURRENT_COMP); 2014 concentric=0; 2015 } 2016 2017 strncpy(VarSqw.compname, NAME_CURRENT_COMP, 256); 2018 VarSqw.T2E =(1/11.605); /* Kelvin to meV = 1000*KB/e */ 2019 VarSqw.sqSE2K = (V2K*SE2V)*(V2K*SE2V); 2020 VarSqw.sqw_threshold = (threshold > 0 ? threshold : 0); 2021 VarSqw.s_abs = sigma_abs; 2022 VarSqw.s_coh = sigma_coh; 2023 VarSqw.s_inc = sigma_inc; /* s_scatt member initialized in Sqw_init */ 2024 VarSqw.maxloop = 100; /* atempts to close triangle */ 2025 VarSqw.minevents = 100; /* minimal # of events required to get dynamical range */ 2026 VarSqw.neutron_removed = 0; 2027 VarSqw.neutron_enter = 0; 2028 VarSqw.neutron_pmult = 0; 2029 VarSqw.neutron_exit = 0; 2030 VarSqw.mat_rho = rho; 2031 VarSqw.sqw_norm = norm; 2032 VarSqw.mean_scatt= 0; 2033 VarSqw.mean_abs = 0; 2034 VarSqw.psum_scatt= 0; 2035 VarSqw.single_coh= 0; 2036 VarSqw.single_inc= 0; 2037 VarSqw.multi = 0; 2038 VarSqw.barns = powder_barns; 2039 VarSqw.sqw_classical = classical; 2040 VarSqw.lookup_length=100; 2041 VarSqw.mat_weight = weight; 2042 VarSqw.mat_density = density; 2043 if (quantum_correction && strlen(quantum_correction)) 2044 strncpy(VarSqw.Q_correction, quantum_correction, 256); 2045 else 2046 strncpy(VarSqw.Q_correction, "default", 256); 2047 2048 /* PowderN compatibility members */ 2049 VarSqw.Dd = powder_Dd; 2050 VarSqw.DWfactor = powder_DW; 2051 VarSqw.Temperature= T; 2052 for (i=0; i< 9; i++) VarSqw.column_order[i] = columns[i]; 2053 VarSqw.column_order[8] = (VarSqw.column_order[0] >= 0 ? 0 : 2); 2054 2055 /* optional ways to define rho */ 2056 if (!VarSqw.mat_rho && powder_Vc > 0) 2057 VarSqw.mat_rho = 1/powder_Vc; 2058 /* import the data files ================================================== */ 2059 if (!Sqw_init(&VarSqw, Sqw_coh, Sqw_inc)) { 2060 printf("Isotropic_Sqw: %s: ERROR importing data files (Sqw_init coh=%s inc=%s).\n",NAME_CURRENT_COMP, Sqw_coh, Sqw_inc); 2061 } 2062 if ( VarSqw.s_coh < 0) VarSqw.s_coh=0; 2063 if ( VarSqw.s_inc < 0) VarSqw.s_inc=0; 2064 if ( VarSqw.s_abs < 0) VarSqw.s_abs=0; 2065 if ((VarSqw.s_coh > 0 || VarSqw.s_inc > 0) && VarSqw.mat_rho <= 0) { 2066 printf("Isotropic_Sqw: %s: WARNING: Null density (V_rho). Unactivating component.\n",NAME_CURRENT_COMP); 2067 VarSqw.s_coh=VarSqw.s_inc=0; 2068 } 2069 /* 100: convert from barns to fm^2 */ 2070 VarSqw.my_a_v =(VarSqw.mat_rho*100*VarSqw.s_abs*2200); 2071 VarSqw.my_s =(VarSqw.mat_rho*100*(VarSqw.s_coh>0 ? VarSqw.s_coh : 0 2072 +VarSqw.s_inc>0 ? VarSqw.s_inc : 0)); 2073 MPI_MASTER( 2074 if ((VarSqw.s_coh > 0 || VarSqw.s_inc > 0) && !VarSqw.Temperature 2075 && (VarSqw.Data_coh.intensity || VarSqw.Data_inc.intensity) 2076 && VarSqw.verbose_output) 2077 printf("Isotropic_Sqw: %s: Sample temperature not defined (T=0).\n" 2078 "Warning Disabling detailed balance.\n", NAME_CURRENT_COMP); 2079 if (VarSqw.s_coh<=0 && VarSqw.s_inc<=0) { 2080 printf("Isotropic_Sqw: %s: Scattering cross section is zero\n" 2081 "ERROR (sigma_coh, sigma_inc).\n",NAME_CURRENT_COMP); 2082 } 2083 ); 2084 if (d_phi) d_phi = fabs(d_phi)*DEG2RAD; 2085 2086 if (d_phi > PI) d_phi = 0; /* V_scatt on 4*PI */ 2087 2088 if (d_phi && order != 1) { 2089 printf("Isotropic_Sqw: %s: Focusing can only apply for single\n" 2090 " scattering. Setting to order=1.\n", 2091 NAME_CURRENT_COMP); 2092 order = 1; 2093 } 2094 2095 /* request statistics */ 2096 if (VarSqw.verbose_output > 1) { 2097 Sqw_diagnosis(&VarSqw, &VarSqw.Data_coh); 2098 Sqw_diagnosis(&VarSqw, &VarSqw.Data_inc); 2099 } 2100 2101 for (i=0; i < 2; i++) { 2102 struct Sqw_Data_struct Data_sqw; 2103 Data_sqw = (i == 0 ? VarSqw.Data_coh : VarSqw.Data_inc); 2104 Table_Free(&(Data_sqw.Sqw)); 2105 } 2106 2107/* end INITIALIZE */ 2108%} 2109 2110/*****************************************************************************/ 2111/*****************************************************************************/ 2112TRACE 2113%{ 2114 2115int intersect=0; /* flag to continue/stop */ 2116double t0, t1, t2, t3; /* times for intersections */ 2117double dt0, dt1, dt2, dt; /* time intervals */ 2118double k=0, Ei=0; 2119double v=0, vf=0; 2120double d_path; /* total path length for straight trajectory */ 2121double my_a; /* absorption cross-section scaled to velocity (2200) */ 2122double ws, p_scatt; /* probability for scattering/absorption and for */ 2123 /* interaction along d_path */ 2124double tmp_rand; /* temporary var */ 2125double ratio_w=0, ratio_q=0; /* variables for bilinear interpolation */ 2126double q11, q21, q22, q12; 2127double omega=0; /* energy transfer */ 2128double q=0; /* wavevector transfer */ 2129long index_w; /* energy index for table look-up SW */ 2130long index_q; /* Q index for table look-up P(Q|w) */ 2131double theta=0, costheta=0; /* for the choice of kf direction */ 2132double u1x,u1y,u1z; 2133double u2x,u2y,u2z; 2134double u0x,u0y,u0z; 2135int index_counter; 2136int flag=0; 2137int flag_concentric=0; 2138int flag_ishollow=0; 2139double solid_angle=0; 2140double my_t=0; 2141double p_mult=1; 2142double mc_trans, p_trans, mc_scatt; 2143double coh=0, inc=0; 2144struct Sqw_Data_struct Data_sqw; 2145 2146 2147/* Store Initial neutron state */ 2148 2149VarSqw.ki_x = V2K*vx; 2150VarSqw.ki_y = V2K*vy; 2151VarSqw.ki_z = V2K*vz; 2152VarSqw.ti = t; 2153VarSqw.vi = 0; 2154VarSqw.ki = 0; 2155VarSqw.type = '\0'; 2156 2157do { /* Main interaction loop. Ends with intersect=0 */ 2158 2159 /* Intersection neutron trajectory / sample (sample surface) */ 2160 if (VarSqw.s_coh > 0 || VarSqw.s_inc > 0) { 2161 if (thickness >= 0) { 2162 if (VarSqw.shape==0) 2163 intersect=cylinder_intersect(&t0,&t3, x,y,z,vx,vy,vz, radius,yheight); 2164 else if (VarSqw.shape==1) 2165 intersect=box_intersect (&t0,&t3, x,y,z,vx,vy,vz, xwidth,yheight,zdepth); 2166 else if (VarSqw.shape==2) 2167 intersect=sphere_intersect (&t0,&t3, x,y,z,vx,vy,vz, radius); 2168 else if (VarSqw.shape == 3) 2169 intersect=off_intersect(&t0, &t3, NULL, NULL, x, y, z, vx, vy, vz, offdata ); 2170 } else { 2171 if (VarSqw.shape==0) 2172 intersect=cylinder_intersect(&t0,&t3, x,y,z,vx,vy,vz, radius-thickness, 2173 yheight-2*thickness > 0 ? yheight-2*thickness : yheight); 2174 else if (VarSqw.shape==1) 2175 intersect=box_intersect (&t0,&t3, x,y,z,vx,vy,vz, 2176 xwidth-2*thickness > 0 ? xwidth-2*thickness : xwidth, 2177 yheight-2*thickness > 0 ? yheight-2*thickness : yheight, 2178 zdepth-2*thickness > 0 ? zdepth-2*thickness : zdepth); 2179 else if (VarSqw.shape==2) 2180 intersect=sphere_intersect (&t0,&t3, x,y,z,vx,vy,vz, radius-thickness); 2181 else if (VarSqw.shape == 3) 2182 intersect=off_intersect(&t0, &t3, NULL, NULL, x, y, z, vx, vy, vz, offdata ); 2183 } 2184 } else intersect=0; 2185 2186 /* Computing the intermediate times */ 2187 if (intersect) { 2188 flag_ishollow = 0; 2189 if (thickness > 0) { 2190 if (VarSqw.shape==0 && cylinder_intersect(&t1,&t2, x,y,z,vx,vy,vz, radius-thickness, 2191 yheight-2*thickness > 0 ? yheight-2*thickness : yheight)) 2192 flag_ishollow=1; 2193 else if (VarSqw.shape==2 && sphere_intersect (&t1,&t2, x,y,z,vx,vy,vz, radius-thickness)) 2194 flag_ishollow=1; 2195 else if (VarSqw.shape==1 && box_intersect(&t1,&t2, x,y,z,vx,vy,vz, 2196 xwidth-2*thickness > 0 ? xwidth-2*thickness : xwidth, 2197 yheight-2*thickness > 0 ? yheight-2*thickness : yheight, 2198 zdepth-2*thickness > 0 ? zdepth-2*thickness : zdepth)) 2199 flag_ishollow=1; 2200 } else if (thickness<0) { 2201 if (VarSqw.shape==0 && cylinder_intersect(&t1,&t2, x,y,z,vx,vy,vz, radius,yheight)) 2202 flag_ishollow=1; 2203 else if (VarSqw.shape==2 && sphere_intersect (&t1,&t2, x,y,z,vx,vy,vz, radius)) 2204 flag_ishollow=1; 2205 else if (VarSqw.shape==1 && box_intersect(&t1,&t2, x,y,z,vx,vy,vz, xwidth, yheight, zdepth)) 2206 flag_ishollow=1; 2207 } 2208 if (!flag_ishollow) t1 = t2 = t3; /* no empty space inside */ 2209 } else break; /* neutron does not hit sample: transmitted */ 2210 2211 if (intersect) { /* the neutron hits the sample */ 2212 2213 if (t0 > 0) { /* we are before the sample */ 2214 PROP_DT(t0); /* propagates neutron to the entry of the sample */ 2215 } else if (t1 > 0 && t1 > t0) { /* we are inside first part of the sample */ 2216 /* no propagation, stay inside */ 2217 } else if (t2 > 0 && t2 > t1) { /* we are in the hole */ 2218 PROP_DT(t2); /* propagate to inner surface of 2nd part of sample */ 2219 } else if (t3 > 0 && t3 > t2) { /* we are in the 2nd part of sample */ 2220 /* no propagation, stay inside */ 2221 } 2222 2223 dt0=t1-(t0 > 0 ? t0 : 0); /* Time in first part of hollow/cylinder/box */ 2224 dt1=t2-(t1 > 0 ? t1 : 0); /* Time in hole */ 2225 dt2=t3-(t2 > 0 ? t2 : 0); /* Time in 2nd part of hollow cylinder */ 2226 2227 if (dt0 < 0) dt0 = 0; 2228 if (dt1 < 0) dt1 = 0; 2229 if (dt2 < 0) dt2 = 0; 2230 2231 /* initialize concentric mode */ 2232 if (concentric && !flag_concentric && t0 >= 0 2233 && VarSqw.shape==0 && thickness) { 2234 flag_concentric=1; 2235 } 2236 2237 if (flag_concentric == 1) { 2238 dt1=dt2=0; /* force exit when reaching hole/2nd part */ 2239 } 2240 2241 if (!dt0 && !dt2) { 2242 intersect = 0; /* the sample was passed entirely */ 2243 break; 2244 } 2245 2246 VarSqw.neutron_enter++; 2247 p_mult = 1; 2248 if (!v) { 2249 v = vx*vx+vy*vy+vz*vz; 2250 v = sqrt(v); 2251 } 2252 k = V2K*v; 2253 Ei = VS2E*v*v; 2254 2255 if (!VarSqw.vi) VarSqw.vi = v; 2256 if (!VarSqw.ki) VarSqw.ki = k; 2257 2258 if (v <= 0) { 2259 printf("Isotropic_Sqw: %s: ERROR: Null velocity !\n",NAME_CURRENT_COMP); 2260 VarSqw.neutron_removed++; 2261 ABSORB; /* should never occur */ 2262 } 2263 2264 /* check for scattering event */ 2265 my_a = VarSqw.my_a_v / v; /* absorption 'mu' */ 2266 /* compute total scattering X section */ 2267 /* \int q S(q) dq /2 /ki^2 sigma OR bare Xsection*/ 2268 /* contains the 4*PI*kf/ki factor */ 2269 coh = VarSqw.s_coh; 2270 inc = VarSqw.s_inc; 2271 if (k && VarSqw.s_coh>0 && VarSqw.Data_coh.intensity) { 2272 double Ei = VS2E*v*v; 2273 double index_Ei = Ei / (VarSqw.Data_coh.Ei_max/VarSqw.Data_coh.iqSq_length); 2274 coh = Table_Value2d(VarSqw.Data_coh.iqSq, index_Ei, 0); 2275 } 2276 if (k && VarSqw.s_inc>0 && VarSqw.Data_inc.intensity) { 2277 double Ei = VS2E*v*v; 2278 double index_Ei = Ei / (VarSqw.Data_inc.Ei_max/VarSqw.Data_inc.iqSq_length); 2279 inc = Table_Value2d(VarSqw.Data_inc.iqSq, index_Ei, 0); 2280 } 2281 if (coh<0) coh=0; 2282 if (inc<0) inc=0; 2283 VarSqw.my_s =(VarSqw.mat_rho*100*(coh + inc)); 2284 2285 my_t = my_a + VarSqw.my_s; /* total scattering Xsect */ 2286 if (my_t <= 0) { 2287 if (VarSqw.neutron_removed<VarSqw.maxloop) printf("Isotropic_Sqw: %s: ERROR: Null total cross section %g. Removing event.\n", 2288 NAME_CURRENT_COMP, my_t); 2289 VarSqw.neutron_removed++; 2290 ABSORB; /* should never occur */ 2291 } else if (VarSqw.my_s <= 0) { 2292 if (VarSqw.verbose_output > 1 && VarSqw.neutron_removed<VarSqw.maxloop) 2293 printf("Isotropic_Sqw: %s: Warning: Null scattering cross section %g. Ignoring.\n", 2294 NAME_CURRENT_COMP, VarSqw.my_s); 2295 VarSqw.my_s = 0; 2296 } 2297 2298 /* Proba of scattering vs absorption (integrating along the whole trajectory) */ 2299 ws = VarSqw.my_s/my_t; /* (inc+coh)/(inc+coh+abs) */ 2300 d_path = v*( dt0 +dt2 ); /* total path lenght in sample */ 2301 /* Proba of transmission/interaction along length d_path */ 2302 p_trans = exp(-my_t*d_path); 2303 p_scatt = 1 - p_trans; /* portion of beam which scatters */ 2304 2305 flag = 0; /* flag used for propagation to exit point before ending */ 2306 2307 /* are we next to the exit ? probably no scattering (avoid rounding errors) */ 2308 if (VarSqw.my_s*d_path <= 4e-7) { 2309 flag = 1; /* No interaction before the exit */ 2310 } 2311 /* force a given fraction of the beam to scatter */ 2312 if (p_interact>0 && p_interact<=1) { 2313 /* we force a portion of the beam to interact */ 2314 /* This is used to improve statistics on single scattering (and multiple) */ 2315 if (!SCATTERED) mc_trans = 1-p_interact; 2316 else mc_trans = 1-p_interact/(4*SCATTERED+1); /* reduce effect on multi scatt */ 2317 } else { 2318 mc_trans = p_trans; /* 1 - p_scatt */ 2319 } 2320 mc_scatt = 1 - mc_trans; /* portion of beam to scatter (or force to) */ 2321 if (mc_scatt <= 0 || mc_scatt>1) flag=1; 2322 /* MC choice: Interaction or transmission ? */ 2323 if (!flag && mc_scatt > 0 && (mc_scatt >= 1 || rand01() < mc_scatt)) { /* Interaction neutron/sample */ 2324 p_mult *= ws; /* Update weight ; account for absorption and retain scattered fraction */ 2325 /* we have chosen portion mc_scatt of beam instead of p_scatt, so we compensate */ 2326 if (!mc_scatt) ABSORB; 2327 p_mult *= fabs(p_scatt/mc_scatt); /* lower than 1 */ 2328 } else { 2329 flag = 1; /* Transmission : no interaction neutron/sample */ 2330 if (!VarSqw.type) VarSqw.type = 't'; 2331 if (!mc_trans) ABSORB; 2332 p_mult *= fabs(p_trans/mc_trans); /* attenuate beam by portion which is scattered (and left along) */ 2333 } 2334 2335 if (flag) { /* propagate to exit of sample and finish */ 2336 intersect = 0; 2337 p *= p_mult; /* apply absorption correction */ 2338 PROP_DT(dt0+dt2); 2339 break; /* exit main multi scatt while loop */ 2340 } 2341 } /* end if intersect the neutron hits the sample */ 2342 else break; 2343 2344 if (intersect) { /* scattering event */ 2345 double kf=0, kf1, kf2; 2346 /* mean scattering probability and absorption fraction */ 2347 VarSqw.mean_scatt += (1-exp(-VarSqw.my_s*d_path))*p; 2348 VarSqw.mean_abs += (1-ws)*p; 2349 VarSqw.psum_scatt += p; 2350 2351 /* Decaying exponential distribution of the path length before scattering */ 2352 /* Select a point at which to scatter the neutron, taking 2353 secondary extinction into account. */ 2354 if (my_t*d_path < 1e-6) 2355 /* For very weak scattering, use simple uniform sampling of scattering 2356 point to avoid rounding errors. */ 2357 dt = rand0max(d_path); /* length */ 2358 else 2359 dt = -log(1 - rand0max((1 - exp(-my_t*d_path)))) / my_t; /* length */ 2360 dt /= v; /* Time from present position to scattering point */ 2361 2362 /* If t0 is in hole, propagate to next part of the hollow cylinder */ 2363 if (dt1 > 0 && dt0 > 0 && dt > dt0) dt += dt1; 2364 2365 /* Neutron propagation to the scattering point */ 2366 PROP_DT(dt); 2367 2368 /* choice between coherent/incoherent scattering */ 2369 tmp_rand = rand01(); 2370 /* local description at the scattering point (scat probability for atom) */ 2371 tmp_rand *= (coh+inc); 2372 2373 flag=0; 2374 if (VarSqw.s_inc>0 && tmp_rand < inc) { 2375 /* CASE 1: incoherent case */ 2376 if (!VarSqw.Data_inc.intensity) { 2377 /* CASE 1a: no incoherent Sqw from file, use isotropic V-like */ 2378 if (d_phi && order == 1) { 2379 randvec_target_rect_angular(&u1x, &u1y, &u1z, &solid_angle, 2380 vx, vy, vz, 2*PI, d_phi, ROT_A_CURRENT_COMP); 2381 p_mult *= solid_angle/4/PI; /* weighted by focused range to total range */ 2382 } else 2383 randvec_target_circle(&u1x, &u1y, &u1z, NULL, vx, vy, vz, 0); 2384 2385 vx = u1x; vy = u1y; vz = u1z; 2386 vf = v; kf = k; 2387 if (!VarSqw.type) VarSqw.type = 'v'; 2388 SCATTER; 2389 } else { 2390 /* CASE 1b: incoherent Sqw from file */ 2391 if (VarSqw.Data_inc.intensity) { 2392 Data_sqw = VarSqw.Data_inc; 2393 if (!VarSqw.type) VarSqw.type = 'i'; 2394 flag = 1; 2395 } 2396 } 2397 } else if (VarSqw.s_coh>0 && tmp_rand > VarSqw.s_inc) { 2398 if (VarSqw.Data_coh.intensity) { 2399 /* CASE2: coherent case */ 2400 Data_sqw = VarSqw.Data_coh; 2401 if (!VarSqw.type) VarSqw.type = 'c'; 2402 flag = 1; 2403 } 2404 } 2405 2406 if (flag) { /* true when S(q,w) table exists (Data_sqw) */ 2407 2408 double alpha=0, alpha0; 2409 /* give us a limited number of tries for scattering: choose W then Q */ 2410 for (index_counter=VarSqw.maxloop; index_counter > 0 ; index_counter--) { 2411 double randmax; 2412 /* MC choice: energy transfer w=Ei-Ef in the S(w) = SW */ 2413 omega = 0; 2414 /* limit energy transfer range in -w_max:Ei */ 2415 index_w = (long)floor((1+Ei/Data_sqw.w_max)/2*Data_sqw.w_bins); 2416 if (index_w >= Data_sqw.w_bins) index_w = Data_sqw.w_bins-1; 2417 randmax = Data_sqw.SW[index_w].cumul_proba; 2418 tmp_rand = rand0max(randmax < 1 ? randmax : 1); 2419 2420 /* energy index for rand > cumul SW */ 2421 index_w = Sqw_search_SW(Data_sqw, tmp_rand); 2422 VarSqw.rw = (double)index_w; 2423 if (index_w >= 0 && &(Data_sqw.SW[index_w]) != NULL) { 2424 if (Data_sqw.w_bins > 1) { 2425 double w1, w2; 2426 if (index_w > 0) { /* interpolate linearly energy */ 2427 ratio_w = (tmp_rand - Data_sqw.SW[index_w-1].cumul_proba) 2428 /(Data_sqw.SW[index_w].cumul_proba - Data_sqw.SW[index_w-1].cumul_proba); 2429 /* ratio_w=0 omega[index_w-1], ratio=1 omega[index] */ 2430 w1 = Data_sqw.SW[index_w-1].omega; w2 = Data_sqw.SW[index_w].omega; 2431 } else { /* index_w = 0 interpolate to 0 energy */ 2432 /* ratio_w=0 omega=0, ratio=1 omega[index] */ 2433 w1 = Data_sqw.SW[index_w].omega; w2= Data_sqw.SW[index_w+1].omega; 2434 if (!w2 && index_w+1 < Data_sqw.w_bins) 2435 w2= Data_sqw.SW[index_w+1].omega; 2436 if (Data_sqw.w_bins && Data_sqw.SW[index_w].cumul_proba) { 2437 ratio_w = tmp_rand/Data_sqw.SW[index_w].cumul_proba; 2438 } else ratio_w=0; 2439 } 2440 if (ratio_w<0) ratio_w=0; else if (ratio_w>1) ratio_w=1; 2441 omega = (1-ratio_w)*w1 + ratio_w*w2; 2442 } else { 2443 ratio_w = 0; 2444 omega = Data_sqw.SW[index_w].omega; 2445 } 2446 } else { 2447 if (VarSqw.verbose_output >= 3 && VarSqw.neutron_removed<VarSqw.maxloop) 2448 printf("Isotropic_Sqw: %s: Warning: No suitable w transfer for index_w=%li.\n", 2449 NAME_CURRENT_COMP, index_w); 2450 continue; /* no W value: try again with an other energy transfer */ 2451 } 2452 2453 /* MC choice: momentum transfer Q in P(Q|w) */ 2454 /* limit momentum choice to 0:Q1=SE2V*V2K*(sqrt(Ei)+sqrt(Ei+w)) ? */ 2455 if (omega <= Ei) 2456 index_q = (long)floor( 2457 (SE2V*V2K*(sqrt(Ei)+sqrt(Ei-omega)) 2458 /Data_sqw.q_max)*Data_sqw.q_bins); 2459 else 2460 index_q = (long)floor(2*k/Data_sqw.q_max*Data_sqw.q_bins); 2461 2462 randmax = Data_sqw.SQW[index_w][index_q].cumul_proba; 2463 tmp_rand = rand0max(randmax < 1 ? randmax : 1); 2464 2465 /* momentum index for rand > cumul SQ|W */ 2466 index_q = Sqw_search_Q_proba_per_w(Data_sqw, tmp_rand, index_w); 2467 VarSqw.rq = (double)index_q; 2468 2469 if (index_q >= 0 && &(Data_sqw.SQW[index_w]) != NULL) { 2470 if (Data_sqw.q_bins > 1 && index_q > 0) { 2471 if (index_w > 0 && Data_sqw.w_bins > 1) { 2472 /* bilinear interpolation on - side: index_w > 0, index_q > 0 */ 2473 ratio_q = (tmp_rand - Data_sqw.SQW[index_w][index_q-1].cumul_proba) 2474 /(Data_sqw.SQW[index_w][index_q].cumul_proba 2475 - Data_sqw.SQW[index_w][index_q-1].cumul_proba); 2476 q22 = Data_sqw.SQW[index_w] [index_q].Q; 2477 q11 = Data_sqw.SQW[index_w-1][index_q-1].Q; 2478 q21 = Data_sqw.SQW[index_w] [index_q-1].Q; 2479 q12 = Data_sqw.SQW[index_w-1][index_q].Q; 2480 if (ratio_q<0) ratio_q=0; else if (ratio_q>1) ratio_q=1; 2481 q = (1-ratio_w)*(1-ratio_q)*q11+ratio_w*(1-ratio_q)*q21 2482 + ratio_w*ratio_q*q22 +(1-ratio_w)*ratio_q*q12; 2483 } else { /* bilinear interpolation on + side: index_w=0, index_q > 0 */ 2484 ratio_q = (tmp_rand - Data_sqw.SQW[index_w][index_q-1].cumul_proba) 2485 /(Data_sqw.SQW[index_w][index_q].cumul_proba 2486 - Data_sqw.SQW[index_w][index_q-1].cumul_proba); 2487 q11 = Data_sqw.SQW[index_w] [index_q-1].Q; 2488 q12 = Data_sqw.SQW[index_w] [index_q].Q; 2489 if (ratio_q<0) ratio_q=0; else if (ratio_q>1) ratio_q=1; 2490 if (index_w < Data_sqw.w_bins-1 && Data_sqw.w_bins > 1) { 2491 q22 = Data_sqw.SQW[index_w+1][index_q].Q; 2492 q21 = Data_sqw.SQW[index_w+1][index_q-1].Q; 2493 q = (1-ratio_w)*(1-ratio_q)*q11+ratio_w*(1-ratio_q)*q21 2494 + ratio_w*ratio_q*q22 +(1-ratio_w)*ratio_q*q12; 2495 } else { 2496 q = (1-ratio_q)*q11 + ratio_q*q12; 2497 } 2498 } 2499 } else { 2500 q = Data_sqw.SQW[index_w][index_q].Q; 2501 } 2502 } else { 2503 if (VarSqw.verbose_output >= 3 && VarSqw.neutron_removed<VarSqw.maxloop) 2504 printf("Isotropic_Sqw: %s: Warning: No suitable q transfer\n" 2505 " for w=%g\n", 2506 NAME_CURRENT_COMP, omega); 2507 VarSqw.neutron_removed++; 2508 continue; /* no Q value for this w choice */ 2509 } 2510 2511 /* Search for length of final wave vector kf */ 2512 /* kf is such that : hbar*w = hbar*hbar/2/m*(k*k - kf*kf) */ 2513 /* acceptable values for kf are kf1 and kf2 */ 2514 if (!solve_2nd_order(&kf1, &kf2, 1, 0, -k*k + VarSqw.sqSE2K*omega)) { 2515 if (VarSqw.verbose_output >= 3 && VarSqw.neutron_removed<VarSqw.maxloop) 2516 printf("Isotropic_Sqw: %s: Warning: imaginary root for w=%g q=%g Ei=%g (triangle can not close)\n", 2517 NAME_CURRENT_COMP, omega, q, Ei); 2518 VarSqw.neutron_removed++; 2519 continue; /* all roots are imaginary */ 2520 } 2521 2522 /* kf1 and kf2 are opposite */ 2523 kf = fabs(kf1); 2524 vf = K2V*kf; 2525 2526 /* Search of the direction of kf such that : q = ki - kf */ 2527 /* cos theta = (ki2+kf2-q2)/(2ki kf) */ 2528 2529 costheta= (k*k+kf*kf-q*q)/(2*kf*k); /* this is cos(theta) */ 2530 2531 if (-1 < costheta && costheta < 1) { 2532 break; /* satisfies q momentum conservation */ 2533 } 2534/* else continue; */ 2535 2536 /* exit for loop on success */ 2537 } /* end for index_counter */ 2538 2539 if (!index_counter) { /* for loop ended: failure for scattering */ 2540 intersect=0; /* Could not scatter: finish multiple scattering loop */ 2541 if (VarSqw.verbose_output >= 2 && VarSqw.neutron_removed<VarSqw.maxloop) 2542 printf("Isotropic_Sqw: %s: Warning: No scattering [q,w] conditions\n" 2543 " last try (%i): type=%c w=%g q=%g cos(theta)=%g k=%g\n", 2544 NAME_CURRENT_COMP, VarSqw.maxloop, (VarSqw.type ? VarSqw.type : '-'), omega, q, costheta, k); 2545 VarSqw.neutron_removed++; 2546 if (order && SCATTERED != order) ABSORB; 2547 break; /* finish multiple scattering loop */ 2548 } 2549 2550 /* scattering angle from ki to DS cone */ 2551 theta = acos(costheta); 2552 2553 /* Choose point on Debye-Scherrer cone */ 2554 if (order == 1 && d_phi) 2555 { /* relate height of detector to the height on DS cone */ 2556 double cone_focus; 2557 cone_focus = sin(d_phi/2)/sin(theta); 2558 /* If full Debye-Scherrer cone is within d_phi, don't focus */ 2559 if (cone_focus < -1 || cone_focus > 1) d_phi = 0; 2560 /* Otherwise, determine alpha to rotate from scattering plane 2561 into d_phi focusing area*/ 2562 else alpha = 2*asin(cone_focus); 2563 if (d_phi) p_mult *= alpha/PI; 2564 } 2565 if (d_phi) { 2566 /* Focusing */ 2567 alpha = fabs(alpha); 2568 /* Trick to get scattering for pos/neg theta's */ 2569 alpha0= 2*rand01()*alpha; 2570 if (alpha0 > alpha) { 2571 alpha0=PI+(alpha0-1.5*alpha); 2572 } else { 2573 alpha0=alpha0-0.5*alpha; 2574 } 2575 } 2576 else 2577 alpha0 = PI*randpm1(); 2578 2579 /* now find a nearly vertical rotation axis (u1) : 2580 * Either 2581 * (v along Z) x (X axis) -> nearly Y axis 2582 * Or 2583 * (v along X) x (Z axis) -> nearly Y axis 2584 */ 2585 if (fabs(scalar_prod(1,0,0,vx/v,vy/v,vz/v)) < fabs(scalar_prod(0,0,1,vx/v,vy/v,vz/v))) { 2586 u1x = 1; u1y = u1z = 0; 2587 } else { 2588 u1x = u1y = 0; u1z = 1; 2589 } 2590 vec_prod(u2x,u2y,u2z, vx,vy,vz, u1x,u1y,u1z); 2591 2592 /* handle case where v and aim are parallel */ 2593 if (!u2x && !u2y && !u2z) { u2x=u2z=0; u2y=1; } 2594 2595 /* u1 = rotate 'v' by theta around u2: DS scattering angle, nearly in horz plane */ 2596 rotate(u1x,u1y,u1z, vx,vy,vz, theta, u2x,u2y,u2z); 2597 2598 /* u0 = rotate u1 by alpha0 around 'v' (Debye-Scherrer cone) */ 2599 rotate(u0x,u0y,u0z, u1x,u1y,u1z, alpha0, vx, vy, vz); 2600 NORM(u0x,u0y,u0z); 2601 vx = u0x*vf; 2602 vy = u0y*vf; 2603 vz = u0z*vf; 2604 2605 SCATTER; 2606 2607 v = vf; k = kf; /* for next iteration */ 2608 2609 } /* end if (flag) */ 2610 2611 VarSqw.neutron_exit++; 2612 p *= p_mult; 2613 if (p_mult > 1) VarSqw.neutron_pmult++; 2614 2615 /* test for a given multiple order */ 2616 if (order && SCATTERED >= order) { 2617 intersect=0; /* reached required number of SCATTERing */ 2618 break; /* finish multiple scattering loop */ 2619 } 2620 2621 } /* end if (intersect) scattering event */ 2622 2623} while (intersect); /* end do (intersect) (multiple scattering loop) */ 2624 2625/* Store Final neutron state */ 2626VarSqw.kf_x = V2K*vx; 2627VarSqw.kf_y = V2K*vy; 2628VarSqw.kf_z = V2K*vz; 2629VarSqw.tf = t; 2630VarSqw.vf = v; 2631VarSqw.kf = k; 2632VarSqw.theta= theta; 2633 2634if (SCATTERED) { 2635 2636 2637 2638 if (SCATTERED == 1) { 2639 if (VarSqw.type == 'c') VarSqw.single_coh += p; 2640 else VarSqw.single_inc += p; 2641 VarSqw.dq = sqrt((VarSqw.kf_x-VarSqw.ki_x)*(VarSqw.kf_x-VarSqw.ki_x) 2642 +(VarSqw.kf_y-VarSqw.ki_y)*(VarSqw.kf_y-VarSqw.ki_y) 2643 +(VarSqw.kf_z-VarSqw.ki_z)*(VarSqw.kf_z-VarSqw.ki_z)); 2644 VarSqw.dw = VS2E*(VarSqw.vf*VarSqw.vf - VarSqw.vi*VarSqw.vi); 2645 } else VarSqw.multi += p; 2646 2647} else VarSqw.dq=VarSqw.dw=0; 2648 2649/* end TRACE */ 2650%} 2651 2652FINALLY 2653%{ 2654 int k; 2655 2656 if (VarSqw.s_coh > 0 || VarSqw.s_inc > 0) 2657 for (k=0; k < 2; k++) { 2658 struct Sqw_Data_struct Data_sqw; 2659 2660 Data_sqw = (k == 0 ? VarSqw.Data_coh : VarSqw.Data_inc); 2661 /* Data_sqw->Sqw has already been freed at end of INIT */ 2662 Table_Free(&(Data_sqw.iqSq)); 2663 2664 if (Data_sqw.SW) free(Data_sqw.SW); 2665 if (Data_sqw.SQW) free(Data_sqw.SQW); 2666 if (Data_sqw.SW_lookup) free(Data_sqw.SW_lookup); 2667 if (Data_sqw.QW_lookup) free(Data_sqw.QW_lookup); 2668 } /* end for */ 2669 2670#ifdef USE_MPI 2671 if (mpi_node_count > 1) { 2672 double tmp; 2673 tmp = (double)VarSqw.neutron_removed; mc_MPI_Sum(&tmp, 1); VarSqw.neutron_removed=(long)tmp; 2674 tmp = (double)VarSqw.neutron_exit; mc_MPI_Sum(&tmp, 1); VarSqw.neutron_exit=(long)tmp; 2675 tmp = (double)VarSqw.neutron_pmult; mc_MPI_Sum(&tmp, 1); VarSqw.neutron_pmult=(long)tmp; 2676 mc_MPI_Sum(&VarSqw.mean_scatt, 1); 2677 mc_MPI_Sum(&VarSqw.psum_scatt, 1); 2678 mc_MPI_Sum(&VarSqw.mean_abs, 1); 2679 mc_MPI_Sum(&VarSqw.single_coh, 1); 2680 mc_MPI_Sum(&VarSqw.single_inc, 1); 2681 mc_MPI_Sum(&VarSqw.multi, 1); 2682 } 2683#endif 2684 MPI_MASTER( 2685 if (VarSqw.neutron_removed) 2686 printf("Isotropic_Sqw: %s: %li neutron events (out of %li) that should have\n" 2687 " scattered were transmitted because scattering conditions\n" 2688 "WARNING could not be satisfied after %i tries.\n", 2689 NAME_CURRENT_COMP, VarSqw.neutron_removed, 2690 VarSqw.neutron_exit+VarSqw.neutron_removed, VarSqw.maxloop); 2691 if (VarSqw.neutron_pmult) 2692 printf("Isotropic_Sqw: %s: %li neutron events (out of %li) reached\n" 2693 "WARNING unrealistic weight. The S(q,w) norm might be too high.\n", 2694 NAME_CURRENT_COMP, VarSqw.neutron_pmult, VarSqw.neutron_exit); 2695 2696 if (VarSqw.verbose_output >= 1 && VarSqw.psum_scatt > 0) { 2697 printf("Isotropic_Sqw: %s: Scattering fraction=%g of incoming intensity\n" 2698 " Absorption fraction =%g\n", 2699 NAME_CURRENT_COMP, 2700 VarSqw.mean_scatt/VarSqw.psum_scatt, VarSqw.mean_abs/VarSqw.psum_scatt); 2701 printf(" Single scattering intensity =%g (coh=%g inc=%g)\n" 2702 " Multiple scattering intensity =%g\n", 2703 VarSqw.single_coh+VarSqw.single_inc, VarSqw.single_coh, VarSqw.single_inc, VarSqw.multi); 2704 ); 2705 } 2706 2707/* end FINALLY */ 2708%} 2709/*****************************************************************************/ 2710/*****************************************************************************/ 2711 2712MCDISPLAY 2713%{ 2714 if (VarSqw.s_coh > 0 || VarSqw.s_inc > 0) { 2715 2716 if(VarSqw.shape==1) 2717 { 2718 double xmin = -0.5*xwidth; 2719 double xmax = 0.5*xwidth; 2720 double ymin = -0.5*yheight; 2721 double ymax = 0.5*yheight; 2722 double zmin = -0.5*zdepth; 2723 double zmax = 0.5*zdepth; 2724 multiline(5, xmin, ymin, zmin, 2725 xmax, ymin, zmin, 2726 xmax, ymax, zmin, 2727 xmin, ymax, zmin, 2728 xmin, ymin, zmin); 2729 multiline(5, xmin, ymin, zmax, 2730 xmax, ymin, zmax, 2731 xmax, ymax, zmax, 2732 xmin, ymax, zmax, 2733 xmin, ymin, zmax); 2734 line(xmin, ymin, zmin, xmin, ymin, zmax); 2735 line(xmax, ymin, zmin, xmax, ymin, zmax); 2736 line(xmin, ymax, zmin, xmin, ymax, zmax); 2737 line(xmax, ymax, zmin, xmax, ymax, zmax); 2738 2739 if (thickness) { 2740 xmin = -0.5*xwidth+thickness; 2741 xmax = -xmin; 2742 ymin = -0.5*yheight+thickness; 2743 ymax = -ymin; 2744 zmin = -0.5*zdepth+thickness; 2745 zmax = -zmin; 2746 multiline(5, xmin, ymin, zmin, 2747 xmax, ymin, zmin, 2748 xmax, ymax, zmin, 2749 xmin, ymax, zmin, 2750 xmin, ymin, zmin); 2751 multiline(5, xmin, ymin, zmax, 2752 xmax, ymin, zmax, 2753 xmax, ymax, zmax, 2754 xmin, ymax, zmax, 2755 xmin, ymin, zmax); 2756 line(xmin, ymin, zmin, xmin, ymin, zmax); 2757 line(xmax, ymin, zmin, xmax, ymin, zmax); 2758 line(xmin, ymax, zmin, xmin, ymax, zmax); 2759 line(xmax, ymax, zmin, xmax, ymax, zmax); 2760 } 2761 } 2762 else if(VarSqw.shape==0) 2763 { 2764 circle("xz", 0, yheight/2.0, 0, radius); 2765 circle("xz", 0, -yheight/2.0, 0, radius); 2766 line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0); 2767 line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0); 2768 line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius); 2769 line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius); 2770 if (thickness) { 2771 double radius_i=radius-thickness; 2772 circle("xz", 0, yheight/2.0, 0, radius_i); 2773 circle("xz", 0, -yheight/2.0, 0, radius_i); 2774 line(-radius_i, -yheight/2.0, 0, -radius_i, +yheight/2.0, 0); 2775 line(+radius_i, -yheight/2.0, 0, +radius_i, +yheight/2.0, 0); 2776 line(0, -yheight/2.0, -radius_i, 0, +yheight/2.0, -radius_i); 2777 line(0, -yheight/2.0, +radius_i, 0, +yheight/2.0, +radius_i); 2778 } 2779 } else if(VarSqw.shape==2) { 2780 if (thickness) { 2781 double radius_i=radius-thickness; 2782 circle("xy",0,0,0,radius_i); 2783 circle("xz",0,0,0,radius_i); 2784 circle("yz",0,0,0,radius_i); 2785 } 2786 circle("xy",0,0,0,radius); 2787 circle("xz",0,0,0,radius); 2788 circle("yz",0,0,0,radius); 2789 } else if (VarSqw.shape == 3) { /* OFF file */ 2790 off_display(offdata); 2791 } 2792 } 2793/* end MCDISPLAY */ 2794%} 2795 2796/*****************************************************************************/ 2797/*****************************************************************************/ 2798 2799END 2800