PEARL Procedures  rev-distro-3.0.0-0-gfa24916-dirty
Igor procedures for the analysis of PEARL data
pearl-fitfuncs.ipf
Go to the documentation of this file.
1 #pragma TextEncoding = "UTF-8"
2 #pragma rtGlobals=3 // Use modern global access method and strict wave access.
3 #pragma IgorVersion = 6.2
4 #pragma ModuleName = PearlFitFuncs
5 #pragma version = 1.02
6 #include "mm-physconst"
7 
21 
26 
27 
28 //------------------------------------------------------------------------------
29 // Gaussian shapes
30 //------------------------------------------------------------------------------
31 
45 threadsafe function MultiGaussLinBG(w,x) : FitFunc
46  wave w
47  variable x
48 
49  variable np = numpnts(w)
50  variable ip
51 
52  variable v = w[0] + x * w[1]
53  for (ip = 2; ip < np; ip += 3)
54  v += w[ip] * exp( -( (x - w[ip+1]) / w[ip+2] )^2 )
55  endfor
56 
57  return v
58 end
59 
80 threadsafe function MultiGaussLinBG_AO(pw, yw, xw) : FitFunc
81  wave pw
82  wave yw
83  wave xw
84 
85  variable np = numpnts(pw)
86  variable ip
87 
88  yw = pw[0] + xw[p] * pw[1]
89  for (ip = 2; ip < np; ip += 3)
90  yw += pw[ip] * exp( -( (xw[p] - pw[ip+1]) / pw[ip+2] )^2 )
91  endfor
92 end
93 
118 threadsafe function DoubletGaussLinBG_AO(pw, yw, xw) : FitFunc
119  wave pw
120  wave yw
121  wave xw
122 
123  yw = pw[0] + xw[p] * pw[1]
124  yw += pw[2] * exp( -( (xw[p] - pw[4] - pw[5] /2) / pw[6] )^2 )
125  yw += pw[2] * pw[3] * exp( -( (xw[p] - pw[4] + pw[5] /2) / pw[6] / pw[7] )^2 )
126 end
127 
157 threadsafe function DblDoubletGaussLinBG_AO(pw, yw, xw) : FitFunc
158  wave pw
159  wave yw
160  wave xw
161 
162  yw = pw[0] + xw[p] * pw[1]
163  // peak 1
164  yw += pw[2] * exp( -( (xw[p] - pw[6]) / pw[9] )^2 )
165  // peak 2
166  yw += pw[3] * exp( -( (xw[p] - pw[6] - pw[8]) / pw[10] )^2 )
167  // peak 3
168  yw += pw[2] * pw[4] * exp( -( (xw[p] - pw[6] - pw[7]) / pw[9] )^2 )
169  // peak 4
170  yw += pw[3] * pw[5] * exp( -( (xw[p] - pw[6] - pw[7] - pw[8]) / pw[10] )^2 )
171 end
172 
173 //------------------------------------------------------------------------------
174 // Voigt shapes
175 //------------------------------------------------------------------------------
176 
190 threadsafe function MultiVoigtLinBG(w,x) : FitFunc
191  wave w
192  variable x
193 
194  variable np = numpnts(w)
195  variable ip
196 
197  variable v = w[0] + x * w[1]
198  for (ip = 2; ip < np; ip += 4)
199  v += w[ip] * VoigtFunc((x - w[ip+1]) / w[ip+2], w[ip+3])
200  endfor
201 
202  return v
203 end
204 
224 threadsafe function MultiVoigtLinBG_AO(pw, yw, xw) : FitFunc
225  wave pw
226  wave yw
227  wave xw
228 
229  variable np = numpnts(pw)
230  variable ip
231 
232  yw = pw[0] + xw[p] * pw[1]
233  for (ip = 2; ip < np; ip += 4)
234  yw += pw[ip] * VoigtFunc((xw[p] - pw[ip+1]) / pw[ip+2], pw[ip+3])
235  endfor
236 end
237 
238 
239 //------------------------------------------------------------------------------
240 // Doniach-Sunjic shapes
241 //------------------------------------------------------------------------------
242 
253 threadsafe function DoniachSunjic(x, amp, pos, sing, fwhm)
254  variable x
255  variable amp
256  variable pos
257  variable sing
258  variable fwhm
259 
260  variable nom, denom
261  nom = cos(pi * sing / 2 + (1 - sing) * atan((x - pos) / fwhm * 2))
262  denom = ((x - pos)^2 + fwhm^2 / 4)^((1 - sing) / 2)
263 
264  return amp * nom / denom * fwhm / 2
265 end
266 
280 threadsafe function MultiDoniachSunjicLinBG(w,x) : FitFunc
281  wave w
282  variable x
283 
284  variable np = numpnts(w)
285  variable ip
286 
287  variable v = w[0] + x * w[1]
288  for (ip = 2; ip < np; ip += 4)
289  v += DoniachSunjic(x, w[ip], w[ip+1], w[ip+3], w[ip+2])
290  endfor
291 
292  return v
293 end
294 
295 
296 threadsafe function ds1_bg(w, x): FitFunc
297  // Doniach-Sunjic fit function
298  // 0 <= sing < 1
299  wave w // coefficients - see below
300  variable x // independent variable
301 
302  //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
303  //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
304  //CurveFitDialog/ Equation:
305  //CurveFitDialog/ f(x) = DoniachSunjic(x, amp, pos, sing, fwhm) + bg
306  //CurveFitDialog/ End of Equation
307  //CurveFitDialog/ Independent Variables 1
308  //CurveFitDialog/ x
309  //CurveFitDialog/ Coefficients 5
310  //CurveFitDialog/ w[0] = bg
311  //CurveFitDialog/ w[1] = amp
312  //CurveFitDialog/ w[2] = pos
313  //CurveFitDialog/ w[3] = sing
314  //CurveFitDialog/ w[4] = FWHM
315 
316  return DoniachSunjic(x, w[1], w[2], w[3], w[4]) + w[0]
317 end
318 
319 threadsafe function ds2_bg(w,x) : FitFunc
320  Wave w
321  Variable x
322 
323  //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
324  //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
325  //CurveFitDialog/ Equation:
326  //CurveFitDialog/ f(x) = w_0+( w_1*cos(pi*w_3/2+(1-w_3)*atan((x-w_2)/w_4)))/(((x-w_2)^2+w_4^2)^((1-w_3)/2)) +(w_5*cos(pi*w_7/2+(1-w_7)*atan((x-(w_6))/w_8)))/(((x-w_6)^2+w_8^2)^((1-w_7)/2))
327  //CurveFitDialog/ End of Equation
328  //CurveFitDialog/ Independent Variables 1
329  //CurveFitDialog/ x
330  //CurveFitDialog/ Coefficients 9
331  //CurveFitDialog/ w[0] = bg
332  //CurveFitDialog/ w[1] = amp1
333  //CurveFitDialog/ w[2] = pos1
334  //CurveFitDialog/ w[3] = sing1
335  //CurveFitDialog/ w[4] = wid1
336  //CurveFitDialog/ w[5] = amp2
337  //CurveFitDialog/ w[6] = pos2
338  //CurveFitDialog/ w[7] = sing2
339  //CurveFitDialog/ w[8] = wid2
340 
341  variable ds1 = DoniachSunjic(x, w[1], w[2], w[3], w[4])
342  variable ds2 = DoniachSunjic(x, w[5], w[6], w[7], w[8])
343 
344  return w[0] + ds1 + ds2
345 End
346 
347 Function ds4_bg(w,x) : FitFunc
348  Wave w
349  Variable x
350 
351  //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
352  //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
353  //CurveFitDialog/ Equation:
354  //CurveFitDialog/ f(x) = w_0+( w_1*cos(pi*w_3/2+(1-w_3)*atan((x-w_2)/w_4)))/(((x-w_2)^2+w_4^2)^((1-w_3)/2)) +(w_5*cos(pi*w_7/2+(1-w_7)*atan((x-(w_6))/w_8)))/(((x-w_6)^2+w_8^2)^((1-w_7)/2)) +( w_9*cos(pi*w_11/2+(1-w_11)*atan((x-w_10)/w_12)))/(((x-w_10)^2+w_12^2)^((1-w_11)/2)) +( w_13*cos(pi*w_15/2+(1-w_15)*atan((x-w_14)/w_16)))/(((x-w_14)^2+w_16^2)^((1-w_15)/2))
355  //CurveFitDialog/ End of Equation
356  //CurveFitDialog/ Independent Variables 1
357  //CurveFitDialog/ x
358  //CurveFitDialog/ Coefficients 17
359  //CurveFitDialog/ w[0] = w_0
360  //CurveFitDialog/ w[1] = w_11
361  //CurveFitDialog/ w[2] = w_12
362  //CurveFitDialog/ w[3] = w_13
363  //CurveFitDialog/ w[4] = w_14
364  //CurveFitDialog/ w[5] = w_21
365  //CurveFitDialog/ w[6] = w_22
366  //CurveFitDialog/ w[7] = w_23
367  //CurveFitDialog/ w[8] = w_24
368  //CurveFitDialog/ w[9] = w_31
369  //CurveFitDialog/ w[10] = w_32
370  //CurveFitDialog/ w[11] = w_33
371  //CurveFitDialog/ w[12] = w_34
372  //CurveFitDialog/ w[13] = w_41
373  //CurveFitDialog/ w[14] = w_42
374  //CurveFitDialog/ w[15] = w_43
375  //CurveFitDialog/ w[16] = w_44
376  Variable ds1, ds2, ds3, ds4
377  ds1=( w[1]*cos(pi*w[3]/2+(1-w[3])*atan((x-w[2])/w[4])))/(((x-w[2])^2+w[4]^2)^((1-w[3])/2))
378  ds2=( w[5]*cos(pi*w[7]/2+(1-w[7])*atan((x-w[6])/w[8])))/(((x-w[6])^2+w[8]^2)^((1-w[7])/2))
379  ds3=( w[9]*cos(pi*w[11]/2+(1-w[11])*atan((x-w[10])/w[12])))/(((x-w[10])^2+w[12]^2)^((1-w[11])/2))
380  ds4=( w[13]*cos(pi*w[15]/2+(1-w[15])*atan((x-w[14])/w[16])))/(((x-w[14])^2+w[16]^2)^((1-w[15])/2))
381 
382 
383  return w[0]+ds1+ds2+ds3+ds4
384 
385 
386 End
387 
388 Function ds6_bg(w,x) : FitFunc
389  Wave w
390  Variable x
391 
392  //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
393  //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
394  //CurveFitDialog/ Equation:
395  //CurveFitDialog/
396  //CurveFitDialog/ Variable g, ds1, ds2, ds3, ds4, ds5, ds6
397  //CurveFitDialog/ ds1=( w_11*cos(pi*w_13/2+(1-w_13)*atan((x-w_12)/w_14)))/(((x-w_12)^2+w_14^2)^((1-w_13)/2))
398  //CurveFitDialog/ ds2=( w_21*cos(pi*w_23/2+(1-w_23)*atan((x-w_22)/w_24)))/(((x-w_22)^2+w_24^2)^((1-w_23)/2))
399  //CurveFitDialog/ ds3=( w_31*cos(pi*w_33/2+(1-w_33)*atan((x-w_32)/w_34)))/(((x-w_32)^2+w_34^2)^((1-w_33)/2))
400  //CurveFitDialog/ ds4=( w_41*cos(pi*w_43/2+(1-w_43)*atan((x-w_42)/w_44)))/(((x-w_42)^2+w_44^2)^((1-w_43)/2))
401  //CurveFitDialog/ ds5=( w_51*cos(pi*w_53/2+(1-w_53)*atan((x-w_52)/w_54)))/(((x-w_52)^2+w_54^2)^((1-w_53)/2))
402  //CurveFitDialog/ ds6=( w_61*cos(pi*w_63/2+(1-w_63)*atan((x-w_62)/w_64)))/(((x-w_62)^2+w_64^2)^((1-w_63)/2))
403  //CurveFitDialog/
404  //CurveFitDialog/ f(x) =w_0+ds1+ds2+ds3+ds4+ds5+ds6
405  //CurveFitDialog/
406  //CurveFitDialog/ End of Equation
407  //CurveFitDialog/ Independent Variables 1
408  //CurveFitDialog/ x
409  //CurveFitDialog/ Coefficients 25
410  //CurveFitDialog/ w[0] = w_0
411  //CurveFitDialog/ w[1] = w_11
412  //CurveFitDialog/ w[2] = w_12
413  //CurveFitDialog/ w[3] = w_13
414  //CurveFitDialog/ w[4] = w_14
415  //CurveFitDialog/ w[5] = w_21
416  //CurveFitDialog/ w[6] = w_22
417  //CurveFitDialog/ w[7] = w_23
418  //CurveFitDialog/ w[8] = w_24
419  //CurveFitDialog/ w[9] = w_31
420  //CurveFitDialog/ w[10] = w_32
421  //CurveFitDialog/ w[11] = w_33
422  //CurveFitDialog/ w[12] = w_34
423  //CurveFitDialog/ w[13] = w_41
424  //CurveFitDialog/ w[14] = w_42
425  //CurveFitDialog/ w[15] = w_43
426  //CurveFitDialog/ w[16] = w_44
427  //CurveFitDialog/ w[17] = w_51
428  //CurveFitDialog/ w[18] = w_52
429  //CurveFitDialog/ w[19] = w_53
430  //CurveFitDialog/ w[20] = w_54
431  //CurveFitDialog/ w[21] = w_61
432  //CurveFitDialog/ w[22] = w_62
433  //CurveFitDialog/ w[23] = w_63
434  //CurveFitDialog/ w[24] = w_64
435 
436 
437  Variable ds1, ds2, ds3, ds4, ds5, ds6
438  ds1=( w[1]*cos(pi*w[3]/2+(1-w[3])*atan((x-w[2])/w[4])))/(((x-w[2])^2+w[4]^2)^((1-w[3])/2))
439  ds2=( w[5]*cos(pi*w[7]/2+(1-w[7])*atan((x-w[6])/w[8])))/(((x-w[6])^2+w[8]^2)^((1-w[7])/2))
440  ds3=( w[9]*cos(pi*w[11]/2+(1-w[11])*atan((x-w[10])/w[12])))/(((x-w[10])^2+w[12]^2)^((1-w[11])/2))
441  ds4=( w[13]*cos(pi*w[15]/2+(1-w[15])*atan((x-w[14])/w[16])))/(((x-w[14])^2+w[16]^2)^((1-w[15])/2))
442  ds5=( w[17]*cos(pi*w[19]/2+(1-w[19])*atan((x-w[18])/w[20])))/(((x-w[18])^2+w[20]^2)^((1-w[19])/2))
443  ds6=( w[21]*cos(pi*w[23]/2+(1-w[23])*atan((x-w[22])/w[24])))/(((x-w[22])^2+w[24]^2)^((1-w[23])/2))
444 
445  return w[0]+ds1+ds2+ds3+ds4+ds5+ds6
446 
447 End
448 
450  // data structure for DoniachSunjicBroadS structural function fit
451 
452  // waves populated by the FuncFit operation
453  wave pw
454  wave yw
455  wave xw
456 
457  // convolution parameters to be set upon creation of the structure
458  variable precision
459  variable oversampling
460 
461  // auxiliary fields used internally by DoniachSunjicBroadS
462  // do not touch these
463  wave xdw
464  wave model
465  wave broadening
466  wave convolution
467 EndStructure
468 
469 //------------------------------------------------------------------------------
470 threadsafe function DoniachSunjicBroadS(s) : FitFunc
471 //------------------------------------------------------------------------------
472  // Two-peak (bulk + surface) Doniach-Sunjic line shape with Gaussian broadening (convolution).
473  // Hold parameter 5 at 0 to fit just one peak.
474 
475  // Structural fit function for efficient fitting in procedures.
476  // Calculating the convolution requires auxiliary waves and additional, non-fitting parameters.
477  // To eliminate the time-consuming overhead of creating and killing the auxiliary waves,
478  // these waves are held in the fitting structure.
479 
480  // Caution: The function on its own is thread-safe.
481  // However, since FuncFit uses the same structure in all threads, the fitting cannot run in parallel.
482  // Set /NTHR=1.
483 
484  // See also Fit_DoniachSunjicBroad (example), DoniachSunjicBroad (conventional fit function)
485  Struct DoniachSunjicStruct &s
486 
487  // pw[0] = bulk amplitude
488  // pw[1] = bulk position
489  // pw[2] = Lorentzian FWHM
490  // pw[3] = Donjach-Sunjic singularity index (0..1)
491  // pw[4] = surface shift
492  // pw[5] = surface/bulk ratio
493  // pw[6] = Gaussian FWHM
494  // pw[7] = constant background
495  // pw[8] = linear background
496 
497  wave xw = s.xw
498  wave yw = s.yw
499  wave pw = s.pw
500 
501  variable precision = s.precision
502  variable oversampling = s.oversampling
503 
504  if (WaveExists(s.xdw))
505  wave xdw = s.xdw
506  wave model = s.model
507  wave broadening = s.broadening
508  wave convolution = s.convolution
509  else
510  make /n=0 /free xdw, model, broadening, convolution
511  redimension /d xdw, model, broadening, convolution
512  wave fs.xdw = xdw
513  wave fs.model = model
514  wave fs.broadening = broadening
515  wave fs.convolution = convolution
516  endif
517 
518  // calculate wave spacing based on minimum spacing of desired x points
519  differentiate /p xw /d=xdw
520  xdw = abs(xdw)
521  variable xd = wavemin(xdw) / oversampling
522 
523  // calculate broadening wave size based on width and precision variable
524  variable x0b = pw[6] * precision
525  variable nb = 2 * floor(x0b / xd) + 1
526 
527  // calculate model size based on desired range for yw
528  variable x0m = max(abs(wavemax(xw) - pw[1]), abs(wavemin(xw) - pw[1])) + x0b
529  variable nm = 2 * floor(x0m / xd) + 1
530  nb = min(nb, nm * 10) // limit wave size to avoid runtime errors for unphysically large parameter
531 
532  // create and calculate initial waves, normalize exponential
533  redimension /n=(nb) broadening
534  redimension /n=(nm) model
535  setscale/i x -x0b, x0b, "", broadening
536  setscale/i x -x0m, x0m, "", model
537 
538  broadening = exp( - (x / pw[6])^2 * 4 * ln(2))
539  variable nrm = area(broadening)
540  broadening /= nrm
541  model = DoniachSunjic(x, 1, 0, pw[3], pw[2]) // bulk
542  model += DoniachSunjic(x, pw[5], pw[4], pw[3], pw[2]) // surface
543 
544  // calculate the convolution
545  Convolve /a broadening, model
546  variable scale = pw[0] / wavemax(model)
547  model *= scale
548 
549  // prepare output
550  nm = numpnts(model)
551  x0m = xd * (nm - 1) / 2
552  setscale/i x -x0m, x0m, "", model
553 
554  yw = model(xw[p] - pw[1]) + pw[7] + pw[8] * xw[p]
555  yw = numtype(yw) ? 0 : yw
556 end
557 
558 //------------------------------------------------------------------------------
559 function DoniachSunjicBroad(pw, yw, xw) : FitFunc
560 //------------------------------------------------------------------------------
561  // Two-peak (bulk + surface) Doniach-Sunjic line shape with Gaussian broadening (convolution).
562  // Hold parameter 5 at 0 to fit just one peak.
563  // Conventional fit function for use with the curve-fitting dialog.
564  // Compared to DoniachSunjicBroadS this function incurs extra overhead
565  // because auxiliary waves are created and killed between function calls.
566  // See also DoniachSunjicBroadS (optimized structural fit function)
567  Wave pw
568  Wave yw
569  Wave xw
570 
571  // pw[0] = bulk amplitude
572  // pw[1] = bulk position
573  // pw[2] = Lorentzian FWHM
574  // pw[3] = Donjach-Sunjic singularity index (0..1)
575  // pw[4] = surface shift
576  // pw[5] = surface/bulk ratio
577  // pw[6] = Gaussian FWHM
578  // pw[7] = constant background
579  // pw[8] = linear background
580 
581  // set up data structure
582  struct DoniachSunjicStruct fs
583  fs.precision = 5
584  fs.oversampling = 4
585 
586  wave fs.pw = pw
587  wave fs.xw = xw
588  wave fs.yw = yw
589 
590  // create temporary calculation waves in a global folder
591  dfref df = root:packages:pearl_fitfuncs:doniach_sunjic
592  if (DataFolderRefStatus(df) == 0)
593  newdatafolder root:packages:pearl_fitfuncs:doniach_sunjic
594  dfref df = root:packages:pearl_fitfuncs:doniach_sunjic
595  endif
596 
597  wave /z /sdfr=df fs.xdw = xdw
598  wave /z /sdfr=df fs.model = model
599  wave /z /sdfr=df fs.broadening = broadening
600  wave /z /sdfr=df fs.convolution = convolution
601 
602  if (WaveExists(fs.xdw) == 0)
603  dfref savedf = GetDataFolderDFR()
604  setdatafolder df
605  make /n=0 /d xdw, model, broadening, convolution
606  wave fs.xdw = xdw
607  wave fs.model = model
608  wave fs.broadening = broadening
609  wave fs.convolution = convolution
610  setdatafolder savedf
611  endif
612 
613  // calculate
615 
616  yw = fs.yw
617 end
618 
619 //------------------------------------------------------------------------------
620 function Calc_DoniachSunjicBroad(pw, yw)
621 //------------------------------------------------------------------------------
622  // Calculate the DoniachSunjicBroadS line shape
623  Wave pw // coefficient wave
624  Wave yw // output wave, correct x-scaling required on input
625 
626  struct DoniachSunjicStruct fs
627  fs.precision = 5
628  fs.oversampling = 4
629 
630  duplicate /free pw, fs.pw
631  duplicate /free yw, fs.xw
632  fs.xw = x
633  duplicate /free yw, fs.yw
634 
636 
637  yw = fs.yw
638 end
639 
640 //------------------------------------------------------------------------------
641 Function Fit_DoniachSunjicBroad(pw, yw, xw, ww)
642 //------------------------------------------------------------------------------
643  // Fit the DoniachSunjicBroadS line shape.
644  // The function applies constraints which assume that the energy scale is in eV.
645  // Returns chi^2.
646  wave pw // coefficient wave- pre-load it with initial guess
647  wave yw
648  wave /z xw
649  wave /z ww // weights (standard deviation)
650 
651  struct DoniachSunjicStruct fs
652  fs.precision = 5
653  fs.oversampling = 4
654 
655  duplicate /free pw, fs.pw
656  if (WaveExists(xw))
657  duplicate /free xw, fs.xw
658  else
659  duplicate /free yw, fs.xw
660  fs.xw = x
661  endif
662  duplicate /free yw, fs.yw
663 
664  variable v_chisq = nan
665  variable V_FitMaxIters = 100
666  make /n=1 /t /free constraints = {"K0 >= 0", "K2 > 0", "K2 < 10", "K3 >= 0", "K3 < 1", "K4 >= -10", "K4 <= 10", "K5 >= 0", "K5 <= 1", "K6 >= 0", "K6 < 10"}
667  // note: only single thread allowed
668  FuncFit /NTHR=1 DoniachSunjicBroadS, pw, yw /X=xw /D /STRC=fs /C=constraints /NWOK /I=1 /W=ww
669 
670  return v_chisq
671 End
672 
673 //------------------------------------------------------------------------------
674 // peak-specific fit functions
675 //------------------------------------------------------------------------------
676 
677 function Au4f(w, x): fitfunc
678  // fit function for a nitrogen 1s-pi* absorption spectrum
679  // modelled as multiple Voigt shapes on a constant background
680  // similar to the Igor VoigtFit function
681  // but with a constant gaussian width (instrumental broadening) for all peaks
682  // gaussian and lorentzian widths are specified as FWHM
683  wave w // parameters
684  // w[0] constant background
685  // w[1] linear background
686  // w[2] global gaussian FWHM
687  // w[3 + 0 + (n-1) * 3] peak n area
688  // w[3 + 1 + (n-1) * 3] peak n position
689  // w[3 + 2 + (n-1) * 3] peak n lorentzian FWHM
690  // length of wave defines number of peaks
691 
692  // for compatibility with older code the linear background term can be omitted.
693  // if the number of parameters divides by 3, the linear background term is added,
694  // otherwise only the constant background.
695  variable x
696 
697  variable np = numpnts(w)
698  variable ip, ip0
699 
700  variable bg = w[0]
701  variable v = bg
702  if (mod(np, 3) == 0)
703  v += w[1] * x
704  ip0 = 3
705  else
706  ip0 = 2
707  endif
708 
709  variable vc1, vc2, vc3, vc4
710  vc2 = 2 * sqrt(ln(2)) / w[ip0-1]
711  for (ip = ip0; ip < np; ip += 3)
712  vc1 = w[ip] / sqrt(pi) * vc2
713  vc3 = w[ip+1]
714  vc4 = vc2 * w[ip+2] / 2
715  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
716  endfor
717 
718  return v
719 
720 end
721 
722 function Au4f_2p2(w, x): fitfunc
723  // Au 4f 5/2 and 7/2 2-component Voigt fit with a common gaussian width
724  // gaussian and lorentzian widths are specified as FWHM
725  wave w // parameters
726  // w[0] constant background
727  // w[1] linear background
728  // w[2] global gaussian FWHM
729  // w[3] 5/2 bulk area
730  // w[4] 5/2 bulk position
731  // w[5] 5/2 lorentzian FWHM
732  // w[6] 7/2 bulk area
733  // w[7] 7/2 bulk position
734  // w[8] 7/2 lorentzian FWHM
735  // w[9] surface/bulk area ratio
736  // w[10] surface core level shift
737  variable x
738 
739  variable bg = w[0] + w[1] * x
740  variable v = bg
741 
742  variable vc1 // amplitude
743  variable vc2 // width
744  variable vc3 // position
745  variable vc4 // shape
746  vc2 = 2 * sqrt(ln(2)) / w[2]
747 
748  // 5/2 bulk
749  vc1 = w[3] / sqrt(pi) * vc2
750  vc3 = w[4]
751  vc4 = vc2 * w[5] / 2
752  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
753 
754  // 5/2 surface
755  vc1 = w[3] / sqrt(pi) * vc2 * w[9]
756  vc3 = w[4] + w[10]
757  vc4 = vc2 * w[5] / 2
758  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
759 
760  // 7/2 bulk
761  vc1 = w[6] / sqrt(pi) * vc2
762  vc3 = w[7]
763  vc4 = vc2 * w[8] / 2
764  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
765 
766  // 7/2 surface
767  vc1 = w[6] / sqrt(pi) * vc2 * w[9]
768  vc3 = w[7] + w[10]
769  vc4 = vc2 * w[8] / 2
770  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
771 
772  return v
773 
774 end
775 
776 function ShowComponents_Au4f_2p2(coef_wave, fit_wave)
777  wave coef_wave
778  wave fit_wave
779 
780  duplicate /free coef_wave, coef1, coef2
781  coef1[9] = 0
782  coef2[3] *= coef_wave[9]
783  coef2[4] += coef_wave[10]
784  coef2[6] *= coef_wave[9]
785  coef2[7] += coef_wave[10]
786  coef2[9] = 0
787 
788  string s_fit_wave = NameOfWave(fit_wave)
789  string s_fit_p1 = s_fit_wave + "_p1"
790  string s_fit_p2 = s_fit_wave + "_p2"
791  duplicate /o fit_wave, $(s_fit_p1) /wave=fit_p1
792  duplicate /o fit_wave, $(s_fit_p2) /wave=fit_p2
793 
794  fit_p1 = Au4f_2p2(coef1, x)
795  fit_p2 = Au4f_2p2(coef2, x)
796 
797  string traces = TraceNameList("", ";", 1)
798  if ((WhichListItem(s_fit_wave, traces, ";") >= 0) && (WhichListItem(s_fit_p1, traces, ";") < 0))
799  appendtograph fit_p1, fit_p2
800  ModifyGraph lstyle($s_fit_p1)=2
801  ModifyGraph lstyle($s_fit_p2)=2
802  ModifyGraph rgb($s_fit_p1)=(0,0,65280)
803  ModifyGraph rgb($s_fit_p2)=(0,0,65280)
804  endif
805 end
806 
807 function Au4f_2p3(w, x): fitfunc
808  // Au 4f 5/2 and 7/2 3-component Voigt fit with a common gaussian width
809  // gaussian and lorentzian widths are specified as FWHM
810  wave w // parameters
811  // w[0] constant background
812  // w[1] linear background
813  // w[2] global gaussian FWHM
814  // w[3] 5/2 bulk area
815  // w[4] 5/2 bulk position
816  // w[5] 5/2 lorentzian FWHM
817  // w[6] 7/2 bulk area
818  // w[7] 7/2 bulk position
819  // w[8] 7/2 lorentzian FWHM
820  // w[9] surface/bulk area ratio
821  // w[10] surface core level shift
822  // w[11] 2nd layer/bulk area ratio
823  // w[12] 2nd layer core level shift
824  variable x
825 
826  variable bg = w[0] + w[1] * x
827  variable v = bg
828 
829  variable vc1 // amplitude
830  variable vc2 // width
831  variable vc3 // position
832  variable vc4 // shape
833  vc2 = 2 * sqrt(ln(2)) / w[2]
834 
835  // 5/2 bulk
836  vc1 = w[3] / sqrt(pi) * vc2
837  vc3 = w[4]
838  vc4 = vc2 * w[5] / 2
839  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
840 
841  // 5/2 surface
842  vc1 = w[3] / sqrt(pi) * vc2 * w[9]
843  vc3 = w[4] + w[10]
844  vc4 = vc2 * w[5] / 2
845  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
846 
847  // 5/2 2nd layer
848  vc1 = w[3] / sqrt(pi) * vc2 * w[11]
849  vc3 = w[4] + w[12]
850  vc4 = vc2 * w[5] / 2
851  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
852 
853  // 7/2 bulk
854  vc1 = w[6] / sqrt(pi) * vc2
855  vc3 = w[7]
856  vc4 = vc2 * w[8] / 2
857  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
858 
859  // 7/2 surface
860  vc1 = w[6] / sqrt(pi) * vc2 * w[9]
861  vc3 = w[7] + w[10]
862  vc4 = vc2 * w[8] / 2
863  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
864 
865  // 7/2 2nd layer
866  vc1 = w[6] / sqrt(pi) * vc2 * w[11]
867  vc3 = w[7] + w[12]
868  vc4 = vc2 * w[8] / 2
869  v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
870 
871  return v
872 
873 end
874 
875 function ShowComponents_Au4f_2p3(coef_wave, fit_wave)
876  wave coef_wave
877  wave fit_wave
878 
879  duplicate /free coef_wave, coef1, coef2, coef3
880  coef1[9] = 0
881  coef1[11] = 0
882 
883  coef2[3] *= coef_wave[9]
884  coef2[4] += coef_wave[10]
885  coef2[6] *= coef_wave[9]
886  coef2[7] += coef_wave[10]
887  coef2[9] = 0
888  coef2[11] = 0
889 
890  coef3[3] *= coef_wave[11]
891  coef3[4] += coef_wave[12]
892  coef3[6] *= coef_wave[11]
893  coef3[7] += coef_wave[12]
894  coef3[9] = 0
895  coef3[11] = 0
896 
897  string s_fit_wave = NameOfWave(fit_wave)
898  string s_fit_p1 = s_fit_wave + "_p1"
899  string s_fit_p2 = s_fit_wave + "_p2"
900  string s_fit_p3 = s_fit_wave + "_p3"
901  duplicate /o fit_wave, $(s_fit_p1) /wave=fit_p1
902  duplicate /o fit_wave, $(s_fit_p2) /wave=fit_p2
903  duplicate /o fit_wave, $(s_fit_p3) /wave=fit_p3
904 
905  fit_p1 = Au4f_2p2(coef1, x)
906  fit_p2 = Au4f_2p2(coef2, x)
907  fit_p3 = Au4f_2p2(coef3, x)
908 
909  string traces = TraceNameList("", ";", 1)
910  if ((WhichListItem(s_fit_wave, traces, ";") >= 0) && (WhichListItem(s_fit_p1, traces, ";") < 0))
911  appendtograph fit_p1, fit_p2, fit_p3
912  ModifyGraph lstyle($s_fit_p1)=2
913  ModifyGraph lstyle($s_fit_p2)=2
914  ModifyGraph lstyle($s_fit_p3)=2
915  ModifyGraph rgb($s_fit_p1)=(0,0,65280)
916  ModifyGraph rgb($s_fit_p2)=(0,0,65280)
917  ModifyGraph rgb($s_fit_p3)=(0,0,65280)
918  endif
919 end
920 
930 function FermiGaussConv(pw, yw, xw) : FitFunc
931  WAVE pw, yw, xw
932 
933  // half width of temporary gaussian wave is pw[5] multiplied by this factor (may be fractional)
934  variable precision_g = 5
935  variable oversampling = 4
936 
937  // calculate wave spacing based on minimum spacing of desired x points
938  duplicate /free xw, xdw
939  differentiate /p xw /d=xdw
940  xdw = abs(xdw)
941  variable xd = wavemin(xdw) / oversampling
942 
943  // calculate gausswave size based on pw[5] and precision variable
944  variable x0g = abs(pw[5]) * precision_g
945  variable ng = 2 * floor(x0g / xd) + 1
946 
947  // calculate fermiwave size based on desired range for yw
948  variable emax = wavemax(xw)
949  variable emin = wavemin(xw)
950  variable x0f = max(abs(emax - pw[3]), abs(emin - pw[3])) + x0g
951  variable ne = 2 * floor(x0f / xd) + 1
952 
953  // create and calculate initial waves, normalize exponential
954  make /d /n=(ng) /free gausswave
955  make /d /n=(ne) /free fermiwave
956  setscale/i x -x0g, x0g, "", gausswave
957  setscale/i x -x0f, x0f, "", fermiwave
958 
959  gausswave = exp( - (x / pw[5] )^2 )
960  fermiwave = 1 / (exp( x / (kBoltzmann * pw[4])) + 1.0 )
961 
962  // calculate the convolution
963  duplicate /free fermiwave, resultwave
964  Convolve /a gausswave, resultwave
965  variable rmax = wavemax(resultwave)
966  resultwave /= rmax
967 
968  // prepare output
969  ng = numpnts(resultwave)
970  x0g = xd * (ng - 1) / 2
971  setscale/i x -x0g, x0g, "", resultwave
972 
973  yw = pw[2] * resultwave(xw[p] - pw[3]) + pw[0] + pw[1] * xw[p]
974 end
975 
976 
980 function ShirleyBG(w, bg, p1, p2)
981  wave w
982  wave bg
983  variable p1, p2
984 
985  duplicate /o w, bg
986  integrate /meth=1 w /d=bg
987 
988  variable bg1 = bg[p1]
989  variable bg2 = bg[p2]
990  bg -= bg1
991  bg /= bg2 - bg1
992  bg *= w[p2] - w[p1]
993  bg += w[p1]
994 end
DoniachSunjicBroadS
threadsafe variable DoniachSunjicBroadS(DoniachSunjicStruct *s)
Definition: pearl-fitfuncs.ipf:470
DoniachSunjicBroad
variable DoniachSunjicBroad(wave pw, wave yw, wave xw)
Definition: pearl-fitfuncs.ipf:559
DblDoubletGaussLinBG_AO
threadsafe variable DblDoubletGaussLinBG_AO(wave pw, wave yw, wave xw)
two doublet gaussian peaks on a linear background fit function (all at once).
Definition: pearl-fitfuncs.ipf:157
Calc_DoniachSunjicBroad
variable Calc_DoniachSunjicBroad(wave pw, wave yw)
Definition: pearl-fitfuncs.ipf:620
MultiVoigtLinBG_AO
threadsafe variable MultiVoigtLinBG_AO(wave pw, wave yw, wave xw)
multiple voigt peaks on a linear background fit function.
Definition: pearl-fitfuncs.ipf:224
MultiGaussLinBG_AO
threadsafe variable MultiGaussLinBG_AO(wave pw, wave yw, wave xw)
multiple gaussian peaks on a linear background fit function (all at once).
Definition: pearl-fitfuncs.ipf:80
DoniachSunjicStruct
Definition: pearl-fitfuncs.ipf:449
DoubletGaussLinBG_AO
threadsafe variable DoubletGaussLinBG_AO(wave pw, wave yw, wave xw)
doublet gaussian peaks on a linear background fit function (all at once).
Definition: pearl-fitfuncs.ipf:118
MultiVoigtLinBG
threadsafe variable MultiVoigtLinBG(wave w, variable x)
multiple voigt peaks on a linear background fit function.
Definition: pearl-fitfuncs.ipf:190
ds2_bg
threadsafe variable ds2_bg(wave w, variable x)
Definition: pearl-fitfuncs.ipf:319
ShowComponents_Au4f_2p3
variable ShowComponents_Au4f_2p3(wave coef_wave, wave fit_wave)
Definition: pearl-fitfuncs.ipf:875
Fit_DoniachSunjicBroad
variable Fit_DoniachSunjicBroad(wave pw, wave yw, wave xw, wave ww)
Definition: pearl-fitfuncs.ipf:641
Au4f_2p3
variable Au4f_2p3(wave w, variable x)
Definition: pearl-fitfuncs.ipf:807
ds6_bg
variable ds6_bg(wave w, variable x)
Definition: pearl-fitfuncs.ipf:388
Au4f
variable Au4f(wave w, variable x)
Definition: pearl-fitfuncs.ipf:677
Au4f_2p2
variable Au4f_2p2(wave w, variable x)
Definition: pearl-fitfuncs.ipf:722
ShowComponents_Au4f_2p2
variable ShowComponents_Au4f_2p2(wave coef_wave, wave fit_wave)
Definition: pearl-fitfuncs.ipf:776
ds1_bg
threadsafe variable ds1_bg(wave w, variable x)
Definition: pearl-fitfuncs.ipf:296
ShirleyBG
variable ShirleyBG(wave w, wave bg, variable p1, variable p2)
calculate the shirley background
Definition: pearl-fitfuncs.ipf:980
ds4_bg
variable ds4_bg(wave w, variable x)
Definition: pearl-fitfuncs.ipf:347
MultiDoniachSunjicLinBG
threadsafe variable MultiDoniachSunjicLinBG(wave w, variable x)
multiple doniach-sunjic peaks on a linear background fit function.
Definition: pearl-fitfuncs.ipf:280
DoniachSunjic
threadsafe variable DoniachSunjic(variable x, variable amp, variable pos, variable sing, variable fwhm)
Doniach-Sunjic line shape.
Definition: pearl-fitfuncs.ipf:253
MultiGaussLinBG
threadsafe variable MultiGaussLinBG(wave w, variable x)
multiple gaussian peaks on a linear background fit function.
Definition: pearl-fitfuncs.ipf:45
FermiGaussConv
variable FermiGaussConv(wave pw, wave yw, wave xw)
convolution of Fermi-Dirac distribution and a Gaussian.
Definition: pearl-fitfuncs.ipf:930