FOX/ObjCryst++  1.10.X (development)
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),10)
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  VFN_DEBUG_EXIT("ReflectionProfileDoubleExponentialPseudoVoigt::GetFullProfileWidth()",5)
1239 }
1240 
1241 bool ReflectionProfileDoubleExponentialPseudoVoigt
1242  ::IsAnisotropic()const{return false;}
1243 
1244 void ReflectionProfileDoubleExponentialPseudoVoigt
1245  ::XMLOutput(ostream &os,int indent)const
1246 {
1247  VFN_DEBUG_ENTRY("ReflectionProfileDoubleExponentialPseudoVoigt::XMLOutput():"<<this->GetName(),5)
1248  for(int i=0;i<indent;i++) os << " " ;
1249  XMLCrystTag tag("ReflectionProfileDoubleExponentialPseudoVoigt");
1250  os <<tag<<endl;
1251  indent++;
1252 
1253  this->GetPar(&mInstrumentAlpha0).XMLOutput(os,"Alpha0",indent);
1254  os <<endl;
1255 
1256  this->GetPar(&mInstrumentAlpha1).XMLOutput(os,"Alpha1",indent);
1257  os <<endl;
1258 
1259  this->GetPar(&mInstrumentBeta0).XMLOutput(os,"Beta0",indent);
1260  os <<endl;
1261 
1262  this->GetPar(&mInstrumentBeta1).XMLOutput(os,"Beta1",indent);
1263  os <<endl;
1264 
1265  this->GetPar(&mGaussianSigma0).XMLOutput(os,"GaussianSigma0",indent);
1266  os <<endl;
1267 
1268  this->GetPar(&mGaussianSigma1).XMLOutput(os,"GaussianSigma1",indent);
1269  os <<endl;
1270 
1271  this->GetPar(&mGaussianSigma2).XMLOutput(os,"GaussianSigma2",indent);
1272  os <<endl;
1273 
1274  this->GetPar(&mLorentzianGamma0).XMLOutput(os,"LorentzianGamma0",indent);
1275  os <<endl;
1276 
1277  this->GetPar(&mLorentzianGamma1).XMLOutput(os,"LorentzianGamma1",indent);
1278  os <<endl;
1279 
1280  this->GetPar(&mLorentzianGamma2).XMLOutput(os,"LorentzianGamma2",indent);
1281  os <<endl;
1282 
1283  indent--;
1284  tag.SetIsEndTag(true);
1285  for(int i=0;i<indent;i++) os << " " ;
1286  os <<tag<<endl;
1287  VFN_DEBUG_EXIT("ReflectionProfileDoubleExponentialPseudoVoigt::XMLOutput():"<<this->GetName(),5)
1288 }
1289 
1290 void ReflectionProfileDoubleExponentialPseudoVoigt
1291  ::XMLInput(istream &is,const XMLCrystTag &tagg)
1292 {
1293  VFN_DEBUG_ENTRY("ReflectionProfileDoubleExponentialPseudoVoigt::XMLInput():"<<this->GetName(),5)
1294  for(unsigned int i=0;i<tagg.GetNbAttribute();i++)
1295  {
1296  if("Name"==tagg.GetAttributeName(i)) this->SetName(tagg.GetAttributeValue(i));
1297  }
1298  while(true)
1299  {
1300  XMLCrystTag tag(is);
1301  if(("ReflectionProfileDoubleExponentialPseudoVoigt"==tag.GetName())&&tag.IsEndTag())
1302  {
1303  this->UpdateDisplay();
1304  VFN_DEBUG_EXIT("ReflectionProfileDoubleExponentialPseudoVoigt::Exit():"<<this->GetName(),5)
1305  return;
1306  }
1307  if("Par"==tag.GetName())
1308  {
1309  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
1310  {
1311  if("Name"==tag.GetAttributeName(i))
1312  {
1313  this->GetPar(tag.GetAttributeValue(i)).XMLInput(is,tag);
1314  }
1315  }
1316  continue;
1317  }
1318  if("Option"==tag.GetName())
1319  {
1320  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
1321  if("Name"==tag.GetAttributeName(i))
1322  mOptionRegistry.GetObj(tag.GetAttributeValue(i)).XMLInput(is,tag);
1323  continue;
1324  }
1325  }
1326 }
1327 
1328 void ReflectionProfileDoubleExponentialPseudoVoigt::SetUnitCell(const UnitCell &cell)
1329 {
1330  VFN_DEBUG_MESSAGE("ReflectionProfileDoubleExponentialPseudoVoigt::SetUnitCell()",10)
1331  mpCell=&cell;
1332 }
1333 
1334 void ReflectionProfileDoubleExponentialPseudoVoigt
1335  ::InitParameters()
1336 {
1337  {
1338  RefinablePar tmp("Alpha0",&mInstrumentAlpha0,0,1e6,gpRefParTypeScattDataProfile,
1339  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1340  tmp.AssignClock(mClockMaster);
1341  tmp.SetDerivStep(1e-4);
1342  this->AddPar(tmp);
1343  }
1344  {
1345  RefinablePar tmp("Alpha1",&mInstrumentAlpha1,0,1e6,gpRefParTypeScattDataProfile,
1346  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1347  tmp.AssignClock(mClockMaster);
1348  tmp.SetDerivStep(1e-6);
1349  this->AddPar(tmp);
1350  }
1351  {
1352  RefinablePar tmp("Beta0",&mInstrumentBeta0,0,1e6,gpRefParTypeScattDataProfile,
1353  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1354  tmp.AssignClock(mClockMaster);
1355  tmp.SetDerivStep(1e-6);
1356  this->AddPar(tmp);
1357  }
1358  {
1359  RefinablePar tmp("Beta1",&mInstrumentBeta1,0,1e6,gpRefParTypeScattDataProfile,
1360  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1361  tmp.AssignClock(mClockMaster);
1362  tmp.SetDerivStep(1e-6);
1363  this->AddPar(tmp);
1364  }
1365  {
1366  RefinablePar tmp("GaussianSigma0",&mGaussianSigma0,0,1e6,gpRefParTypeScattDataProfileWidth,
1367  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1368  tmp.AssignClock(mClockMaster);
1369  tmp.SetDerivStep(1e-4);
1370  this->AddPar(tmp);
1371  }
1372  {
1373  RefinablePar tmp("GaussianSigma1",&mGaussianSigma1,0,1e6,gpRefParTypeScattDataProfileWidth,
1374  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1375  tmp.AssignClock(mClockMaster);
1376  tmp.SetDerivStep(1e-4);
1377  this->AddPar(tmp);
1378  }
1379  {
1380  RefinablePar tmp("GaussianSigma2",&mGaussianSigma2,0,1e6,gpRefParTypeScattDataProfileWidth,
1381  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1382  tmp.AssignClock(mClockMaster);
1383  tmp.SetDerivStep(1e-4);
1384  this->AddPar(tmp);
1385  }
1386  {
1387  RefinablePar tmp("LorentzianGamma0",&mLorentzianGamma0,0,1e6,gpRefParTypeScattDataProfileWidth,
1388  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1389  tmp.AssignClock(mClockMaster);
1390  tmp.SetDerivStep(1e-4);
1391  this->AddPar(tmp);
1392  }
1393  {
1394  RefinablePar tmp("LorentzianGamma1",&mLorentzianGamma1,0,1e6,gpRefParTypeScattDataProfileWidth,
1395  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1396  tmp.AssignClock(mClockMaster);
1397  tmp.SetDerivStep(1e-4);
1398  this->AddPar(tmp);
1399  }
1400  {
1401  RefinablePar tmp("LorentzianGamma2",&mLorentzianGamma2,0,1e6,gpRefParTypeScattDataProfileWidth,
1402  REFPAR_DERIV_STEP_ABSOLUTE,true,true,true,false);
1403  tmp.AssignClock(mClockMaster);
1404  tmp.SetDerivStep(1e-4);
1405  this->AddPar(tmp);
1406  }
1407 }
1408 
1409 #ifdef __WX__CRYST__
1410 WXCrystObjBasic* ReflectionProfileDoubleExponentialPseudoVoigt::WXCreate(wxWindow* parent)
1411 {
1412  if(mpWXCrystObj==0)
1413  mpWXCrystObj=new WXProfileDoubleExponentialPseudoVoigt(parent,this);
1414  return mpWXCrystObj;
1415 }
1416 #endif
1417 
1418 //######################################################################
1419 // Basic PROFILE FUNCTIONS
1420 //######################################################################
1421 
1422 CrystVector_REAL PowderProfileGauss (const CrystVector_REAL ttheta,const REAL fw,
1423  const REAL center, const REAL asym)
1424 {
1425  TAU_PROFILE("PowderProfileGauss()","Vector (Vector,REAL)",TAU_DEFAULT);
1426  REAL fwhm=fw;
1427  if(fwhm<=0) fwhm=1e-6;
1428  const long nbPoints=ttheta.numElements();
1429  CrystVector_REAL result(nbPoints);
1430  result=ttheta;
1431  result+= -center;
1432  result *= result;
1433  REAL *p;
1434  if(false)// fabs(asym-1.) < 1e-5)
1435  {
1436  //reference: IUCr Monographs on Crystallo 5 - The Rietveld Method (ed RA Young)
1437  result *= -4.*log(2.)/fwhm/fwhm;
1438  }
1439  else
1440  { // Adapted from Toraya J. Appl. Cryst 23(1990),485-491
1441  const REAL c1= -(1.+asym)/asym*(1.+asym)/asym*log(2.)/fwhm/fwhm;
1442  const REAL c2= -(1.+asym) *(1.+asym) *log(2.)/fwhm/fwhm;
1443  long i;
1444  p=result.data();
1445  const REAL *pt=ttheta.data();
1446  for(i=0;i<nbPoints;i++){ *p++ *= c1;if(*pt++>center) break;}
1447  i++;
1448  for( ;i<nbPoints;i++) *p++ *= c2;
1449  }
1450  p=result.data();
1451  #ifdef _MSC_VER
1452  // Bug from Hell (in MSVC++) !
1453  // The *last* point ends up sometimes with an arbitrary large value...
1454  for(long i=0;i<nbPoints;i++) { *p = pow((float)2.71828182846,(float)*p) ; p++ ;}
1455  #else
1456  long i=nbPoints;
1457  for(;i>3;i-=4)
1458  {
1459  #ifdef HAVE_SSE_MATHFUN
1460  v4sf x=_mm_loadu_ps(p);
1461  _mm_storeu_ps(p,exp_ps(x));
1462  p+=4;
1463  #else
1464  for(unsigned int j=0;j<4;++j)
1465  {// Fixed-length loop enables vectorization
1466  *p = exp(*p) ;
1467  p++ ;
1468  }
1469  #endif
1470  }
1471  for(;i>0;i--) { *p = exp(*p) ; p++ ;}
1472  #endif
1473 
1474 #if 0
1475  #if 1 //def _MSC_VER
1476  // Bug from Hell (in MSVC++) !
1477  // The *last* point ends up sometimes with an arbitrary large value...
1478  long i=0;
1479  for(;i<nbPoints;i+=4)
1480  for(unsigned int j=0;j<4;++j)
1481  {// Fixed-length loop enables vectorization
1482  *p = pow((float)2.71828182846,(float)*p) ;
1483  p++ ;
1484  }
1485  #else
1486  long i=0;
1487  for(;i<nbPoints;i+=1)
1488  {
1489  //for(unsigned int j=0;j<4;++j)
1490  {// Fixed-length loop enables vectorization
1491  *p = exp(*p) ;
1492  p++ ;
1493  }
1494  }
1495  #endif
1496 #endif
1497  result *= 2. / fwhm * sqrt(log(2.)/M_PI);
1498  return result;
1499 }
1500 
1501 CrystVector_REAL PowderProfileLorentz(const CrystVector_REAL ttheta,const REAL fw,
1502  const REAL center, const REAL asym)
1503 {
1504  TAU_PROFILE("PowderProfileLorentz()","Vector (Vector,REAL)",TAU_DEFAULT);
1505  REAL fwhm=fw;
1506  if(fwhm<=0) fwhm=1e-6;
1507  const long nbPoints=ttheta.numElements();
1508  CrystVector_REAL result(nbPoints);
1509  result=ttheta;
1510  result+= -center;
1511  result *= result;
1512  REAL *p;
1513  if(false)// fabs(asym-1.) < 1e-5)
1514  {
1515  //reference: IUCr Monographs on Crystallo 5 - The Rietveld Method (ed RA Young)
1516  result *= 4./fwhm/fwhm;
1517  }
1518  else
1519  { // Adapted from Toraya J. Appl. Cryst 23(1990),485-491
1520  const REAL c1= (1+asym)/asym*(1+asym)/asym/fwhm/fwhm;
1521  const REAL c2= (1+asym) *(1+asym) /fwhm/fwhm;
1522  long i;
1523  p=result.data();
1524  const REAL *pt=ttheta.data();
1525  for(i=0;i<nbPoints;i++){ *p++ *= c1;if(*pt++>center) break;}
1526  i++;
1527  for( ;i<nbPoints;i++) *p++ *= c2 ;
1528  }
1529  p=result.data();
1530  result += 1. ;
1531  for(long i=0;i<nbPoints;i++) { *p = 1/(*p) ; p++ ;}
1532  result *= 2./M_PI/fwhm;
1533  return result;
1534 }
1535 
1536 CrystVector_REAL AsymmetryBerarBaldinozzi(const CrystVector_REAL x,
1537  const REAL fw, const REAL center,
1538  const REAL a0, const REAL a1,
1539  const REAL b0, const REAL b1)
1540 {
1541  TAU_PROFILE("AsymmetryBerarBaldinozzi()","Vector (Vector,REAL)",TAU_DEFAULT);
1542  REAL fwhm=fw;
1543  if(fwhm<=0) fwhm=1e-6;
1544  const long nbPoints=x.numElements();
1545  CrystVector_REAL result(nbPoints);
1546  result=x;
1547  result+= -center;
1548  result *= 1/fwhm;
1549  REAL *p=result.data();
1550  const REAL a=a0/tan(center/2)+a1/tan(center);
1551  const REAL b=b0/tan(center/2)+b1/tan(center);
1552  for(long i=0;i<nbPoints;i++)
1553  {
1554  *p = 1+*p * exp(-*p * *p)*(2*a+b*(8* *p * *p-12));
1555  p++ ;
1556  }
1557  return result;
1558 }
1559 /*
1560 from python:
1561 E1(1)= 0.219383934396 (0.219383934396+0j)
1562 E1(1j)= (-0.337403922901-0.624713256428j)
1563 E1(1+1j)= (0.000281624451981-0.179324535039j)
1564 E1(100+1j)= (1.95936883899e-46-3.11904399563e-46j)
1565 E1(10+20j)= (-1.20141500252e-06-1.58298052926e-06j)
1566 
1567 this code (REAL=float)
1568 CE1(1.000000000000+0.000000000000j) = 0.219383955002+0.000000000000j
1569 CE1(0.000000000000+1.000000000000j) = -0.337403953075+-0.624713361263j
1570 CE1(1.000000000000+1.000000000000j) = 0.000281602144+-0.179324567318j
1571 CE1(100.000000000000+1.000000000000j) = 0.000000000000+-0.000000000000j
1572 CE1(10.000000000000+20.000000000000j) = -0.000001201415+-0.000001582981j
1573  {
1574  complex<REAL>z(1.0,0.0);
1575  complex<REAL>ce1=ExponentialIntegral1(z);
1576  cout<<"CE1("<<z.real()<<"+"<<z.imag()<<"j) = "<<ce1.real()<<"+"<<ce1.imag()<<"j"<<endl;
1577  }
1578  {
1579  complex<REAL>z(0.0,1.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(1.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(100.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(10.0,20.0);
1595  complex<REAL>ce1=ExponentialIntegral1(z);
1596  cout<<"CE1("<<z.real()<<"+"<<z.imag()<<"j) = "<<ce1.real()<<"+"<<ce1.imag()<<"j"<<endl;
1597  }
1598  exit(0);
1599 */
1600 
1601 template <class T>std::complex<T>ExponentialIntegral1(const complex<T> z)
1602 {
1603  return exp(-z)*ExponentialIntegral1_ExpZ(z);
1604 }
1605 
1606 template <class T>std::complex<T>ExponentialIntegral1_ExpZ(const complex<T> z)
1607 {
1608  const T zr=z.real();
1609  const T zn=abs(z);
1610  complex<T> ce1;
1611  if(zn==0.0) return 1e100;// Should return an error ? std::numeric_limits::quiet_NaN() ?
1612  if((zn<10.0)||((zr<0.0)&&(zn<20.0)))// Abramowitz & Stegun 5.1.11
1613  {
1614  ce1=complex<T>(1,0);
1615  complex<T> y(1,0);
1616  for(unsigned int i=1;i<=150;i++)
1617  {
1618  y=-y*(T)i*z / (T)((i+1)*(i+1));
1619  ce1+=y;
1620  if(abs(y)<=abs(ce1)*1e-15) break;
1621  }
1622  static const T EulerMascheroni=0.5772156649015328606065120900;
1623  return exp(z)*(z*ce1-EulerMascheroni-log(z));// Euler-Mascheroni constant
1624  }
1625  else// Abramowitz & Stegun 5.1.51
1626  {
1627  if(zn>500) return 1.0/z;
1628  complex<T> y(0.0,0.0);
1629  for(unsigned int i=120;i>=1;i--) y=(T)i/((T)1+(T)i/(z+y));
1630  ce1/=(z+y);
1631  if((zr<0)&&(z.imag()==0)) ce1 -= complex<T>(0.0,M_PI)*exp(z);
1632  return ce1;
1633  }
1634 }
1635 
1636 }
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::GetClassName
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
Definition: ReflectionProfile.cpp:1057
ObjCryst::RefinablePar::SetDerivStep
void SetDerivStep(const REAL)
Fixed step to use to compute numerical derivative.
Definition: RefinableObj.cpp:589
ObjCryst::RefinableObj::AddPar
void AddPar(const RefinablePar &newRefPar)
Add a refinable parameter.
Definition: RefinableObj.cpp:1416
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::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:1917
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt
Double-Exponential Pseudo-Voigt profile for TOF.
Definition: ReflectionProfile.h:253
ObjCryst::ReflectionProfilePseudoVoigt::GetClassName
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
Definition: ReflectionProfile.cpp:131
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::IsAnisotropic
bool IsAnisotropic() const
Is the profile anisotropic ?
Definition: ReflectionProfile.cpp:642
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::ReflectionProfileDoubleExponentialPseudoVoigt::SetUnitCell
void SetUnitCell(const UnitCell &cell)
Set unit cell.
Definition: ReflectionProfile.cpp:1328
ObjCryst::ReflectionProfilePseudoVoigt::IsAnisotropic
bool IsAnisotropic() const
Is the profile anisotropic ?
Definition: ReflectionProfile.cpp:185
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::IsAnisotropic
bool IsAnisotropic() const
Is the profile anisotropic ?
Definition: ReflectionProfile.cpp:1242
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:284
ObjCryst::RefinableObjClock::Click
void Click()
Record an event for this clock (generally, the 'time' an object has been modified, or some computation has been made)
Definition: RefinableObj.cpp:130
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::XMLInput
virtual void XMLInput(istream &is, const XMLCrystTag &tag)
Input From stream.
Definition: ReflectionProfile.cpp:712
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::RefinablePar::XMLInput
void XMLInput(istream &is, const XMLCrystTag &tag)
XMLInput From stream.
Definition: RefinableObj/IO.cpp:229
ObjCryst::RefinableObj::GetPar
RefinablePar & GetPar(const long i)
Access all parameters in the order they were inputted.
Definition: RefinableObj.cpp:1352
ObjCryst::WXProfilePseudoVoigtAnisotropic
Class to display a Powder Pattern Pseudo-Voigt Profile with Anisotropic broadening.
Definition: wxPowderPattern.h:296
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::GetClassName
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
Definition: ReflectionProfile.cpp:517
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::XMLInput
virtual void XMLInput(istream &is, const XMLCrystTag &tag)
Input From stream.
Definition: ReflectionProfile.cpp:1291
ObjCryst::gpRefParTypeScattDataProfile
const RefParType * gpRefParTypeScattDataProfile
Type for reflection profile.
Definition: ScatteringData.cpp:57
FormatVertVector
output one or several vectors as (a) column(s):
Definition: VFNStreamFormat.h:96
ObjCryst::gpRefParTypeScattDataProfileAsym
const RefParType * gpRefParTypeScattDataProfileAsym
Type for reflection profile asymmetry.
Definition: ScatteringData.cpp:60
ObjCryst::ReflectionProfilePseudoVoigt::mAsymBerarBaldinozziA0
REAL mAsymBerarBaldinozziA0
Asymmetry parameters, following the Bérar & Baldinozzi approach ( Bérar & baldinozzi, J.
Definition: ReflectionProfile.h:141
ObjCryst::ReflectionProfilePseudoVoigt::InitParameters
void InitParameters()
Initialize parameters.
Definition: ReflectionProfile.cpp:229
ObjCryst::RefinableObj::mClockMaster
RefinableObjClock mClockMaster
Master clock, which is changed whenever the object has been altered.
Definition: RefinableObj.h:1291
ObjCryst::ReflectionProfilePseudoVoigt::mAsym0
REAL mAsym0
Asymmetry parameters, following the analytical function for asymmetric pseudo-voigt given by Toraya i...
Definition: ReflectionProfile.h:153
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::mPseudoVoigtEta0
REAL mPseudoVoigtEta0
Pseudo-Voigt mixing parameter : eta=1 -> pure Lorentzian ; eta=0 -> pure Gaussian.
Definition: ReflectionProfile.h:231
ObjCryst::WXCrystObjBasic
Abstract base class for all objects in wxCryst.
Definition: wxCryst.h:127
ObjCryst::ReflectionProfile
Abstract base class for reflection profiles.
Definition: ReflectionProfile.h:64
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::ReflectionProfilePseudoVoigt::XMLOutput
virtual void XMLOutput(ostream &os, int indent=0) const
Output to stream in well-formed XML.
Definition: ReflectionProfile.cpp:322
ObjCryst::RefinablePar::AssignClock
void AssignClock(RefinableObjClock &clock)
Definition: RefinableObj.cpp:663
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::gpRefParTypeScattDataProfileWidth
const RefParType * gpRefParTypeScattDataProfileWidth
Type for reflection profile width.
Definition: ScatteringData.cpp:59
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::XMLOutput
virtual void XMLOutput(ostream &os, int indent=0) const
Output to stream in well-formed XML.
Definition: ReflectionProfile.cpp:1245
ObjCryst::WXProfilePseudoVoigt
Class to display a Powder Pattern Pseudo-Voigt Profile.
Definition: wxPowderPattern.h:283
ObjCryst::gpRefParTypeScattDataProfileType
const RefParType * gpRefParTypeScattDataProfileType
Type for reflection profiles type (e.g. gaussian/lorentzian mix)
Definition: ScatteringData.cpp:58
ObjCryst::ExponentialIntegral1
std::complex< T > ExponentialIntegral1(const complex< T > z)
Complex exponential integral E1(z) (Abramowitz & Stegun, chap.
Definition: ReflectionProfile.cpp:1601
ObjCryst::ReflectionProfilePseudoVoigtAnisotropic::mAsym0
REAL mAsym0
Asymmetry parameters, following the analytical function for asymmetric pseudo-voigt given by Toraya i...
Definition: ReflectionProfile.h:242
ObjCryst::ReflectionProfilePseudoVoigt::XMLInput
virtual void XMLInput(istream &is, const XMLCrystTag &tag)
Input From stream.
Definition: ReflectionProfile.cpp:372
ObjCryst
The namespace which includes all objects (crystallographic and algorithmic) in ObjCryst++.
Definition: Atom.cpp:47
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::UnitCell
Unit Cell class: Unit cell with spacegroup information.
Definition: UnitCell.h:71
ObjCryst::RefinablePar
Generic class for parameters of refinable objects.
Definition: RefinableObj.h:223
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::ReflectionProfilePseudoVoigt::mCagliotiU
REAL mCagliotiU
FWHM parameters, following Caglioti's law.
Definition: ReflectionProfile.h:132
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::ReflectionProfilePseudoVoigtAnisotropic::XMLOutput
virtual void XMLOutput(ostream &os, int indent=0) const
Output to stream in well-formed XML.
Definition: ReflectionProfile.cpp:647
ObjCryst::RefinableObj::GetName
virtual const string & GetName() const
Name of the object.
Definition: RefinableObj.cpp:1249
ObjCryst::WXProfileDoubleExponentialPseudoVoigt
Class to display a Powder Pattern Pseudo-Voigt Profile.
Definition: wxPowderPattern.h:309
ObjCryst::ReflectionProfilePseudoVoigt
Pseudo-Voigt reflection profile.
Definition: ReflectionProfile.h:100
ObjCryst::XMLCrystTag
class to input or output a well-formatted xml beginning or ending tag.
Definition: RefinableObj/IO.h:72
ObjCryst::RefinableObj::mOptionRegistry
ObjRegistry< RefObjOpt > mOptionRegistry
List of options for this object.
Definition: RefinableObj.h:1264
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:1536
ObjCryst::ReflectionProfileDoubleExponentialPseudoVoigt::ReflectionProfileDoubleExponentialPseudoVoigt
ReflectionProfileDoubleExponentialPseudoVoigt()
Constructor, without unit cell.
Definition: ReflectionProfile.cpp:984
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::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::RefinableObj::SetName
virtual void SetName(const string &name)
Name of the object.
Definition: RefinableObj.cpp:1250
ObjCryst::ExponentialIntegral1_ExpZ
std::complex< T > ExponentialIntegral1_ExpZ(const complex< T > z)
E1(z)*exp(z)
Definition: ReflectionProfile.cpp:1606
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::gReflectionProfileRegistry
ObjRegistry< ReflectionProfile > gReflectionProfileRegistry("List of all ReflectionProfile types")
Global registry for all ReflectionProfile objects.
Definition: ReflectionProfile.h:311