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