| Author |
 Topic  |
|
|
gingko
Posts |
 Posted - 07/01/2006 : 03:49:20 AM
|
Origin Version (Select Help-->About Origin): 7.0
Operating System:xp
Hello to all,
I have a problem with user-defined fitting functions. What I would like to do is to use the square of Voigt function to fit the peak profile. But I have a problem writing the integration function from negative infinity to infinity. How do I edit the function in Code Builder?
The equation of Voigt function is as follows:
Thanks in advance
Gingko
Edited by - gingko on 07/01/2006 03:50:18 AM [code]
Edited by - gingko on 07/04/2006 02:11:27 AM
Edited by - gingko on 07/04/2006 02:12:01 AM |
|
|
obauer
Germany
15 Posts |
Posted - 05/22/2009 : 11:01:08 AM
|
Origin Version: OriginPro 8 SR4
Operating System: Windows XP
Hello!
I am also trying to modify the Voigt function (VOIGT5.FDF) for fitting since I need to fit experimental data which is a convolution of two functions (one of which is a Gaussian). As the Voigt function is a convolution of Gaussian and Lorentzian by definition I had the idea to start from this built-in function. It makes use of the followinf formula:
y = y0 + (A*2*ln(2)*wL)/(Pi^1.5*wG^2) * integ( exp(-t^2) / ((sqrt(ln(2))*wL/wG)^2+(sqrt(4*ln(2))*(x-xc)/wG-t)^2) )
Yet as soon as I duplicate the Voigt function in the NLSF assistant and try to use it - without any further modification - I get the following error message:
"Error, function or variable integ not found"
What is the reason for this? How comes the Voigt function is computed but not its one-to-one copy? Is there any other chance to create a user-defined fitting function which includes a convolution and which then can be called from an Origin C routine? I have read the posts about using the trapezoidal rule but can´t I just use the "integ" command like in VOIGT5.FDF ... ? The computation of the trapezoidal rule is very time-consuming....
Has anyone an idea how to sort this out? I would be more than glad to hear your opinion and experience on this. Thanks in advance!
Cheers, Oliver |
 |
