FOX/ObjCryst++  2022
ReflectionProfile.cpp
1 /* ObjCryst++ Object-Oriented Crystallographic Library
2  (c) 2000-2005 Vincent Favre-Nicolin vincefn@users.sourceforge.net
3  2000-2001 University of Geneva (Switzerland)
4 
5  This program is free software; you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 */
19 /*
20 * source file for LibCryst++ ReflectionProfile and derived classes
21 *
22 */
23 #include <limits>
24 #include "ObjCryst/ObjCryst/ReflectionProfile.h"
25 #include "ObjCryst/Quirks/VFNStreamFormat.h"
26 #ifdef __WX__CRYST__
27  #include "ObjCryst/wxCryst/wxPowderPattern.h"
28 #endif
29 
30 #ifdef HAVE_SSE_MATHFUN
31 #include "ObjCryst/Quirks/sse_mathfun.h"
32 #endif
33 
34 namespace ObjCryst
35 {
36 #if defined(_MSC_VER) || defined(__BORLANDC__)
37 #undef min // Predefined macros.... (wx?)
38 #undef max
39 
40 double erfc(const double x)// in C99, but not in VC++....
41 {
42  if(x<0.0) return 2.0-erfc(-x);
43  if(x<3.8)
44  { // Series, Abramowitz & Stegun 7.1.6
45  double y=x,y0=x;
46  for(int i=1;i<=50;i++)
47  {
48  y0*=2*x*x/(2*i+1.0);
49  y+=y0;
50  }
51  static const double spi=2/sqrt(M_PI);
52  return 1-spi*exp(-x*x)*y;
53  }
54  double y=1.0,y0=1.0;
55  for(int i=1;i<=10;i++)
56  {// Asymptotic, Abramowitz & Stegun 7.1.23
57  y0*=-(2*i-1)/(2*x*x);
58  y+=y0;
59  }
60  static const double invsqrtpi=1.0/sqrt(M_PI);
61  return invsqrtpi*exp(-x*x)/x*y;
62 }
63 
64 #endif
65 extern const RefParType *gpRefParTypeScattDataProfile;
66 extern const RefParType *gpRefParTypeScattDataProfileType;
67 extern const RefParType *gpRefParTypeScattDataProfileWidth;
68 extern const RefParType *gpRefParTypeScattDataProfileAsym;
69 
70 ObjRegistry<ReflectionProfile>
71  gReflectionProfileRegistry("List of all ReflectionProfile types");;
73 //
74 // ReflectionProfile
75 //
77 ReflectionProfile::ReflectionProfile():
78 RefinableObj()
79 {}
80 ReflectionProfile::ReflectionProfile(const ReflectionProfile &old)
81 {}
82 ReflectionProfile::~ReflectionProfile()
83 {}
84 bool ReflectionProfile::IsAnisotropic()const
85 {return false;}
87 //
88 // ReflectionProfilePseudoVoigt
89 //
91 ReflectionProfilePseudoVoigt::ReflectionProfilePseudoVoigt():
92 ReflectionProfile(),
93 mCagliotiU(0),mCagliotiV(0),mCagliotiW(.01*DEG2RAD*DEG2RAD),
94 mPseudoVoigtEta0(0.5),mPseudoVoigtEta1(0.0),
95 mAsymBerarBaldinozziA0(0.0),mAsymBerarBaldinozziA1(0.0),
96 mAsymBerarBaldinozziB0(0.0),mAsymBerarBaldinozziB1(0.0),
97 mAsym0(1.0),mAsym1(0.0),mAsym2(0.0)
98 {
99  this->InitParameters();
100 }
101 
102 ReflectionProfilePseudoVoigt::ReflectionProfilePseudoVoigt
103  (const ReflectionProfilePseudoVoigt &old):
104 mCagliotiU(old.mCagliotiU),mCagliotiV(old.mCagliotiV),mCagliotiW(old.mCagliotiW),
105 mPseudoVoigtEta0(old.mPseudoVoigtEta0),mPseudoVoigtEta1(old.mPseudoVoigtEta1),
106 mAsymBerarBaldinozziA0(old.mAsymBerarBaldinozziA0),
107 mAsymBerarBaldinozziA1(old.mAsymBerarBaldinozziA1),
108 mAsymBerarBaldinozziB0(old.mAsymBerarBaldinozziB0),
109 mAsymBerarBaldinozziB1(old.mAsymBerarBaldinozziB1),
110 mAsym0(old.mAsym0),mAsym1(old.mAsym1),mAsym2(old.mAsym2)
111 {
112  this->InitParameters();
113 }
114 
115 ReflectionProfilePseudoVoigt::~ReflectionProfilePseudoVoigt()
116 {
117  #ifdef __WX__CRYST__
118  if(mpWXCrystObj!=0)
119  {
120  delete mpWXCrystObj;
121  mpWXCrystObj=0;
122  }
123  #endif
124 }
125 
126 ReflectionProfilePseudoVoigt* ReflectionProfilePseudoVoigt::CreateCopy()const
127 {
128  return new ReflectionProfilePseudoVoigt(*this);
129 }
130 
131 const string& ReflectionProfilePseudoVoigt::GetClassName()const
132 {
133  static string className="ReflectionProfilePseudoVoigt";
134  return className;
135 }
136 
137 CrystVector_REAL ReflectionProfilePseudoVoigt::GetProfile(const CrystVector_REAL &x,
138  const REAL center,const REAL h, const REAL k, const REAL l)const
139 {
140  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigt::GetProfile(),c="<<center,2)
141  REAL fwhm= mCagliotiW
142  +mCagliotiV*tan(center/2.0)
143  +mCagliotiU*pow(tan(center/2.0),2);
144  if(fwhm<=0)
145  {
146  VFN_DEBUG_MESSAGE("ReflectionProfilePseudoVoigt::GetProfile(): fwhm**2<0 ! "
147  <<h<<","<<k<<","<<l<<":"<<center<<","<<mCagliotiU<<","<<mCagliotiV<<","<<","<<mCagliotiW<<":"<<fwhm,10);
148  fwhm=1e-6;
149  }
150  else fwhm=sqrt(fwhm);
151  CrystVector_REAL profile,tmpV;
152  REAL asym=mAsym0+mAsym1/sin(center)+mAsym2/pow((REAL)sin(center),(REAL)2.0);
153  if(asym<0.05) asym=0.05;
154  if(asym>20.0) asym=20.0;
155  profile=PowderProfileGauss(x,fwhm,center,asym);
156 
157  // Eta for gaussian/lorentzian mix. Make sure 0<=eta<=1, else profiles could be <0 !
158  REAL eta=mPseudoVoigtEta0+center*mPseudoVoigtEta1;
159  if(eta>1) eta=1;
160  if(eta<0) eta=0;
161 
162  profile *= 1-eta;
163  tmpV=PowderProfileLorentz(x,fwhm,center,asym);
164  tmpV *= eta;
165  profile += tmpV;
166  //profile *= AsymmetryBerarBaldinozzi(x,fwhm,center,
167  // mAsymBerarBaldinozziA0,mAsymBerarBaldinozziA1,
168  // mAsymBerarBaldinozziB0,mAsymBerarBaldinozziB1);
169  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigt::GetProfile()",2)
170  return profile;
171 }
172 
173 void ReflectionProfilePseudoVoigt::SetProfilePar(const REAL fwhmCagliotiW,
174  const REAL fwhmCagliotiU,
175  const REAL fwhmCagliotiV,
176  const REAL eta0,
177  const REAL eta1)
178 {
179  mCagliotiU=fwhmCagliotiU;
180  mCagliotiV=fwhmCagliotiV;
181  mCagliotiW=fwhmCagliotiW;
182  mPseudoVoigtEta0=eta0;
183  mPseudoVoigtEta1=eta1;
184  mClockMaster.Click();
185 }
186 
187 bool ReflectionProfilePseudoVoigt::IsAnisotropic()const{return false;}
188 REAL ReflectionProfilePseudoVoigt::GetFullProfileWidth(const REAL relativeIntensity,
189  const REAL center,const REAL h, const REAL k, const REAL l)
190 {
191  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigt::GetFullProfileWidth()",2)
192  const int nb=100;
193  const int halfnb=nb/2;
194  CrystVector_REAL x(nb);
195  REAL n=5.0;
196  REAL fwhm= mCagliotiW
197  +mCagliotiV*tan(center/2.0)
198  +mCagliotiU*pow(tan(center/2.0),2);
199  if(fwhm<=0) fwhm=1e-6;
200  else fwhm=sqrt(fwhm);
201  CrystVector_REAL prof;
202  while(true)
203  {
204  //Create an X array with 100 elements reaching +/- n*FWHM/2
205  REAL *p=x.data();
206  const REAL tmp=fwhm*n/nb;
207  for(int i=0;i<nb;i++) *p++ = tmp*(i-halfnb);
208  x+=center;
209 
210  prof=this->GetProfile(x,center,0,0,0);
211  const REAL max=prof.max();
212  const REAL test=max*relativeIntensity;
213  int n1=0,n2=0;
214  if((prof(0)<test)&&(prof(nb-1)<test))
215  {
216  p=prof.data();
217  while(*p<test){ p++; n1++;n2++;}
218  n1--;
219  while(*p>test){ p++; n2++;}
220  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigt::GetFullProfileWidth():"<<x(n2)-x(n1),2)
221  return x(n2)-x(n1);
222  }
223  VFN_DEBUG_MESSAGE("ReflectionProfilePseudoVoigt::GetFullProfileWidth():"<<relativeIntensity<<","
224  <<fwhm<<","<<center<<","<<h<<","<<k<<","<<l<<","<<max<<","<<test,2)
225  VFN_DEBUG_MESSAGE(FormatVertVector<REAL>(x,prof),2)
226  n*=2.0;
227  //if(n>200) exit(0);
228  }
229 }
230 
231 void ReflectionProfilePseudoVoigt::InitParameters()
232 {
233  {
234  RefinablePar tmp("U",&mCagliotiU,-1/RAD2DEG/RAD2DEG,1./RAD2DEG/RAD2DEG,
235  gpRefParTypeScattDataProfileWidth,
236  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG*RAD2DEG);
237  tmp.AssignClock(mClockMaster);
238  tmp.SetDerivStep(1e-9);
239  this->AddPar(tmp);
240  }
241  {
242  RefinablePar tmp("V",&mCagliotiV,-1/RAD2DEG/RAD2DEG,1./RAD2DEG/RAD2DEG,
243  gpRefParTypeScattDataProfileWidth,
244  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG*RAD2DEG);
245  tmp.AssignClock(mClockMaster);
246  tmp.SetDerivStep(1e-9);
247  this->AddPar(tmp);
248  }
249  {
250  RefinablePar tmp("W",&mCagliotiW,0,1./RAD2DEG/RAD2DEG,
251  gpRefParTypeScattDataProfileWidth,
252  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG*RAD2DEG);
253  tmp.AssignClock(mClockMaster);
254  tmp.SetDerivStep(1e-9);
255  this->AddPar(tmp);
256  }
257  {
258  RefinablePar tmp("Eta0",&mPseudoVoigtEta0,0,1.,gpRefParTypeScattDataProfileType,
259  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
260  tmp.AssignClock(mClockMaster);
261  tmp.SetDerivStep(1e-4);
262  this->AddPar(tmp);
263  }
264  {
265  RefinablePar tmp("Eta1",&mPseudoVoigtEta1,-1,1.,gpRefParTypeScattDataProfileType,
266  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
267  tmp.AssignClock(mClockMaster);
268  tmp.SetDerivStep(1e-4);
269  this->AddPar(tmp);
270  }
271  #if 0
272  {
273  RefinablePar tmp("AsymA0",&mAsymBerarBaldinozziA0,-0.05,0.05,gpRefParTypeScattDataProfileAsym,
274  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
275  tmp.AssignClock(mClockMaster);
276  tmp.SetDerivStep(1e-4);
277  this->AddPar(tmp);
278  }
279  {
280  RefinablePar tmp("AsymA1",&mAsymBerarBaldinozziA1,-0.05,0.05,gpRefParTypeScattDataProfileAsym,
281  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
282  tmp.AssignClock(mClockMaster);
283  tmp.SetDerivStep(1e-4);
284  this->AddPar(tmp);
285  }
286  {
287  RefinablePar tmp("AsymB0",&mAsymBerarBaldinozziB0,-0.01,0.01,gpRefParTypeScattDataProfileAsym,
288  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
289  tmp.AssignClock(mClockMaster);
290  tmp.SetDerivStep(1e-4);
291  this->AddPar(tmp);
292  }
293  {
294  RefinablePar tmp("AsymB1",&mAsymBerarBaldinozziB1,-0.01,0.01,gpRefParTypeScattDataProfileAsym,
295  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
296  tmp.AssignClock(mClockMaster);
297  tmp.SetDerivStep(1e-4);
298  this->AddPar(tmp);
299  }
300  #endif
301  {
302  RefinablePar tmp("Asym0",&mAsym0,0.01,10.0,gpRefParTypeScattDataProfileAsym,
303  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
304  tmp.AssignClock(mClockMaster);
305  tmp.SetDerivStep(1e-4);
306  this->AddPar(tmp);
307  }
308  {
309  RefinablePar tmp("Asym1",&mAsym1,-1.0,1.0,gpRefParTypeScattDataProfileAsym,
310  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
311  tmp.AssignClock(mClockMaster);
312  tmp.SetDerivStep(1e-4);
313  this->AddPar(tmp);
314  }
315  {
316  RefinablePar tmp("Asym2",&mAsym2,-1.0,1.0,gpRefParTypeScattDataProfileAsym,
317  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
318  tmp.AssignClock(mClockMaster);
319  tmp.SetDerivStep(1e-4);
320  this->AddPar(tmp);
321  }
322 }
323 
324 void ReflectionProfilePseudoVoigt::XMLOutput(ostream &os,int indent)const
325 {
326  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigt::XMLOutput():"<<this->GetName(),5)
327  for(int i=0;i<indent;i++) os << " " ;
328  XMLCrystTag tag("ReflectionProfilePseudoVoigt");
329  os <<tag<<endl;
330  indent++;
331 
332  this->GetPar(&mCagliotiU).XMLOutput(os,"U",indent);
333  os <<endl;
334 
335  this->GetPar(&mCagliotiV).XMLOutput(os,"V",indent);
336  os <<endl;
337 
338  this->GetPar(&mCagliotiW).XMLOutput(os,"W",indent);
339  os <<endl;
340 
341  this->GetPar(&mPseudoVoigtEta0).XMLOutput(os,"Eta0",indent);
342  os <<endl;
343 
344  this->GetPar(&mPseudoVoigtEta1).XMLOutput(os,"Eta1",indent);
345  os <<endl;
346 
347  this->GetPar(&mAsym0).XMLOutput(os,"Asym0",indent);
348  os <<endl;
349 
350  this->GetPar(&mAsym1).XMLOutput(os,"Asym1",indent);
351  os <<endl;
352 
353  this->GetPar(&mAsym2).XMLOutput(os,"Asym2",indent);
354  os <<endl;
355  #if 0
356  this->GetPar(&mAsymBerarBaldinozziA0).XMLOutput(os,"AsymA0",indent);
357  os <<endl;
358 
359  this->GetPar(&mAsymBerarBaldinozziA1).XMLOutput(os,"AsymA1",indent);
360  os <<endl;
361 
362  this->GetPar(&mAsymBerarBaldinozziB0).XMLOutput(os,"AsymB0",indent);
363  os <<endl;
364 
365  this->GetPar(&mAsymBerarBaldinozziB1).XMLOutput(os,"AsymB1",indent);
366  os <<endl;
367  #endif
368  indent--;
369  tag.SetIsEndTag(true);
370  for(int i=0;i<indent;i++) os << " " ;
371  os <<tag<<endl;
372  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigt::XMLOutput():"<<this->GetName(),5)
373 }
374 void ReflectionProfilePseudoVoigt::XMLInput(istream &is,const XMLCrystTag &tagg)
375 {
376  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigt::XMLInput():"<<this->GetName(),5)
377  for(unsigned int i=0;i<tagg.GetNbAttribute();i++)
378  {
379  if("Name"==tagg.GetAttributeName(i)) this->SetName(tagg.GetAttributeValue(i));
380  }
381  while(true)
382  {
383  XMLCrystTag tag(is);
384  if(("ReflectionProfilePseudoVoigt"==tag.GetName())&&tag.IsEndTag())
385  {
386  this->UpdateDisplay();
387  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigt::Exit():"<<this->GetName(),5)
388  return;
389  }
390  if("Par"==tag.GetName())
391  {
392  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
393  {
394  if("Name"==tag.GetAttributeName(i))
395  {
396  if("U"==tag.GetAttributeValue(i))
397  {
398  this->GetPar(&mCagliotiU).XMLInput(is,tag);
399  break;
400  }
401  if("V"==tag.GetAttributeValue(i))
402  {
403  this->GetPar(&mCagliotiV).XMLInput(is,tag);
404  break;
405  }
406  if("W"==tag.GetAttributeValue(i))
407  {
408  this->GetPar(&mCagliotiW).XMLInput(is,tag);
409  break;
410  }
411  if("Eta0"==tag.GetAttributeValue(i))
412  {
413  this->GetPar(&mPseudoVoigtEta0).XMLInput(is,tag);
414  break;
415  }
416  if("Eta1"==tag.GetAttributeValue(i))
417  {
418  this->GetPar(&mPseudoVoigtEta1).XMLInput(is,tag);
419  break;
420  }
421  if("Asym0"==tag.GetAttributeValue(i))
422  {
423  this->GetPar(&mAsym0).XMLInput(is,tag);
424  break;
425  }
426  if("Asym1"==tag.GetAttributeValue(i))
427  {
428  this->GetPar(&mAsym1).XMLInput(is,tag);
429  break;
430  }
431  if("Asym2"==tag.GetAttributeValue(i))
432  {
433  this->GetPar(&mAsym2).XMLInput(is,tag);
434  break;
435  }
436  #if 0
437  if("AsymA0"==tag.GetAttributeValue(i))
438  {
439  this->GetPar(&mAsymBerarBaldinozziA0).XMLInput(is,tag);
440  break;
441  }
442  if("AsymA1"==tag.GetAttributeValue(i))
443  {
444  this->GetPar(&mAsymBerarBaldinozziA1).XMLInput(is,tag);
445  break;
446  }
447  if("AsymB0"==tag.GetAttributeValue(i))
448  {
449  this->GetPar(&mAsymBerarBaldinozziB0).XMLInput(is,tag);
450  break;
451  }
452  if("AsymB1"==tag.GetAttributeValue(i))
453  {
454  this->GetPar(&mAsymBerarBaldinozziB1).XMLInput(is,tag);
455  break;
456  }
457  #endif
458  }
459  }
460  continue;
461  }
462  if("Option"==tag.GetName())
463  {
464  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
465  if("Name"==tag.GetAttributeName(i))
466  mOptionRegistry.GetObj(tag.GetAttributeValue(i)).XMLInput(is,tag);
467  continue;
468  }
469  }
470 }
471 #ifdef __WX__CRYST__
472 WXCrystObjBasic* ReflectionProfilePseudoVoigt::WXCreate(wxWindow* parent)
473 {
474  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigt::WXCreate()",6)
475  if(mpWXCrystObj==0)
476  mpWXCrystObj=new WXProfilePseudoVoigt(parent,this);
477  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigt::WXCreate()",6)
478  return mpWXCrystObj;
479 }
480 #endif
481 
483 //
484 // ReflectionProfilePseudoVoigtAnisotropic
485 //
487 
488 ReflectionProfilePseudoVoigtAnisotropic::ReflectionProfilePseudoVoigtAnisotropic():
489 mCagliotiU(0),mCagliotiV(0),mCagliotiW(.01*DEG2RAD*DEG2RAD),mScherrerP(0),mLorentzX(0),mLorentzY(0),
490 mLorentzGammaHH(0),mLorentzGammaKK(0),mLorentzGammaLL(0),mLorentzGammaHK(0),mLorentzGammaHL(0),mLorentzGammaKL(0),
491 mPseudoVoigtEta0(0.5),mPseudoVoigtEta1(0),mAsym0(1.0),mAsym1(0),mAsym2(0)
492 {
493  this->InitParameters();
494 }
495 
496 ReflectionProfilePseudoVoigtAnisotropic::ReflectionProfilePseudoVoigtAnisotropic(const ReflectionProfilePseudoVoigtAnisotropic &old):
497 mCagliotiU(old.mCagliotiU),mCagliotiV(old.mCagliotiV),mCagliotiW(old.mCagliotiW),mScherrerP(old.mScherrerP),mLorentzX(old.mLorentzX),mLorentzY(old.mLorentzY),
498 mLorentzGammaHH(old.mLorentzGammaHH),mLorentzGammaKK(old.mLorentzGammaKK),mLorentzGammaLL(old.mLorentzGammaLL),mLorentzGammaHK(old.mLorentzGammaHK),mLorentzGammaHL(old.mLorentzGammaHL),mLorentzGammaKL(old.mLorentzGammaKL),
499 mPseudoVoigtEta0(old.mPseudoVoigtEta0),mPseudoVoigtEta1(old.mPseudoVoigtEta1),mAsym0(old.mAsym0),mAsym1(old.mAsym1),mAsym2(old.mAsym2)
500 {
501  this->InitParameters();
502 }
503 ReflectionProfilePseudoVoigtAnisotropic::~ReflectionProfilePseudoVoigtAnisotropic()
504 {
505  #ifdef __WX__CRYST__
506  if(mpWXCrystObj!=0)
507  {
508  delete mpWXCrystObj;
509  mpWXCrystObj=0;
510  }
511  #endif
512 }
513 
514 ReflectionProfilePseudoVoigtAnisotropic* ReflectionProfilePseudoVoigtAnisotropic::CreateCopy()const
515 {
516  return new ReflectionProfilePseudoVoigtAnisotropic(*this);
517 }
518 
519 const string& ReflectionProfilePseudoVoigtAnisotropic::GetClassName()const
520 {
521  static string className="ReflectionProfilePseudoVoigtAnisotropic";
522  return className;
523 }
524 
525 CrystVector_REAL ReflectionProfilePseudoVoigtAnisotropic::GetProfile(const CrystVector_REAL &x, const REAL center,
526  const REAL h, const REAL k, const REAL l)const
527 {
528  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigtAnisotropic::GetProfile()",2)
529  const REAL tantheta=tan(center/2.0);
530  const REAL costheta=cos(center/2.0);
531  const REAL sintheta=sin(center/2.0);
532  const REAL fwhmG=sqrt(abs( mCagliotiW+mCagliotiV*tantheta+mCagliotiU*tantheta*tantheta+mScherrerP/(costheta*costheta)));
533  const REAL gam=mLorentzGammaHH*h*h+mLorentzGammaKK*k*k+mLorentzGammaLL*l*l+2*mLorentzGammaHK*h*k+2*mLorentzGammaHL*h*l+2*mLorentzGammaKL*k*l;
534  const REAL fwhmL= mLorentzX/costheta+(mLorentzY+gam/(sintheta*sintheta))*tantheta;
535  // Eta for gaussian/lorentzian mix. Make sure 0<=eta<=1, else profiles could be <0 !
536  REAL eta=mPseudoVoigtEta0+center*mPseudoVoigtEta1;
537  if(eta>1) eta=1;
538  if(eta<0) eta=0;
539 
540  CrystVector_REAL profile(x.numElements()),tmpV(x.numElements());
541  REAL asym=mAsym0+mAsym1/sin(center)+mAsym2/pow((REAL)sin(center),(REAL)2.0);
542  if(asym<0.05) asym=0.05;
543  if(asym>20.0) asym=20.0;
544  VFN_DEBUG_MESSAGE("ReflectionProfilePseudoVoigtAnisotropic::GetProfile():("<<int(h)<<","<<int(k)<<","<<int(l)<<"),fwhmG="<<fwhmG<<",fwhmL="<<fwhmL<<",gam="<<gam<<",asym="<<asym<<",center="<<center<<",eta="<<eta, 2)
545  if(fwhmG>0)
546  {
547  profile=PowderProfileGauss(x,fwhmG,center,asym);
548  profile *= 1-eta;
549  }
550  else profile=0;
551  if(fwhmL>0)
552  {
553  tmpV=PowderProfileLorentz(x,fwhmL,center,asym);
554  tmpV *= eta;
555  profile += tmpV;
556  }
557  VFN_DEBUG_MESSAGE(FormatVertVector<REAL>(x,profile),1)
558  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigtAnisotropic::GetProfile()",2)
559  return profile;
560 }
561 
562 void ReflectionProfilePseudoVoigtAnisotropic::SetProfilePar(const REAL fwhmCagliotiW,
563  const REAL fwhmCagliotiU,
564  const REAL fwhmCagliotiV,
565  const REAL fwhmGaussP,
566  const REAL fwhmLorentzX,
567  const REAL fwhmLorentzY,
568  const REAL fwhmLorentzGammaHH,
569  const REAL fwhmLorentzGammaKK,
570  const REAL fwhmLorentzGammaLL,
571  const REAL fwhmLorentzGammaHK,
572  const REAL fwhmLorentzGammaHL,
573  const REAL fwhmLorentzGammaKL,
574  const REAL pseudoVoigtEta0,
575  const REAL pseudoVoigtEta1,
576  const REAL asymA0,
577  const REAL asymA1,
578  const REAL asymA2
579  )
580 {
581  mCagliotiU=fwhmCagliotiU;
582  mCagliotiV=fwhmCagliotiV;
583  mCagliotiW=fwhmCagliotiW;
584  mLorentzX=fwhmLorentzX;
585  mLorentzY=fwhmLorentzY;
586  mLorentzGammaHH=fwhmLorentzGammaHH;
587  mLorentzGammaKK=fwhmLorentzGammaKK;
588  mLorentzGammaLL=fwhmLorentzGammaLL;
589  mLorentzGammaHK=fwhmLorentzGammaHK;
590  mLorentzGammaHL=fwhmLorentzGammaHL;
591  mLorentzGammaKL=fwhmLorentzGammaKL;
592  mPseudoVoigtEta0=pseudoVoigtEta0;
593  mPseudoVoigtEta1=pseudoVoigtEta1;
594  mAsym0=asymA0;
595  mAsym1=asymA1;
596  mAsym2=asymA2;
597  mClockMaster.Click();
598 }
599 
600 REAL ReflectionProfilePseudoVoigtAnisotropic::GetFullProfileWidth(const REAL relativeIntensity, const REAL center,
601  const REAL h, const REAL k, const REAL l)
602 {
603  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigt::GetFullProfileWidth()",2)
604  const int nb=100;
605  const int halfnb=nb/2;
606  CrystVector_REAL x(nb);
607  REAL n=5.0;
608  const REAL tantheta=tan(center/2.0);
609  const REAL costheta=cos(center/2.0);
610  const REAL sintheta=sin(center/2.0);
611  const REAL fwhmG=sqrt(abs( mCagliotiW+mCagliotiV*tantheta+mCagliotiU*tantheta*tantheta+mScherrerP/(costheta*costheta)));
612  const REAL gam=mLorentzGammaHH*h*h+mLorentzGammaKK*k*k+mLorentzGammaLL*l*l+2*mLorentzGammaHK*h*k+2*mLorentzGammaHL*h*l+2*mLorentzGammaKL*k*l;
613  const REAL fwhmL= mLorentzX/costheta+(mLorentzY+gam/(sintheta*sintheta))*tantheta;
614  const REAL eta=mPseudoVoigtEta0+mPseudoVoigtEta1*center;
615  // Obviously this is not the REAL FWHM, just a _very_ crude starting approximation
616  REAL fwhm=fwhmL*eta+fwhmG*(1-eta);
617  if(fwhm<=0) fwhm=1e-3;
618  CrystVector_REAL prof;
619  while(true)
620  {
621  //Create an X array with 100 elements reaching +/- n*FWHM/2
622  REAL *p=x.data();
623  const REAL tmp=fwhm*n/nb;
624  for(int i=0;i<nb;i++) *p++ = tmp*(i-halfnb);
625  x+=center;
626  prof=this->GetProfile(x,center,h,k,l);
627  const REAL max=prof.max();
628  const REAL test=max*relativeIntensity;
629  int n1=0,n2=0;
630  if((prof(0)<test)&&(prof(nb-1)<test))
631  {
632  p=prof.data();
633  while(*p<test){ p++; n1++;n2++;}
634  n1--;
635  while(*p>test){ p++; n2++;}
636  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigtAnisotropic::GetFullProfileWidth():"<<x(n2)-x(n1),2)
637  return x(n2)-x(n1);
638  }
639  VFN_DEBUG_MESSAGE("ReflectionProfilePseudoVoigtAnisotropic::GetFullProfileWidth():"<<relativeIntensity<<","
640  <<fwhm<<","<<center<<","<<h<<","<<k<<","<<l<<","<<max<<","<<test,2)
641  VFN_DEBUG_MESSAGE(FormatVertVector<REAL>(x,prof),1)
642  n*=2.0;
643  }
644 }
645 
646 bool ReflectionProfilePseudoVoigtAnisotropic::IsAnisotropic()const
647 {
648  return true;
649 }
650 
651 void ReflectionProfilePseudoVoigtAnisotropic::XMLOutput(ostream &os,int indent)const
652 {
653  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigtAnisotropic::XMLOutput():"<<this->GetName(),5)
654  for(int i=0;i<indent;i++) os << " " ;
655  XMLCrystTag tag("ReflectionProfilePseudoVoigtAnisotropic");
656  os <<tag<<endl;
657  indent++;
658 
659  this->GetPar(&mCagliotiU).XMLOutput(os,"U",indent);
660  os <<endl;
661 
662  this->GetPar(&mCagliotiV).XMLOutput(os,"V",indent);
663  os <<endl;
664 
665  this->GetPar(&mCagliotiW).XMLOutput(os,"W",indent);
666  os <<endl;
667 
668  this->GetPar(&mScherrerP).XMLOutput(os,"P",indent);
669  os <<endl;
670 
671  this->GetPar(&mLorentzX).XMLOutput(os,"X",indent);
672  os <<endl;
673 
674  this->GetPar(&mLorentzY).XMLOutput(os,"Y",indent);
675  os <<endl;
676 
677  this->GetPar(&mLorentzGammaHH).XMLOutput(os,"G_HH",indent);
678  os <<endl;
679 
680  this->GetPar(&mLorentzGammaKK).XMLOutput(os,"G_KK",indent);
681  os <<endl;
682 
683  this->GetPar(&mLorentzGammaLL).XMLOutput(os,"G_LL",indent);
684  os <<endl;
685 
686  this->GetPar(&mLorentzGammaHK).XMLOutput(os,"G_HK",indent);
687  os <<endl;
688 
689  this->GetPar(&mLorentzGammaHL).XMLOutput(os,"G_HL",indent);
690  os <<endl;
691 
692  this->GetPar(&mLorentzGammaKL).XMLOutput(os,"G_KL",indent);
693  os <<endl;
694 
695  this->GetPar(&mPseudoVoigtEta0).XMLOutput(os,"Eta0",indent);
696  os <<endl;
697 
698  this->GetPar(&mPseudoVoigtEta1).XMLOutput(os,"Eta1",indent);
699  os <<endl;
700 
701  this->GetPar(&mAsym0).XMLOutput(os,"Asym0",indent);
702  os <<endl;
703 
704  this->GetPar(&mAsym1).XMLOutput(os,"Asym1",indent);
705  os <<endl;
706 
707  this->GetPar(&mAsym2).XMLOutput(os,"Asym2",indent);
708  os <<endl;
709  indent--;
710  tag.SetIsEndTag(true);
711  for(int i=0;i<indent;i++) os << " " ;
712  os <<tag<<endl;
713  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigtAnisotropic::XMLOutput():"<<this->GetName(),5)
714 }
715 
716 void ReflectionProfilePseudoVoigtAnisotropic::XMLInput(istream &is,const XMLCrystTag &tagg)
717 {
718  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigtAnisotropic::XMLInput():"<<this->GetName(),5)
719  for(unsigned int i=0;i<tagg.GetNbAttribute();i++)
720  {
721  if("Name"==tagg.GetAttributeName(i)) this->SetName(tagg.GetAttributeValue(i));
722  }
723  while(true)
724  {
725  XMLCrystTag tag(is);
726  if(("ReflectionProfilePseudoVoigtAnisotropic"==tag.GetName())&&tag.IsEndTag())
727  {
728  this->UpdateDisplay();
729  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigtAnisotropic::Exit():"<<this->GetName(),5)
730  return;
731  }
732  if("Par"==tag.GetName())
733  {
734  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
735  {
736  if("Name"==tag.GetAttributeName(i))
737  {
738  if("U"==tag.GetAttributeValue(i))
739  {
740  this->GetPar(&mCagliotiU).XMLInput(is,tag);
741  break;
742  }
743  if("V"==tag.GetAttributeValue(i))
744  {
745  this->GetPar(&mCagliotiV).XMLInput(is,tag);
746  break;
747  }
748  if("W"==tag.GetAttributeValue(i))
749  {
750  this->GetPar(&mCagliotiW).XMLInput(is,tag);
751  break;
752  }
753  if("P"==tag.GetAttributeValue(i))
754  {
755  this->GetPar(&mScherrerP).XMLInput(is,tag);
756  break;
757  }
758  if("X"==tag.GetAttributeValue(i))
759  {
760  this->GetPar(&mLorentzX).XMLInput(is,tag);
761  break;
762  }
763  if("Y"==tag.GetAttributeValue(i))
764  {
765  this->GetPar(&mLorentzY).XMLInput(is,tag);
766  break;
767  }
768  if("G_HH"==tag.GetAttributeValue(i))
769  {
770  this->GetPar(&mLorentzGammaHH).XMLInput(is,tag);
771  break;
772  }
773  if("G_KK"==tag.GetAttributeValue(i))
774  {
775  this->GetPar(&mLorentzGammaKK).XMLInput(is,tag);
776  break;
777  }
778  if("G_LL"==tag.GetAttributeValue(i))
779  {
780  this->GetPar(&mLorentzGammaLL).XMLInput(is,tag);
781  break;
782  }
783  if("G_HK"==tag.GetAttributeValue(i))
784  {
785  this->GetPar(&mLorentzGammaHK).XMLInput(is,tag);
786  break;
787  }
788  if("G_HL"==tag.GetAttributeValue(i))
789  {
790  this->GetPar(&mLorentzGammaHL).XMLInput(is,tag);
791  break;
792  }
793  if("G_KL"==tag.GetAttributeValue(i))
794  {
795  this->GetPar(&mLorentzGammaKL).XMLInput(is,tag);
796  break;
797  }
798  if("Eta0"==tag.GetAttributeValue(i))
799  {
800  this->GetPar(&mPseudoVoigtEta0).XMLInput(is,tag);
801  break;
802  }
803  if("Eta1"==tag.GetAttributeValue(i))
804  {
805  this->GetPar(&mPseudoVoigtEta1).XMLInput(is,tag);
806  break;
807  }
808  if("Asym0"==tag.GetAttributeValue(i))
809  {
810  this->GetPar(&mAsym0).XMLInput(is,tag);
811  break;
812  }
813  if("Asym1"==tag.GetAttributeValue(i))
814  {
815  this->GetPar(&mAsym1).XMLInput(is,tag);
816  break;
817  }
818  if("Asym2"==tag.GetAttributeValue(i))
819  {
820  this->GetPar(&mAsym2).XMLInput(is,tag);
821  break;
822  }
823  }
824  }
825  continue;
826  }
827  if("Option"==tag.GetName())
828  {
829  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
830  if("Name"==tag.GetAttributeName(i))
831  mOptionRegistry.GetObj(tag.GetAttributeValue(i)).XMLInput(is,tag);
832  continue;
833  }
834  }
835 }
836 
837 void ReflectionProfilePseudoVoigtAnisotropic::InitParameters()
838 {
839  {
840  RefinablePar tmp("U",&mCagliotiU,-1/RAD2DEG/RAD2DEG,1./RAD2DEG/RAD2DEG,
841  gpRefParTypeScattDataProfileWidth,
842  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG*RAD2DEG);
843  tmp.AssignClock(mClockMaster);
844  tmp.SetDerivStep(1e-9);
845  this->AddPar(tmp);
846  }
847  {
848  RefinablePar tmp("V",&mCagliotiV,-1/RAD2DEG/RAD2DEG,1./RAD2DEG/RAD2DEG,
849  gpRefParTypeScattDataProfileWidth,
850  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG*RAD2DEG);
851  tmp.AssignClock(mClockMaster);
852  tmp.SetDerivStep(1e-9);
853  this->AddPar(tmp);
854  }
855  {
856  RefinablePar tmp("W",&mCagliotiW,0,1./RAD2DEG/RAD2DEG,
857  gpRefParTypeScattDataProfileWidth,
858  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG*RAD2DEG);
859  tmp.AssignClock(mClockMaster);
860  tmp.SetDerivStep(1e-9);
861  this->AddPar(tmp);
862  }
863  {
864  RefinablePar tmp("P",&mScherrerP,-1./RAD2DEG/RAD2DEG,1./RAD2DEG/RAD2DEG,
865  gpRefParTypeScattDataProfileWidth,
866  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG*RAD2DEG);
867  tmp.AssignClock(mClockMaster);
868  tmp.SetDerivStep(1e-9);
869  this->AddPar(tmp);
870  }
871  {
872  RefinablePar tmp("X",&mLorentzX,0,1./RAD2DEG,
873  gpRefParTypeScattDataProfileWidth,
874  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG);
875  tmp.AssignClock(mClockMaster);
876  tmp.SetDerivStep(1e-9);
877  this->AddPar(tmp);
878  }
879  {
880  RefinablePar tmp("Y",&mLorentzY,-1./RAD2DEG,1./RAD2DEG,
881  gpRefParTypeScattDataProfileWidth,
882  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG);
883  tmp.AssignClock(mClockMaster);
884  tmp.SetDerivStep(1e-9);
885  this->AddPar(tmp);
886  }
887  {
888  RefinablePar tmp("G_HH",&mLorentzGammaHH,-1./RAD2DEG,1./RAD2DEG,
889  gpRefParTypeScattDataProfileWidth,
890  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG);
891  tmp.AssignClock(mClockMaster);
892  tmp.SetDerivStep(1e-9);
893  this->AddPar(tmp);
894  }
895  {
896  RefinablePar tmp("G_KK",&mLorentzGammaKK,-1./RAD2DEG,1./RAD2DEG,
897  gpRefParTypeScattDataProfileWidth,
898  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG);
899  tmp.AssignClock(mClockMaster);
900  tmp.SetDerivStep(1e-9);
901  this->AddPar(tmp);
902  }
903  {
904  RefinablePar tmp("G_LL",&mLorentzGammaLL,-1./RAD2DEG,1./RAD2DEG,
905  gpRefParTypeScattDataProfileWidth,
906  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG);
907  tmp.AssignClock(mClockMaster);
908  tmp.SetDerivStep(1e-9);
909  this->AddPar(tmp);
910  }
911  {
912  RefinablePar tmp("G_HK",&mLorentzGammaHK,-1./RAD2DEG,1./RAD2DEG,
913  gpRefParTypeScattDataProfileWidth,
914  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG);
915  tmp.AssignClock(mClockMaster);
916  tmp.SetDerivStep(1e-9);
917  this->AddPar(tmp);
918  }
919  {
920  RefinablePar tmp("G_HL",&mLorentzGammaHL,-1./RAD2DEG,1./RAD2DEG,
921  gpRefParTypeScattDataProfileWidth,
922  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG);
923  tmp.AssignClock(mClockMaster);
924  tmp.SetDerivStep(1e-9);
925  this->AddPar(tmp);
926  }
927  {
928  RefinablePar tmp("G_KL",&mLorentzGammaKL,-1./RAD2DEG,1./RAD2DEG,
929  gpRefParTypeScattDataProfileWidth,
930  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false,RAD2DEG);
931  tmp.AssignClock(mClockMaster);
932  tmp.SetDerivStep(1e-9);
933  this->AddPar(tmp);
934  }
935  {
936  RefinablePar tmp("Eta0",&mPseudoVoigtEta0,0,1.,gpRefParTypeScattDataProfileType,
937  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
938  tmp.AssignClock(mClockMaster);
939  tmp.SetDerivStep(1e-4);
940  this->AddPar(tmp);
941  }
942  {
943  RefinablePar tmp("Eta1",&mPseudoVoigtEta1,-1,1.,gpRefParTypeScattDataProfileType,
944  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
945  tmp.AssignClock(mClockMaster);
946  tmp.SetDerivStep(1e-4);
947  this->AddPar(tmp);
948  }
949  {
950  RefinablePar tmp("Asym0",&mAsym0,0.01,10.0,gpRefParTypeScattDataProfileAsym,
951  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
952  tmp.AssignClock(mClockMaster);
953  tmp.SetDerivStep(1e-4);
954  this->AddPar(tmp);
955  }
956  {
957  RefinablePar tmp("Asym1",&mAsym1,-1.0,1.0,gpRefParTypeScattDataProfileAsym,
958  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
959  tmp.AssignClock(mClockMaster);
960  tmp.SetDerivStep(1e-4);
961  this->AddPar(tmp);
962  }
963  {
964  RefinablePar tmp("Asym2",&mAsym2,-1.0,1.0,gpRefParTypeScattDataProfileAsym,
965  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
966  tmp.AssignClock(mClockMaster);
967  tmp.SetDerivStep(1e-4);
968  this->AddPar(tmp);
969  }
970 }
971 
972 #ifdef __WX__CRYST__
973 WXCrystObjBasic* ReflectionProfilePseudoVoigtAnisotropic::WXCreate(wxWindow* parent)
974 {
975  VFN_DEBUG_ENTRY("ReflectionProfilePseudoVoigt::WXCreate()",6)
976  if(mpWXCrystObj==0)
977  mpWXCrystObj=new WXProfilePseudoVoigtAnisotropic(parent,this);
978  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigt::WXCreate()",6)
979  return mpWXCrystObj;
980 }
981 #endif
982 
984 //
985 // ReflectionProfileDoubleExponentialPseudoVoigt
986 //
988 ReflectionProfileDoubleExponentialPseudoVoigt::ReflectionProfileDoubleExponentialPseudoVoigt():
989 ReflectionProfile(),
990 mInstrumentAlpha0(0.0),
991 mInstrumentAlpha1(0.0952),
992 mInstrumentBeta0(0.0239),
993 mInstrumentBeta1(0.0043),
994 mGaussianSigma0(0.0),
995 mGaussianSigma1(7.0),
996 mGaussianSigma2(0.0),
997 mLorentzianGamma0(0.0),
998 mLorentzianGamma1(0.0),
999 mLorentzianGamma2(0.414),
1000 mpCell(0)
1001 {
1002  VFN_DEBUG_MESSAGE("ReflectionProfileDoubleExponentialPseudoVoigt::ReflectionProfileDoubleExponentialPseudoVoigt()",10)
1003  this->InitParameters();
1004 }
1005 
1006 ReflectionProfileDoubleExponentialPseudoVoigt
1007  ::ReflectionProfileDoubleExponentialPseudoVoigt(const UnitCell &cell):
1008 ReflectionProfile(),
1009 mInstrumentAlpha0(0.0),
1010 mInstrumentAlpha1(0.0952),
1011 mInstrumentBeta0(0.0239),
1012 mInstrumentBeta1(0.0043),
1013 mGaussianSigma0(0.0),
1014 mGaussianSigma1(7.0),
1015 mGaussianSigma2(0.0),
1016 mLorentzianGamma0(0.0),
1017 mLorentzianGamma1(0.0),
1018 mLorentzianGamma2(0.414),
1019 mpCell(&cell)
1020 {
1021  VFN_DEBUG_MESSAGE("ReflectionProfileDoubleExponentialPseudoVoigt::ReflectionProfileDoubleExponentialPseudoVoigt()",10)
1022  this->InitParameters();
1023 }
1024 
1025 ReflectionProfileDoubleExponentialPseudoVoigt::ReflectionProfileDoubleExponentialPseudoVoigt
1026  (const ReflectionProfileDoubleExponentialPseudoVoigt &old):
1027 ReflectionProfile(),
1028 mInstrumentAlpha0(old.mInstrumentAlpha0),
1029 mInstrumentAlpha1(old.mInstrumentAlpha1),
1030 mInstrumentBeta0(old.mInstrumentBeta0),
1031 mInstrumentBeta1(old.mInstrumentBeta1),
1032 mGaussianSigma0(old.mGaussianSigma0),
1033 mGaussianSigma1(old.mGaussianSigma1),
1034 mGaussianSigma2(old.mGaussianSigma2),
1035 mLorentzianGamma0(old.mLorentzianGamma0),
1036 mLorentzianGamma1(old.mLorentzianGamma1),
1037 mLorentzianGamma2(old.mLorentzianGamma2),
1038 mpCell(old.mpCell)
1039 {
1040  VFN_DEBUG_MESSAGE("ReflectionProfileDoubleExponentialPseudoVoigt::ReflectionProfileDoubleExponentialPseudoVoigt()",10)
1041  this->InitParameters();
1042 }
1043 
1044 ReflectionProfileDoubleExponentialPseudoVoigt::~ReflectionProfileDoubleExponentialPseudoVoigt()
1045 {
1046  #ifdef __WX__CRYST__
1047  if(mpWXCrystObj!=0)
1048  {
1049  delete mpWXCrystObj;
1050  mpWXCrystObj=0;
1051  }
1052  #endif
1053 }
1054 
1055 ReflectionProfileDoubleExponentialPseudoVoigt*
1056  ReflectionProfileDoubleExponentialPseudoVoigt::CreateCopy()const
1057 {
1058  return new ReflectionProfileDoubleExponentialPseudoVoigt(*this);
1059 }
1060 
1061 const string& ReflectionProfileDoubleExponentialPseudoVoigt::GetClassName()const
1062 {
1063  static string className="ReflectionProfileDoubleExponentialPseudoVoigt";
1064  return className;
1065 }
1066 
1067 CrystVector_REAL ReflectionProfileDoubleExponentialPseudoVoigt
1068  ::GetProfile(const CrystVector_REAL &x, const REAL center,
1069  const REAL h, const REAL k, const REAL l)const
1070 {
1071  VFN_DEBUG_ENTRY("ReflectionProfileDoubleExponentialPseudoVoigt::GetProfile()",4)
1072  REAL dcenter=0;
1073  if(mpCell!=0)
1074  {
1075  REAL hh=h,kk=k,ll=l;// orthonormal coordinates in reciprocal space
1076  mpCell->MillerToOrthonormalCoords(hh,kk,ll);
1077  dcenter=1.0/sqrt(hh*hh+kk*kk+ll*ll);//d_hkl, in Angstroems
1078  }
1079  const REAL alpha=mInstrumentAlpha0+mInstrumentAlpha1/dcenter;
1080  const REAL beta=mInstrumentBeta0+mInstrumentBeta1/pow(dcenter,4);
1081  const REAL siggauss2= mGaussianSigma0
1082  +mGaussianSigma1*pow(dcenter,2)
1083  +mGaussianSigma2*pow(dcenter,4);
1084  static const REAL log2=log(2.0);
1085  const REAL hg=sqrt(8*siggauss2*log2);
1086  const REAL hl= mLorentzianGamma0
1087  +mLorentzianGamma1*dcenter
1088  +mLorentzianGamma2*dcenter*dcenter;
1089  const REAL hcom=pow(pow(hg,5)+2.69269*pow(hg,4)*hl+2.42843*pow(hg,3)*hl*hl
1090  +4.47163*hg*hg*pow(hl,3)+0.07842*hg*pow(hl,4)+pow(hl,5),0.2);
1091  const REAL sigcom2=hcom*hcom/(8.0*log2);
1092  const REAL eta=1.36603*hl/hcom-0.47719*pow(hl/hcom,2)+0.11116*pow(hl/hcom,3);
1093  const long nbPoints=x.numElements();
1094  VFN_DEBUG_MESSAGE("ReflectionProfileDoubleExponentialPseudoVoigt::GetProfile():alpha="
1095  <<alpha<<",beta="<<beta<<",siggauss2="<<siggauss2
1096  <<",hg="<<hg<<",hl="<<hl<<",hcom="<<hcom<<",sigcom2="<<sigcom2
1097  <<",eta="<<eta,2)
1098  CrystVector_REAL prof;
1099  prof=x;
1100  prof+=-center;
1101  REAL *pp=prof.data();
1102  for(long i=0;i<nbPoints;i++)
1103  {
1104  const double u=alpha/2*(alpha*sigcom2+2* *pp);
1105  const double nu=beta/2*(beta *sigcom2-2* *pp);
1106  const double y=(alpha*sigcom2+*pp)/sqrt(2*sigcom2);
1107  const double z=(beta *sigcom2-*pp)/sqrt(2*sigcom2);
1108  const complex<double> p(alpha* *pp,alpha*hcom/2);
1109  const complex<double> q(-beta* *pp, beta*hcom/2);
1110  const complex<double> e1p=ExponentialIntegral1_ExpZ(p);
1111  const complex<double> e1q=ExponentialIntegral1_ExpZ(q);
1112  VFN_DEBUG_MESSAGE("dt="<<*pp<<", u="<<u<<",nu="<<nu<<",y="<<y<<",z="<<z
1113  <<",p=("<<p.real()<<","<<p.imag()
1114  <<"),q=("<<q.real()<<","<<q.imag()
1115  <<"),e^p*E1(p)=("<<e1p.real()<<","<<e1p.imag()
1116  <<"),e^q*E1(q)=("<<e1q.real()<<","<<e1q.imag(),2)
1117  REAL expnu_erfcz,expu_erfcy;
1118  // Use asymptotic value for erfc(x) = 1/(sqrt(pi)*x*exp(x^2)) [A&S 7.1.23]
1119  if(z>10.0) expnu_erfcz=exp(nu-z*z)/(z*sqrt(M_PI));
1120  else expnu_erfcz=exp(nu)*erfc(z);
1121 
1122  if(y>10.0) expu_erfcy=exp(u-y*y)/(y*sqrt(M_PI));
1123  else expu_erfcy=exp(u)*erfc(y);
1124 
1125  #if 0
1126  double tmp=(1-eta)*alpha*beta/(2*(alpha+beta))*(expu_erfcy+expnu_erfcz)
1127  -eta*alpha*beta/(M_PI*(alpha+beta))*(e1p.imag()+e1q.imag());
1128  if(isnan(*pp))// Is this portable ? Test for numeric_limits<REAL>::quiet_NaN()
1129  {
1130  cout<<"*pp==numeric_limits<REAL>::quiet_NaN()"<<endl;
1131  cout<<"ReflectionProfileDoubleExponentialPseudoVoigt::GetProfile():"<<endl
1132  <<" alpha="<<alpha<<",beta="<<beta<<",siggauss2="<<siggauss2
1133  <<",hg="<<hg<<",hl="<<hl<<",hcom="<<hcom<<",sigcom2="<<sigcom2
1134  <<",eta="<<eta<<endl;
1135  cout<<" dt="<<*pp<<", u="<<u<<",nu="<<nu<<",y="<<y<<",z="<<z
1136  <<",e^u*E1(y)="<<expu_erfcy
1137  <<",e^nu*E1(z)="<<expnu_erfcz
1138  <<endl
1139  <<" p=("<<p.real()<<","<<p.imag()
1140  <<"),q=("<<q.real()<<","<<q.imag()
1141  <<"),e^p*E1(p)=("<<e1p.real()<<","<<e1p.imag()
1142  <<"),e^q*E1(q)=("<<e1q.real()<<","<<e1q.imag()<<endl;
1143  cout<<(1-eta)*alpha*beta/(2*(alpha+beta))*(expu_erfcy+expnu_erfcz)<<endl
1144  <<eta*alpha*beta/(M_PI*(alpha+beta))*(e1p.imag()+e1q.imag())<<endl
1145  << *pp<<endl
1146  << tmp<<endl;
1147  exit(0);
1148  }
1149  if(abs(*pp)==numeric_limits<REAL>::infinity())
1150  {
1151  cout<<"*pp==numeric_limits<REAL>::infinity()"<<endl;
1152  exit(0);
1153  }
1154  //if(*pp>1e30) exit(0);
1155  #endif
1156  *pp++=(1-eta)*alpha*beta/(2*(alpha+beta))*(expu_erfcy+expnu_erfcz)
1157  -eta*alpha*beta/(M_PI*(alpha+beta))*(e1p.imag()+e1q.imag());
1158  }
1159  VFN_DEBUG_EXIT("ReflectionProfileDoubleExponentialPseudoVoigt::GetProfile()",4)
1160  return prof;
1161 }
1162 
1163 void ReflectionProfileDoubleExponentialPseudoVoigt
1164  ::SetProfilePar(const REAL instrumentAlpha0,
1165  const REAL instrumentAlpha1,
1166  const REAL instrumentBeta0,
1167  const REAL instrumentBeta1,
1168  const REAL gaussianSigma0,
1169  const REAL gaussianSigma1,
1170  const REAL gaussianSigma2,
1171  const REAL lorentzianGamma0,
1172  const REAL lorentzianGamma1,
1173  const REAL lorentzianGamma2)
1174 {
1175  mInstrumentAlpha0=instrumentAlpha0;
1176  mInstrumentAlpha1=instrumentAlpha1;
1177  mInstrumentBeta0=instrumentBeta0;
1178  mInstrumentBeta1=instrumentBeta1;
1179  mGaussianSigma0=gaussianSigma0;
1180  mGaussianSigma1=gaussianSigma1;
1181  mGaussianSigma2=gaussianSigma2;
1182  mLorentzianGamma0=lorentzianGamma0;
1183  mLorentzianGamma1=lorentzianGamma1;
1184  mLorentzianGamma2=lorentzianGamma2;
1185  mClockMaster.Click();
1186 }
1187 
1188 REAL ReflectionProfileDoubleExponentialPseudoVoigt
1189  ::GetFullProfileWidth(const REAL relativeIntensity, const REAL center,
1190  const REAL h, const REAL k, const REAL l)
1191 {
1192  VFN_DEBUG_ENTRY("ReflectionProfileDoubleExponentialPseudoVoigt::GetFullProfileWidth()",5)
1193  REAL dcenter=0;
1194  if(mpCell!=0)
1195  {
1196  REAL hh=h,kk=k,ll=l;// orthonormal coordinates in reciprocal space
1197  VFN_DEBUG_MESSAGE("ReflectionProfileDoubleExponentialPseudoVoigt::GetFullProfileWidth(),"<<dcenter<<","<<mpCell->GetName(),5)
1198  mpCell->MillerToOrthonormalCoords(hh,kk,ll);
1199  VFN_DEBUG_MESSAGE("ReflectionProfileDoubleExponentialPseudoVoigt::GetFullProfileWidth(),"<<dcenter,5)
1200  dcenter=sqrt(hh*hh+kk*kk+ll*ll);//1/d
1201  }
1202  VFN_DEBUG_MESSAGE("ReflectionProfileDoubleExponentialPseudoVoigt::GetFullProfileWidth(),"<<dcenter,5)
1203  const int nb=100;
1204  const int halfnb=nb/2;
1205  CrystVector_REAL x(nb);
1206  REAL n=5.0;
1207  const REAL siggauss2= mGaussianSigma0
1208  +mGaussianSigma1*pow(dcenter,2)
1209  +mGaussianSigma2*pow(dcenter,4);
1210  static const REAL log2=log(2.0);
1211  const REAL hg=sqrt(8*siggauss2*log2);
1212  const REAL hl= mLorentzianGamma0
1213  +mLorentzianGamma1*dcenter
1214  +mLorentzianGamma2*dcenter*dcenter;
1215  const REAL fwhm=pow(pow(hg,5)+2.69269*pow(hg,4)*hl+2.42843*pow(hg,3)*hl*hl
1216  +4.47163*hg*hg*pow(hl,3)+0.07842*hg*pow(hl,4)+pow(hl,5),0.2);
1217  CrystVector_REAL prof;
1218  while(true)
1219  {
1220  REAL *p=x.data();
1221  const REAL tmp=fwhm*n/nb;
1222  for(int i=0;i<nb;i++) *p++ = tmp*(i-halfnb);
1223  x+=center;
1224  prof=this->GetProfile(x,center,h,k,l);
1225  const REAL max=prof.max();
1226  const REAL test=max*relativeIntensity;
1227  int n1=0,n2=0;
1228  if((prof(0)<test)&&(prof(nb-1)<test))
1229  {
1230  p=prof.data();
1231  while(*p<test){ p++; n1++;n2++;}
1232  n1--;
1233  while(*p>test){ p++; n2++;}
1234  VFN_DEBUG_EXIT("ReflectionProfilePseudoVoigt::GetFullProfileWidth():"<<x(n2)-x(n1),5)
1235  return abs(x(n2)-x(n1));
1236  }
1237  VFN_DEBUG_MESSAGE("ReflectionProfilePseudoVoigt::GetFullProfileWidth():"<<max<<","<<test
1238  <<endl<<FormatVertVector<REAL>(x,prof),5)
1239  n*=2.0;
1240  //if(n>200) exit(0);
1241  }
1242  // unreachable code.
1243  // VFN_DEBUG_EXIT("ReflectionProfileDoubleExponentialPseudoVoigt::GetFullProfileWidth()",5)
1244 }
1245 
1246 bool ReflectionProfileDoubleExponentialPseudoVoigt
1247  ::IsAnisotropic()const{return false;}
1248 
1249 void ReflectionProfileDoubleExponentialPseudoVoigt
1250  ::XMLOutput(ostream &os,int indent)const
1251 {
1252  VFN_DEBUG_ENTRY("ReflectionProfileDoubleExponentialPseudoVoigt::XMLOutput():"<<this->GetName(),5)
1253  for(int i=0;i<indent;i++) os << " " ;
1254  XMLCrystTag tag("ReflectionProfileDoubleExponentialPseudoVoigt");
1255  os <<tag<<endl;
1256  indent++;
1257 
1258  this->GetPar(&mInstrumentAlpha0).XMLOutput(os,"Alpha0",indent);
1259  os <<endl;
1260 
1261  this->GetPar(&mInstrumentAlpha1).XMLOutput(os,"Alpha1",indent);
1262  os <<endl;
1263 
1264  this->GetPar(&mInstrumentBeta0).XMLOutput(os,"Beta0",indent);
1265  os <<endl;
1266 
1267  this->GetPar(&mInstrumentBeta1).XMLOutput(os,"Beta1",indent);
1268  os <<endl;
1269 
1270  this->GetPar(&mGaussianSigma0).XMLOutput(os,"GaussianSigma0",indent);
1271  os <<endl;
1272 
1273  this->GetPar(&mGaussianSigma1).XMLOutput(os,"GaussianSigma1",indent);
1274  os <<endl;
1275 
1276  this->GetPar(&mGaussianSigma2).XMLOutput(os,"GaussianSigma2",indent);
1277  os <<endl;
1278 
1279  this->GetPar(&mLorentzianGamma0).XMLOutput(os,"LorentzianGamma0",indent);
1280  os <<endl;
1281 
1282  this->GetPar(&mLorentzianGamma1).XMLOutput(os,"LorentzianGamma1",indent);
1283  os <<endl;
1284 
1285  this->GetPar(&mLorentzianGamma2).XMLOutput(os,"LorentzianGamma2",indent);
1286  os <<endl;
1287 
1288  indent--;
1289  tag.SetIsEndTag(true);
1290  for(int i=0;i<indent;i++) os << " " ;
1291  os <<tag<<endl;
1292  VFN_DEBUG_EXIT("ReflectionProfileDoubleExponentialPseudoVoigt::XMLOutput():"<<this->GetName(),5)
1293 }
1294 
1295 void ReflectionProfileDoubleExponentialPseudoVoigt
1296  ::XMLInput(istream &is,const XMLCrystTag &tagg)
1297 {
1298  VFN_DEBUG_ENTRY("ReflectionProfileDoubleExponentialPseudoVoigt::XMLInput():"<<this->GetName(),5)
1299  for(unsigned int i=0;i<tagg.GetNbAttribute();i++)
1300  {
1301  if("Name"==tagg.GetAttributeName(i)) this->SetName(tagg.GetAttributeValue(i));
1302  }
1303  while(true)
1304  {
1305  XMLCrystTag tag(is);
1306  if(("ReflectionProfileDoubleExponentialPseudoVoigt"==tag.GetName())&&tag.IsEndTag())
1307  {
1308  this->UpdateDisplay();
1309  VFN_DEBUG_EXIT("ReflectionProfileDoubleExponentialPseudoVoigt::Exit():"<<this->GetName(),5)
1310  return;
1311  }
1312  if("Par"==tag.GetName())
1313  {
1314  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
1315  {
1316  if("Name"==tag.GetAttributeName(i))
1317  {
1318  this->GetPar(tag.GetAttributeValue(i)).XMLInput(is,tag);
1319  }
1320  }
1321  continue;
1322  }
1323  if("Option"==tag.GetName())
1324  {
1325  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
1326  if("Name"==tag.GetAttributeName(i))
1327  mOptionRegistry.GetObj(tag.GetAttributeValue(i)).XMLInput(is,tag);
1328  continue;
1329  }
1330  }
1331 }
1332 
1333 void ReflectionProfileDoubleExponentialPseudoVoigt::SetUnitCell(const UnitCell &cell)
1334 {
1335  VFN_DEBUG_MESSAGE("ReflectionProfileDoubleExponentialPseudoVoigt::SetUnitCell()",10)
1336  mpCell=&cell;
1337 }
1338 
1339 void ReflectionProfileDoubleExponentialPseudoVoigt
1340  ::InitParameters()
1341 {
1342  {
1343  RefinablePar tmp("Alpha0",&mInstrumentAlpha0,0,1e6,gpRefParTypeScattDataProfile,
1344  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1345  tmp.AssignClock(mClockMaster);
1346  tmp.SetDerivStep(1e-4);
1347  this->AddPar(tmp);
1348  }
1349  {
1350  RefinablePar tmp("Alpha1",&mInstrumentAlpha1,0,1e6,gpRefParTypeScattDataProfile,
1351  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1352  tmp.AssignClock(mClockMaster);
1353  tmp.SetDerivStep(1e-6);
1354  this->AddPar(tmp);
1355  }
1356  {
1357  RefinablePar tmp("Beta0",&mInstrumentBeta0,0,1e6,gpRefParTypeScattDataProfile,
1358  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1359  tmp.AssignClock(mClockMaster);
1360  tmp.SetDerivStep(1e-6);
1361  this->AddPar(tmp);
1362  }
1363  {
1364  RefinablePar tmp("Beta1",&mInstrumentBeta1,0,1e6,gpRefParTypeScattDataProfile,
1365  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1366  tmp.AssignClock(mClockMaster);
1367  tmp.SetDerivStep(1e-6);
1368  this->AddPar(tmp);
1369  }
1370  {
1371  RefinablePar tmp("GaussianSigma0",&mGaussianSigma0,0,1e6,gpRefParTypeScattDataProfileWidth,
1372  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1373  tmp.AssignClock(mClockMaster);
1374  tmp.SetDerivStep(1e-4);
1375  this->AddPar(tmp);
1376  }
1377  {
1378  RefinablePar tmp("GaussianSigma1",&mGaussianSigma1,0,1e6,gpRefParTypeScattDataProfileWidth,
1379  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1380  tmp.AssignClock(mClockMaster);
1381  tmp.SetDerivStep(1e-4);
1382  this->AddPar(tmp);
1383  }
1384  {
1385  RefinablePar tmp("GaussianSigma2",&mGaussianSigma2,0,1e6,gpRefParTypeScattDataProfileWidth,
1386  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1387  tmp.AssignClock(mClockMaster);
1388  tmp.SetDerivStep(1e-4);
1389  this->AddPar(tmp);
1390  }
1391  {
1392  RefinablePar tmp("LorentzianGamma0",&mLorentzianGamma0,0,1e6,gpRefParTypeScattDataProfileWidth,
1393  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1394  tmp.AssignClock(mClockMaster);
1395  tmp.SetDerivStep(1e-4);
1396  this->AddPar(tmp);
1397  }
1398  {
1399  RefinablePar tmp("LorentzianGamma1",&mLorentzianGamma1,0,1e6,gpRefParTypeScattDataProfileWidth,
1400  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1401  tmp.AssignClock(mClockMaster);
1402  tmp.SetDerivStep(1e-4);
1403  this->AddPar(tmp);
1404  }
1405  {
1406  RefinablePar tmp("LorentzianGamma2",&mLorentzianGamma2,0,1e6,gpRefParTypeScattDataProfileWidth,
1407  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1408  tmp.AssignClock(mClockMaster);
1409  tmp.SetDerivStep(1e-4);
1410  this->AddPar(tmp);
1411  }
1412 }
1413 
1414 #ifdef __WX__CRYST__
1415 WXCrystObjBasic* ReflectionProfileDoubleExponentialPseudoVoigt::WXCreate(wxWindow* parent)
1416 {
1417  if(mpWXCrystObj==0)
1418  mpWXCrystObj=new WXProfileDoubleExponentialPseudoVoigt(parent,this);
1419  return mpWXCrystObj;
1420 }
1421 #endif
1422 
1423 //######################################################################
1424 // Basic PROFILE FUNCTIONS
1425 //######################################################################
1426 
1427 CrystVector_REAL PowderProfileGauss (const CrystVector_REAL ttheta,const REAL fw,
1428  const REAL center, const REAL asym)
1429 {
1430  TAU_PROFILE("PowderProfileGauss()","Vector (Vector,REAL)",TAU_DEFAULT);
1431  REAL fwhm=fw;
1432  if(fwhm<=0) fwhm=1e-6;
1433  const long nbPoints=ttheta.numElements();
1434  CrystVector_REAL result(nbPoints);
1435  result=ttheta;
1436  result+= -center;
1437  result *= result;
1438  REAL *p;
1439  if(false)// fabs(asym-1.) < 1e-5)
1440  {
1441  //reference: IUCr Monographs on Crystallo 5 - The Rietveld Method (ed RA Young)
1442  result *= -4.*log(2.)/fwhm/fwhm;
1443  }
1444  else
1445  { // Adapted from Toraya J. Appl. Cryst 23(1990),485-491
1446  const REAL c1= -(1.+asym)/asym*(1.+asym)/asym*log(2.)/fwhm/fwhm;
1447  const REAL c2= -(1.+asym) *(1.+asym) *log(2.)/fwhm/fwhm;
1448  long i;
1449  p=result.data();
1450  const REAL *pt=ttheta.data();
1451  for(i=0;i<nbPoints;i++){ *p++ *= c1;if(*pt++>center) break;}
1452  i++;
1453  for( ;i<nbPoints;i++) *p++ *= c2;
1454  }
1455  p=result.data();
1456  #ifdef _MSC_VER
1457  // Bug from Hell (in MSVC++) !
1458  // The *last* point ends up sometimes with an arbitrary large value...
1459  for(long i=0;i<nbPoints;i++) { *p = pow((float)2.71828182846,(float)*p) ; p++ ;}
1460  #else
1461  long i=nbPoints;
1462  for(;i>3;i-=4)
1463  {
1464  #ifdef HAVE_SSE_MATHFUN
1465  v4sf x=_mm_loadu_ps(p);
1466  _mm_storeu_ps(p,exp_ps(x));
1467  p+=4;
1468  #else
1469  for(unsigned int j=0;j<4;++j)
1470  {// Fixed-length loop enables vectorization
1471  *p = exp(*p) ;
1472  p++ ;
1473  }
1474  #endif
1475  }
1476  for(;i>0;i--) { *p = exp(*p) ; p++ ;}
1477  #endif
1478 
1479 #if 0
1480  #if 1 //def _MSC_VER
1481  // Bug from Hell (in MSVC++) !
1482  // The *last* point ends up sometimes with an arbitrary large value...
1483  long i=0;
1484  for(;i<nbPoints;i+=4)
1485  for(unsigned int j=0;j<4;++j)
1486  {// Fixed-length loop enables vectorization
1487  *p = pow((float)2.71828182846,(float)*p) ;
1488  p++ ;
1489  }
1490  #else
1491  long i=0;
1492  for(;i<nbPoints;i+=1)
1493  {
1494  //for(unsigned int j=0;j<4;++j)
1495  {// Fixed-length loop enables vectorization
1496  *p = exp(*p) ;
1497  p++ ;
1498  }
1499  }
1500  #endif
1501 #endif
1502  result *= 2. / fwhm * sqrt(log(2.)/M_PI);
1503  return result;
1504 }
1505 
1506 CrystVector_REAL PowderProfileLorentz(const CrystVector_REAL ttheta,const REAL fw,
1507  const REAL center, const REAL asym)
1508 {
1509  TAU_PROFILE("PowderProfileLorentz()","Vector (Vector,REAL)",TAU_DEFAULT);
1510  REAL fwhm=fw;
1511  if(fwhm<=0) fwhm=1e-6;
1512  const long nbPoints=ttheta.numElements();
1513  CrystVector_REAL result(nbPoints);
1514  result=ttheta;
1515  result+= -center;
1516  result *= result;
1517  REAL *p;
1518  if(false)// fabs(asym-1.) < 1e-5)
1519  {
1520  //reference: IUCr Monographs on Crystallo 5 - The Rietveld Method (ed RA Young)
1521  result *= 4./fwhm/fwhm;
1522  }
1523  else
1524  { // Adapted from Toraya J. Appl. Cryst 23(1990),485-491
1525  const REAL c1= (1+asym)/asym*(1+asym)/asym/fwhm/fwhm;
1526  const REAL c2= (1+asym) *(1+asym) /fwhm/fwhm;
1527  long i;
1528  p=result.data();
1529  const REAL *pt=ttheta.data();
1530  for(i=0;i<nbPoints;i++){ *p++ *= c1;if(*pt++>center) break;}
1531  i++;
1532  for( ;i<nbPoints;i++) *p++ *= c2 ;
1533  }
1534  p=result.data();
1535  result += 1. ;
1536  for(long i=0;i<nbPoints;i++) { *p = 1/(*p) ; p++ ;}
1537  result *= 2./M_PI/fwhm;
1538  return result;
1539 }
1540 
1541 CrystVector_REAL AsymmetryBerarBaldinozzi(const CrystVector_REAL x,
1542  const REAL fw, const REAL center,
1543  const REAL a0, const REAL a1,
1544  const REAL b0, const REAL b1)
1545 {
1546  TAU_PROFILE("AsymmetryBerarBaldinozzi()","Vector (Vector,REAL)",TAU_DEFAULT);
1547  REAL fwhm=fw;
1548  if(fwhm<=0) fwhm=1e-6;
1549  const long nbPoints=x.numElements();
1550  CrystVector_REAL result(nbPoints);
1551  result=x;
1552  result+= -center;
1553  result *= 1/fwhm;
1554  REAL *p=result.data();
1555  const REAL a=a0/tan(center/2)+a1/tan(center);
1556  const REAL b=b0/tan(center/2)+b1/tan(center);
1557  for(long i=0;i<nbPoints;i++)
1558  {
1559  *p = 1+*p * exp(-*p * *p)*(2*a+b*(8* *p * *p-12));
1560  p++ ;
1561  }
1562  return result;
1563 }
1564 /*
1565 from python:
1566 E1(1)= 0.219383934396 (0.219383934396+0j)
1567 E1(1j)= (-0.337403922901-0.624713256428j)
1568 E1(1+1j)= (0.000281624451981-0.179324535039j)
1569 E1(100+1j)= (1.95936883899e-46-3.11904399563e-46j)
1570 E1(10+20j)= (-1.20141500252e-06-1.58298052926e-06j)
1571 
1572 this code (REAL=float)
1573 CE1(1.000000000000+0.000000000000j) = 0.219383955002+0.000000000000j
1574 CE1(0.000000000000+1.000000000000j) = -0.337403953075+-0.624713361263j
1575 CE1(1.000000000000+1.000000000000j) = 0.000281602144+-0.179324567318j
1576 CE1(100.000000000000+1.000000000000j) = 0.000000000000+-0.000000000000j
1577 CE1(10.000000000000+20.000000000000j) = -0.000001201415+-0.000001582981j
1578  {
1579  complex<REAL>z(1.0,0.0);
1580  complex<REAL>ce1=ExponentialIntegral1(z);
1581  cout<<"CE1("<<z.real()<<"+"<<z.imag()<<"j) = "<<ce1.real()<<"+"<<ce1.imag()<<"j"<<endl;
1582  }
1583  {
1584  complex<REAL>z(0.0,1.0);
1585  complex<REAL>ce1=ExponentialIntegral1(z);
1586  cout<<"CE1("<<z.real()<<"+"<<z.imag()<<"j) = "<<ce1.real()<<"+"<<ce1.imag()<<"j"<<endl;
1587  }
1588  {
1589  complex<REAL>z(1.0,1.0);
1590  complex<REAL>ce1=ExponentialIntegral1(z);
1591  cout<<"CE1("<<z.real()<<"+"<<z.imag()<<"j) = "<<ce1.real()<<"+"<<ce1.imag()<<"j"<<endl;
1592  }
1593  {
1594  complex<REAL>z(100.0,1.0);
1595  complex<REAL>ce1=ExponentialIntegral1(z);
1596  cout<<"CE1("<<z.real()<<"+"<<z.imag()<<"j) = "<<ce1.real()<<"+"<<ce1.imag()<<"j"<<endl;
1597  }
1598  {
1599  complex<REAL>z(10.0,20.0);
1600  complex<REAL>ce1=ExponentialIntegral1(z);
1601  cout<<"CE1("<<z.real()<<"+"<<z.imag()<<"j) = "<<ce1.real()<<"+"<<ce1.imag()<<"j"<<endl;
1602  }
1603  exit(0);
1604 */
1605 
1606 template <class T>std::complex<T>ExponentialIntegral1(const complex<T> z)
1607 {
1608  return exp(-z)*ExponentialIntegral1_ExpZ(z);
1609 }
1610 
1611 template <class T>std::complex<T>ExponentialIntegral1_ExpZ(const complex<T> z)
1612 {
1613  const T zr=z.real();
1614  const T zn=abs(z);
1615  complex<T> ce1;
1616  if(zn==0.0) return 1e100;// Should return an error ? std::numeric_limits::quiet_NaN() ?
1617  if((zn<10.0)||((zr<0.0)&&(zn<20.0)))// Abramowitz & Stegun 5.1.11
1618  {
1619  ce1=complex<T>(1,0);
1620  complex<T> y(1,0);
1621  for(unsigned int i=1;i<=150;i++)
1622  {
1623  y=-y*(T)i*z / (T)((i+1)*(i+1));
1624  ce1+=y;
1625  if(abs(y)<=abs(ce1)*1e-15) break;
1626  }
1627  static const T EulerMascheroni=0.5772156649015328606065120900;
1628  return exp(z)*(z*ce1-EulerMascheroni-log(z));// Euler-Mascheroni constant
1629  }
1630  else// Abramowitz & Stegun 5.1.51
1631  {
1632  if(zn>500) return 1.0/z;
1633  complex<T> y(0.0,0.0);
1634  for(unsigned int i=120;i>=1;i--) y=(T)i/((T)1+(T)i/(z+y));
1635  ce1/=(z+y);
1636  if((zr<0)&&(z.imag()==0)) ce1 -= complex<T>(0.0,M_PI)*exp(z);
1637  return ce1;
1638  }
1639 }
1640 
1641 }
ObjCryst
The namespace which includes all objects (crystallographic and algorithmic) in ObjCryst++.
Definition: doc-main.h:25
ObjCryst::gpRefParTypeScattDataProfileWidth
const RefParType * gpRefParTypeScattDataProfileWidth
Type for reflection profile width.
Definition: ScatteringData.cpp:59
ObjCryst::gReflectionProfileRegistry
ObjRegistry< ReflectionProfile > gReflectionProfileRegistry("List of all ReflectionProfile types")
Global registry for all ReflectionProfile objects.
Definition: ReflectionProfile.h:311
ObjCryst::gpRefParTypeScattDataProfile
const RefParType * gpRefParTypeScattDataProfile
Type for reflection profile.
Definition: ScatteringData.cpp:57
ObjCryst::ExponentialIntegral1_ExpZ
std::complex< T > ExponentialIntegral1_ExpZ(const complex< T > z)
E1(z)*exp(z)
Definition: ReflectionProfile.cpp:1611
ObjCryst::gpRefParTypeScattDataProfileType
const RefParType * gpRefParTypeScattDataProfileType
Type for reflection profiles type (e.g. gaussian/lorentzian mix)
Definition: ScatteringData.cpp:58
ObjCryst::AsymmetryBerarBaldinozzi
CrystVector_REAL AsymmetryBerarBaldinozzi(const CrystVector_REAL x, const REAL fw, const REAL center, const REAL a0, const REAL a1, const REAL b0, const REAL b1)
Asymmetry function [Ref J. Appl. Cryst 26 (1993), 128-129.
Definition: ReflectionProfile.cpp:1541
ObjCryst::ExponentialIntegral1
std::complex< T > ExponentialIntegral1(const complex< T > z)
Complex exponential integral E1(z) (Abramowitz & Stegun, chap.
Definition: ReflectionProfile.cpp:1606
ObjCryst::gpRefParTypeScattDataProfileAsym
const RefParType * gpRefParTypeScattDataProfileAsym
Type for reflection profile asymmetry.
Definition: ScatteringData.cpp:60
ObjCryst::ReflectionProfile
Abstract base class for reflection profiles.
Definition: ReflectionProfile.h:65
ObjCryst::ReflectionProfilePseudoVoigt
Pseudo-Voigt reflection profile.
Definition: ReflectionProfile.h:101
ObjCryst::ReflectionProfilePseudoVoigt::XMLOutput
virtual void XMLOutput(ostream &os, int indent=0) const
Output to stream in well-formed XML.
Definition: ReflectionProfile.cpp:324
ObjCryst::ReflectionProfilePseudoVoigt::GetClassName
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
Definition: ReflectionProfile.cpp:131
ObjCryst::ReflectionProfilePseudoVoigt::InitParameters
void InitParameters()
Initialize parameters.
Definition: ReflectionProfile.cpp:231
ObjCryst::ReflectionProfilePseudoVoigt::SetProfilePar
void SetProfilePar(const REAL fwhmCagliotiW, const REAL fwhmCagliotiU=0, const REAL fwhmCagliotiV=0, const REAL eta0=0.5, const REAL eta1=0.)
Set reflection profile parameters.
Definition: ReflectionProfile.cpp:173
ObjCryst::ReflectionProfilePseudoVoigt::mCagliotiU
REAL mCagliotiU
FWHM parameters, following Caglioti's law.
Definition: ReflectionProfile.h:132
ObjCryst::ReflectionProfilePseudoVoigt::mAsym0
REAL mAsym0
Asymmetry parameters, following the analytical function for asymmetric pseudo-voigt given by Toraya i...
Definition: ReflectionProfile.h:153
ObjCryst::ReflectionProfilePseudoVoigt::GetProfile
CrystVector_REAL GetProfile(const CrystVector_REAL &x, const REAL xcenter, const REAL h, const REAL k, const REAL l) const
Get the reflection profile.
Definition: ReflectionProfile.cpp:137
ObjCryst::ReflectionProfilePseudoVoigt::GetFullProfileWidth
virtual REAL GetFullProfileWidth(const REAL relativeIntensity, const REAL xcenter, const REAL h, const REAL k, const REAL l)
Get the (approximate) full profile width at a given percentage of the profile maximum (e....
Definition: ReflectionProfile.cpp:188
ObjCryst::ReflectionProfilePseudoVoigt::IsAnisotropic
bool IsAnisotropic() const
Is the profile anisotropic ?
Definition: ReflectionProfile.cpp:187
ObjCryst::ReflectionProfilePseudoVoigt::mPseudoVoigtEta0
REAL mPseudoVoigtEta0
Pseudo-Voigt mixing parameter : eta=eta0 +2*theta*eta1 eta=1 -> pure Lorentzian ; eta=0 -> pure Gauss...
Definition: ReflectionProfile.h:135
ObjCryst::ReflectionProfilePseudoVoigt::mAsymBerarBaldinozziA0
REAL mAsymBerarBaldinozziA0
Asymmetry parameters, following the Bérar & Baldinozzi approach ( Bérar & baldinozzi,...
Definition: ReflectionProfile.h:141
ObjCryst::ReflectionProfilePseudoVoigt::XMLInput
virtual void XMLInput(istream &is, const XMLCrystTag &tag)
Input From stream.
Definition: ReflectionProfile.cpp:374
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::GetFullProfileWidth
virtual REAL GetFullProfileWidth(const REAL relativeIntensity, const REAL xcenter, const REAL h, const REAL k, const REAL l)
Get the (approximate) full profile width at a given percentage of the profile maximum (e....
Definition: ReflectionProfile.cpp:600
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::mPseudoVoigtEta0
REAL mPseudoVoigtEta0
Pseudo-Voigt mixing parameter : eta=1 -> pure Lorentzian ; eta=0 -> pure Gaussian.
Definition: ReflectionProfile.h:231
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::GetClassName
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
Definition: ReflectionProfile.cpp:519
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::GetProfile
CrystVector_REAL GetProfile(const CrystVector_REAL &x, const REAL xcenter, const REAL h, const REAL k, const REAL l) const
Get the reflection profile.
Definition: ReflectionProfile.cpp:525
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::XMLInput
virtual void XMLInput(istream &is, const XMLCrystTag &tag)
Input From stream.
Definition: ReflectionProfile.cpp:716
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::XMLOutput
virtual void XMLOutput(ostream &os, int indent=0) const
Output to stream in well-formed XML.
Definition: ReflectionProfile.cpp:651
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::SetProfilePar
void SetProfilePar(const REAL fwhmCagliotiW, const REAL fwhmCagliotiU=0, const REAL fwhmCagliotiV=0, const REAL fwhmGaussP=0, const REAL fwhmLorentzX=0, const REAL fwhmLorentzY=0, const REAL fwhmLorentzGammaHH=0, const REAL fwhmLorentzGammaKK=0, const REAL fwhmLorentzGammaLL=0, const REAL fwhmLorentzGammaHK=0, const REAL fwhmLorentzGammaHL=0, const REAL fwhmLorentzGammaKL=0, const REAL pseudoVoigtEta0=0, const REAL pseudoVoigtEta1=0, const REAL asymA0=0, const REAL asymA1=0, const REAL asymA2=0)
Set reflection profile parameters.
Definition: ReflectionProfile.cpp:562
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::mAsym0
REAL mAsym0
Asymmetry parameters, following the analytical function for asymmetric pseudo-voigt given by Toraya i...
Definition: ReflectionProfile.h:242
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::IsAnisotropic
bool IsAnisotropic() const
Is the profile anisotropic ?
Definition: ReflectionProfile.cpp:646
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt
Double-Exponential Pseudo-Voigt profile for TOF.
Definition: ReflectionProfile.h:254
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::SetProfilePar
void SetProfilePar(const REAL instrumentAlpha0, const REAL instrumentAlpha1, const REAL instrumentBeta0, const REAL instrumentBeta1, const REAL gaussianSigma0, const REAL gaussianSigma1, const REAL gaussianSigma2, const REAL lorentzianGamma0, const REAL lorentzianGamma1, const REAL lorentzianGamma2)
Set reflection profile parameters.
Definition: ReflectionProfile.cpp:1164
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::IsAnisotropic
bool IsAnisotropic() const
Is the profile anisotropic ?
Definition: ReflectionProfile.cpp:1247
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::XMLOutput
virtual void XMLOutput(ostream &os, int indent=0) const
Output to stream in well-formed XML.
Definition: ReflectionProfile.cpp:1250
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::GetProfile
CrystVector_REAL GetProfile(const CrystVector_REAL &x, const REAL xcenter, const REAL h, const REAL k, const REAL l) const
Get the reflection profile.
Definition: ReflectionProfile.cpp:1068
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::XMLInput
virtual void XMLInput(istream &is, const XMLCrystTag &tag)
Input From stream.
Definition: ReflectionProfile.cpp:1296
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::SetUnitCell
void SetUnitCell(const UnitCell &cell)
Set unit cell.
Definition: ReflectionProfile.cpp:1333
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::GetFullProfileWidth
virtual REAL GetFullProfileWidth(const REAL relativeIntensity, const REAL xcenter, const REAL h, const REAL k, const REAL l)
Get the (approximate) full profile width at a given percentage of the profile maximum (e....
Definition: ReflectionProfile.cpp:1189
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::ReflectionProfileDoubleExponentialPseudoVoigt
ReflectionProfileDoubleExponentialPseudoVoigt()
Constructor, without unit cell.
Definition: ReflectionProfile.cpp:988
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::GetClassName
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
Definition: ReflectionProfile.cpp:1061
ObjCryst::UnitCell
Unit Cell class: Unit cell with spacegroup information.
Definition: UnitCell.h:72
ObjCryst::UnitCell::MillerToOrthonormalCoords
void MillerToOrthonormalCoords(REAL &x, REAL &y, REAL &z) const
Get Miller H,K, L indices from orthonormal coordinates in reciprocal space.
Definition: UnitCell.cpp:286
ObjCryst::XMLCrystTag
class to input or output a well-formatted xml beginning or ending tag.
Definition: RefinableObj/IO.h:73
ObjCryst::RefinableObjClock::Click
void Click()
Record an event for this clock (generally, the 'time' an object has been modified,...
Definition: RefinableObj.cpp:131
ObjCryst::RefinablePar
Generic class for parameters of refinable objects.
Definition: RefinableObj.h:225
ObjCryst::RefinablePar::XMLOutput
void XMLOutput(ostream &os, const string &name, int indent=0) const
XMLOutput to stream in well-formed XML.
Definition: RefinableObj/IO.cpp:184
ObjCryst::RefinablePar::SetDerivStep
void SetDerivStep(const REAL)
Fixed step to use to compute numerical derivative.
Definition: RefinableObj.cpp:590
ObjCryst::RefinablePar::AssignClock
void AssignClock(RefinableObjClock &clock)
Definition: RefinableObj.cpp:664
ObjCryst::RefinablePar::XMLInput
void XMLInput(istream &is, const XMLCrystTag &tag)
XMLInput From stream.
Definition: RefinableObj/IO.cpp:229
ObjCryst::RefinableObj::AddPar
void AddPar(const RefinablePar &newRefPar)
Add a refinable parameter.
Definition: RefinableObj.cpp:1455
ObjCryst::RefinableObj::SetName
virtual void SetName(const string &name)
Name of the object.
Definition: RefinableObj.cpp:1289
ObjCryst::RefinableObj::GetPar
RefinablePar & GetPar(const long i)
Access all parameters in the order they were inputted.
Definition: RefinableObj.cpp:1391
ObjCryst::RefinableObj::GetName
virtual const string & GetName() const
Name of the object.
Definition: RefinableObj.cpp:1288
ObjCryst::RefinableObj::mOptionRegistry
ObjRegistry< RefObjOpt > mOptionRegistry
List of options for this object.
Definition: RefinableObj.h:1303
ObjCryst::RefinableObj::UpdateDisplay
virtual void UpdateDisplay() const
If there is an interface, this should be automatically be called each time there is a 'new,...
Definition: RefinableObj.cpp:1986
ObjCryst::RefinableObj::mClockMaster
RefinableObjClock mClockMaster
Master clock, which is changed whenever the object has been altered.
Definition: RefinableObj.h:1330
FormatVertVector
output one or several vectors as (a) column(s):
Definition: VFNStreamFormat.h:97
ObjCryst::WXCrystObjBasic
Abstract base class for all objects in wxCryst.
Definition: wxCryst.h:128
ObjCryst::WXProfilePseudoVoigt
Class to display a Powder Pattern Pseudo-Voigt Profile.
Definition: wxPowderPattern.h:287
ObjCryst::WXProfilePseudoVoigtAnisotropic
Class to display a Powder Pattern Pseudo-Voigt Profile with Anisotropic broadening.
Definition: wxPowderPattern.h:300
ObjCryst::WXProfileDoubleExponentialPseudoVoigt
Class to display a Powder Pattern Pseudo-Voigt Profile.
Definition: wxPowderPattern.h:313