|
|
easwar
USA
1964 Posts |
|
|
obauer
Germany
15 Posts |
Posted - 05/27/2009 : 07:28:51 AM
|
Hello Easwar!
Thank you for your quick reply. Unfortunately I don´t have an analytic expression for the function that takes the convolution into account. This is why I don't see a chance to use the example which you have suggested...
My problem is as follows:
I need to fit an experimental signal which is a convolution of a theorectical signal F convoluted with the instrumental broadening R. The latter is assumed to be of Gaussian type. Yet the theoretical signal is calculated with the use of complex numbers although the exp. data points are real numbers of course. I have no idea how to compute this with an analytical formula in Origin C and NLSF respectively. This is why I evenetually decided to use the trapezoidal rule and I came up with the following code in the NLSF:
******* start of Origin C code *****
[Formula]
double dIntegral = 0.0;
double dInt = 0.0;
double t = 3.032; // the integration interval ranges from 3.032 to 3.038 keV (= exp. data range)
double dPrecision = 1e-6; // the precision to control the integraton to infinity
double dt = 0.0002; // this is the step size for the integration in keV.
//the following lines actually perform a integrate to the function, the trapezoidal rule is used.
double F1, F2, R1, R2; // F = function (diffraction theory), R = response (Gaussian)
double tcL = 3.03441;
double gamma = 0.00056;
do
{
R1 = 1 / (wG * sqrt(2.0 * PI)) * exp( -0.5*((x - xcG) - t)^2/wG^2 );
R2 = 1 / (wG * sqrt(2.0 * PI)) * exp( -0.5*((x - xcG) - (t + dt))^2/wG^2 );
//*************************
double E_Bragg = 3.03489;
double theta_Bragg = 90 * PI/180;
double Sigma = 2.195624e-006;
double F0_real = 151.4992;
double F0_imag = 15.33788;
double FH_real = 105.4824;
double FH_imag = 15.33788;
double DW = 0.98;
double b = -1.0;
double Delta_E1,Delta_E2; // deviation from calculated Bragg energy in eV
double E_photon1,E_photon2; // photon energy in eV = E_Bragg + Delta_E
// definition of complex structure factos F0, FH
complex F0(F0_real,F0_imag);
//out_complex("F0 = ", F0);
complex FH(FH_real,FH_imag);
//out_complex("FH = ", FH);
// temperatue correction of reflected intensity
FH = FH*DW;
//out_complex("FH*DW = ", FH);
// Calculation of eta and Reflectivity = |b| * |eta +- (eta^2 -1)^0.5|^2
complex eta1, eta2;
complex z01,z02; // z0 = ((-2.0*pow(sin(theta_Bragg),2.0)*Delta_E[ii]/E_photon[ii] + Sigma * F0)
complex z11,z12; // z1 = eta^2
complex z21,z22; // z2 = eta^2 - 1
complex z31,z32; // z3 = (eta^2 - 1)^0.5
double d11,d12; // d1 = | eta +- (eta^2 - 1)^0.5 |
double d21,d22; // d2 = Re(eta)
double Reflectivity1,Reflectivity2; // reflectivity = |b| * | eta +- (eta^2 - 1)^0.5 |^2
Delta_E1 = t - E_Bragg; // in eV
E_photon1 = t;
z01 = -2.0*pow(sin(theta_Bragg),2.0)*Delta_E1/E_photon1 + Sigma * F0;
eta1 = z01 / (Sigma * FH);
z11 = cpow(eta1,2.0+0i);
z21 = z11 - 1;
z31 = sqrt(z21);
d21 = eta1.m_re;
if (d21 < 0)
{
d11 = cabs(eta1+z31);
Reflectivity1= abs(b) * pow(d11,2.0);
}
else
{
d11 = cabs(eta1-z31);
Reflectivity1= abs(b) * pow(d11,2.0);
}
Delta_E2 = (t + dt) - E_Bragg; // in eV
E_photon2 = (t + dt);
z02 = -2.0*pow(sin(theta_Bragg),2.0)*Delta_E2/E_photon2 + Sigma * F0;
eta2 = z02 / (Sigma * FH);
z12 = cpow(eta2,2.0+0i);
z22 = z12 - 1;
z32 = sqrt(z22);
d22 = eta2.m_re;
if (d22 < 0)
{
d12 = cabs(eta2+z32);
Reflectivity2= abs(b) * pow(d12,2.0);
}
else
{
d12 = cabs(eta2-z32);
Reflectivity2= abs(b) * pow(d12,2.0);
}
F1 = Reflectivity1;
F2 = Reflectivity2;
//****************************
dInt = 0.5 * A * (R1*F1 + R2*F2) * dt ;
dIntegral += dInt;
t +=dt;
}while (t < 3.038); // ( dInt/dIntegral > dPrecision );
y=dIntegral + y0;
******* end of Origin C code ******
This code nicely fits an exemplary experimental data set so far. Now my questions are:
1. This fdf function is called from another Origin C - based program. Parameter initialisation is externally done in that other code which first of all imports the experimental data into Origin. Also the nine parameters which you find in the beginning of the code (E_Bragg, ..., b) are read in from a parameter file. Yet I did not succeed in handing them over from Origin C to the NLSF. This is why I defined them manually in the above given NLSF routine. That is very inconvenient for me as those parameters may change from experiment to experiment. Is there any chance to transfer those values from a parameter file to the NLSF routine (even via the other Origin C code)?
2. Eventually I will have to compute a convolution of three function and fit it to my experimental data. The experimental broadening R will still be a Gaussian while the theoretical curve will be a convolution of two theoretical curves already (i.e. two functions of type F as above). Can this be done with the trapezoidal rule also? Or are there any other options? Do you see any cance to compute such a convolution of three functions in a user-defined fitting function at all?
Your opinion will be highly accomplished. Thank you so much for dealing with those issues!
Best regards,
Oliver
|
Edited by - obauer on 05/27/2009 07:32:09 AM |
|
|
|
Iris_Bai
China
Posts |
Posted - 05/31/2009 : 06:37:21 AM
|
Hi Oliver,
For 1, for example the parameter E_Bragg in your code above, you can add this variable as function parameter and read value from your parameter file then assign to E_Bragg parameter in Parameter Initialization. The steps are:
1. Open Fitting Function Organizer in Tools menu.
2. Choose your function in left list.
3. In right panel, find Parameter Names and add E_Bragg.
4. Find Parameter Initialization and click Code Builder button to turn to Code Builder, now can add reading value from parameter file and assigned to E_Bragg.
Iris |
|
|
|
obauer
Germany
15 Posts |
Posted - 06/02/2009 : 11:01:59 AM
|
Hi Iris,
thank you for your fast and helpful reply. Meanwhile I have worked out something similar to what you suggested: I also defined E_Bragg ( and the other variables... ) as parameters - just as you suggested - and fixed them to a certian value which I passed on from the ORIGIN C code via the command:
NLSF.p5 = E_Bragg/1000.0; // E_Bragg is definded as double above in my code
NLSF.v5 = 0;
....
... and so on. E_Bragg is the 5th parameeter in my case and it is read in from a parameter file which is opened by the ORIGIN C routine. Yet if I understand you correctly you suggest to open the parameter file directly from the NLSF routine (*.fdf) in the section [Parameter Initialization]. How can I do that? Do you have an example at hand?
best regards and thank you so much,
Oliver |
|
|
|
Iris_Bai
China
Posts |
Posted - 06/03/2009 : 05:55:06 AM
|
Hi Oliver,
I did not know the structure of your file, you know, different file must use the different Origin C code to read. But the following give a sample to show the process about how to add Origin C code for parameter init in user defined function.
1. New a txt file under User File Folder, type "3,4,5" in this file and save with name "myfile.txt".
2. Type F9 to open Fitting Function Organizer, and new a function.
3. In Parameter Names editbox, and give parameters a1, a2, a2.
4. In Function edit box, change to y = x + a1 + a2 + a3;
5. In Parameter Initialization code edit box, click Code Builder Button turn to Code Builder window, copy the following code under "//Code to be executed to initialize parameters" line, and click Compile button, should get compiling Done in Output window, then click Return to Dialog button. Click Save button to save this function and close Fitting Function Organizer dailog.
string strFile = GetAppPath(false) + "myfile.txt";
if( strFile.IsFile() )
{
stdioFile ff(strFile, file::modeRead);
string strVal;
if( ff.ReadString(strVal) )
{
vector vVals;
strVal.GetTokens(vVals, ',');
a1 = vVals[0];
a2 = vVals[1];
a3 = vVals[2];
}
}
6. New a worksheet in Origin, and fill the two columns with row numbers, high light the two columns and open NLFit dialog.
7. Choose this new function, choose Parameter Tab, the value of a1, a2, a3 should be 3,4,5 just the same in myfile.txt file.
Iris |
Edited by - Iris_Bai on 06/03/2009 06:00:06 AM |
|
|
|
Stefan.E.S
Germany
11 Posts |
Posted - 11/17/2009 : 12:12:00 PM
|
Hi Oliver,
the Voigt function can be expressed analytically via the complex error function. I use the fit function below (Origin 7.5) which is much faster than a numerical integration.
Stefan
#include <origin.h>
//my own header, see below
#include "..\fitting\Fitfunctions.h"
//----------------------------------------------------------
//
void _nlsfMyVoigt(
// Fit Parameter(s):
double x0, double wL, double wG, double A,
// Independent Variable(s):
double x,
// Dependent Variable(s):
double& y)
{
// Beginning of editable part
const double sqrtln2 = 0.832554611;
const double sqrtpi = 1.772453851;
double Rez=2.0*sqrtln2*(x0-x)/wG; // the sign matters!
double Imz=sqrtln2*wL/wG;
double Rew, Imw;
complex_error_function(Rez,Imz,Rew,Imw);
y = A*2.0*sqrtln2/sqrtpi/wG*Rew;
// End of editable part
}
The complex error function itself is defined in a separate OriginC file Fitfunctions.c:
////////////////////////////////////////////////////////////////////////////////////
/**
* @brief Complex error function w(z)
*
* Calculates the complex error function (Faddeeva function) with relative error
* less than 10^(-r). r0=1.51*exp(1.144*r) and r1=1.60*exp(0.554*r) can be
* set by the the user subject to the constraints 14.88<r0<460.4 and 4.85<r1<25.5
*
* code taken from R. J. Wells, J. Quant. Spectrosc. Radiat. Transfer 62 (1999) 29-48
*
* @param x real part of the argument
* @param y imaginary part of the argument
* @param Re on exit real part of the complex error function
* @param Im on exit imaginary part of the complex error function
*/
void complex_error_function(double x, double y, double &Re, double &Im)
{
// const double r0 = 146.7, r1 = 14.67; // region boundaries for r=4
const double r0 = 460.4, r1 = 25.53; // region boundaries for r=5
// arguments
// x in input x value
// y in input y value >=0.0
// k out real part of w(z)
// l out imaginary part of w(z)
// constants
const double rrtpi = 0.56418958; // 1/sqrt(pi)
const double y0 = 1.5;
const double y0py0 = 3.0; // = y0+y0
const double y0q = 2.25; //= y0*y0
const double c[6]={1.0117281, -0.75197147, 0.012557727,
0.010022008, -0.00024206814, 0.00000050084806};
const double s[6]={1.393237, 0.23115241, -0.15535147,
0.0062183662, 0.000091908299, -0.00000062752596};
const double t[6]={0.31424038, 0.94778839, 1.5976826,
2.2795071, 3.0206370, 3.8897249};
// local variables
int j; // loop variable
int rg1, rg2, rg3; // y polynomial flags
double abx, xq, yq, yrrtpi; // |x|, x^2, y^2, y/sqrt(pi)
double xlim0, xlim1, xlim2, xlim3, xlim4; // |x| on region boundaries
double a0, d0, d2, e0, e2, e4, h0, h2, h4, h6; // W4 temporary variables
double p0, p2, p4, p6, p8, z0, z2, z4, z6, z8;
double b1, f1, f3, f5, q1, q3, q5, q7;
double xp[6], xm[6], yp[6], ym[6]; // CPF12 temporary values
double mq[6], pq[6], mf[6], pf[6];
double d, yf, ypy0, ypy0q;
//***** start of executable code *****************************************
rg1 = 1; // set flags
rg2 = 1;
rg3 = 1;
yq = y*y; // y^2
yrrtpi = y*rrtpi; // y/sqrt(pi)
// region boundaries when both k and l are required or when r!=4
xlim0 = r0 - y;
xlim1 = r1 - y;
xlim3 = 3.097*y - 0.45;
/* for speed the following 3 lines should replace the 3 above if r=4 and l is not required *
* xlim0 = 15100.0 + y*(40.0 + y*3.6); *
* xlim1 = 164.0 - y*(4.3 + y*1.8); *
* xlim3 = 5.76*yq; *
*/
xlim2 = 6.8 - y;
xlim4 = 18.1*y + 1.65;
if ( y <= 0.000001 ) // when y<10^-6
{
xlim1 = xlim0; // avoid W4 algorithm
xlim2 = xlim0;
}
abx = fabs(x); // |x|
xq = abx*abx; // x^2
if ( abx>xlim0 ) // region 0 algorithm
{
Re = yrrtpi / (xq + yq);
Im = Re*x/y;
}
else if ( abx > xlim1 ) // Humlicek W4 region 1
{
if ( rg1 ) // first point in region 1
{
rg1 = 0;
a0 = yq + 0.5; // region 1 y-dependents
d0 = a0*a0;
d2 = yq + yq - 1.0;
b1 = yq - 0.5;
}
d = rrtpi / (d0 + xq*(d2 + xq));
Re = d*y*(a0 + xq);
Im = d*x*(b1 + xq);
}
else if ( abx > xlim2 ) // Humlicek W4 region 2
{
if ( rg2 ) // first point in region 2
{
rg2 = 0;
h0 = 0.5625 + yq*(4.5 + yq*(10.5 + yq*(6.0 + yq))); // region 2 y-dependents
h2 = -4.5 + yq*(9.0 + yq*( 6.0 + yq* 4.0));
h4 = 10.5 - yq*(6.0 - yq* 6.0);
h6 = -6.0 + yq* 4.0;
e0 = 1.875 + yq*(8.25 + yq*(5.5 + yq));
e2 = 5.25 + yq*(1.0 + yq* 3.0);
e4 = 0.75*h6;
f1 = -1.875 + yq*(5.25 + yq*(4.5 + yq));
f3 = 8.25 - yq*(1.0 - yq* 3.0);
f5 = -5.5 + yq* 3.0;
}
d = rrtpi / (h0 + xq*(h2 + xq*(h4 + xq*(h6 + xq))));
Re = d*y*(e0 + xq*(e2 + xq*(e4 + xq)));
Im = d*x*(f1 + xq*(f3 + xq*(f5 + xq)));
}
else if ( abx < xlim3 ) // Humlicek W4 region 3
{
if ( rg3 ) // first point in region 3
{
rg3 = 0; // region 3 y-dependents
z0 = 272.1014 + y*(1280.829 + y*(2802.870 + y*(3764.966
+ y*(3447.629 + y*(2256.981 + y*(1074.409 + y*(369.1989
+ y*(88.26741 + y*(13.39880 + y)))))))));
z2 = 211.678 + y*(902.3066 + y*(1758.336 + y*(2037.310
+ y*(1549.675 + y*(793.4273 + y*(266.2987
+ y*(53.59518 + y*5.0)))))));
z4 = 78.86585 + y*(308.1852 + y*(497.3014 + y*(479.2576
+ y*(269.2916 + y*(80.39278 + y*10.0)))));
z6 = 22.03523 + y*(55.02933 + y*(92.75679 + y*(53.59518
+ y*10.0)));
z8 = 1.496460 + y*(13.39880 + y*5.0);
p0 = 153.5168 + y*(549.3954 + y*(919.4955 + y*(946.8970
+ y*(662.8097 + y*(328.2151 + y*(115.3772 + y*(27.93941
+ y*(4.264678 + y*0.3183291))))))));
p2 = -34.16955 + y*(-1.322256+ y*(124.5975 + y*(189.7730
+ y*(139.4665 + y*(56.81652 + y*(12.79458
+ y*1.2733163))))));
p4 = 2.584042 + y*(10.46332 + y*(24.01655 + y*(29.81482
+ y*(12.79568 + y*1.9099744))));
p6 = -0.07272979 + y*(0.9377051+ y*(4.266322 + y*1.273316));
p8 = 0.0005480304 + y*0.3183291;
q1 = 173.2355 + y*(508.2585 + y*(685.8378 + y*(557.5178
+ y*(301.3208 + y*(111.0528 + y*(27.62940
+ y*(4.264130 + y*0.3183291)))))));
q3 = 18.97431 + y*(100.7375 + y*(160.4013 + y*(130.8905
+ y*(55.88650 + y*(12.79239+y*1.273316)))));
q5 = 7.985877 + y*(19.83766 + y*(28.88480 + y*(12.79239
+ y*1.909974)));
q7 = 0.6276985 + y*(4.264130 + y*1.273316);
}
d = 1.7724538 / (z0 + xq*(z2 + xq*(z4 + xq*(z6 + xq*(z8+xq)))));
Re = d*(p0 + xq*(p2 + xq*(p4 + xq*(p6 + xq*p8))));
Im = d*x*(q1 + xq*(q3 + xq*(q5 + xq*(q7 + xq*0.3183291))));
}
else // Humlicek cpf12 algorithm
{
ypy0 = y + y0;
ypy0q = ypy0*ypy0;
Re = 0.0;
Im = 0.0;
for (j = 0; j <= 5; j++)
{
d = x - t[j];
mq[j] = d*d;
mf[j] = 1.0 / (mq[j] + ypy0q);
xm[j] = mf[j]*d;
ym[j] = mf[j]*ypy0;
d = x + t[j];
pq[j] = d*d;
pf[j] = 1.0 / (pq[j] + ypy0q);
xp[j] = pf[j]*d;
yp[j] = pf[j]*ypy0;
Im = Im + c[j]*(xm[j]+xp[j]) + s[j]*(ym[j]-yp[j]);
} // end for
if ( abx < xlim4 ) // Humlicek CPF12 region i
{
for (j = 0; j <= 5; j++)
{
Re = Re + c[j]*(ym[j]+yp[j]) - s[j]*(xm[j]-xp[j]);
}
}
else // Humlicek CPF12 region ii
{
yf = y + y0py0;
for (j = 0; j <= 5; j++)
{
Re = Re
+ (c[j]*(mq[j]*mf[j]-y0*ym[j]) + s[j]*yf*xm[j]) / (mq[j]+y0q)
+ (c[j]*(pq[j]*pf[j]-y0*yp[j]) - s[j]*yf*xp[j]) / (pq[j]+y0q);
}
Re = y*Re + exp ( -xq );
}
}
}
The function header is passed in Fitfunctions.h:
/**
* @brief Complex error function w(z)
*/
void complex_error_function(double x, double y, double &Re, double &Im);
|
|
|
| |
Topic |
|
|
|
About editing the square of Voigt function