non-analytical solutions

Share icon

Wed, 01/02/2013 - 06:52 am

I am analyzing a lot of time of flight data and having an issue with the analysis. I am modeling the data with a feature that combines Gaussian and Lorentzian distributions. In addition, I can add a bit of Skew to account for the tail in the data. But basically the function is best described the Intensity as a function of time equals Gaussian + Lorentzian. It is easy to determine the maximum intensity and the time associated with that maximum. However, I need to determine the time at 90%, 50%, 40%, and 10% of that maximum, but I cannot solve the equation for time analytically.

I assume that I can give IGOR the following (a function, the fit parameters, and a given Intensity) and then ask it to ask it to return to me the Time at 0.9*maximum, 0.5*maximum, 0.4*maximum, and 0.1*maximum. Can someone help me?

Here is the function I am using. Also included is the fit parameters for a feature, and the representative feature that I am trying to work with. Thanks!!!!!




Function XES_Fit1(param, xx)        // basis: SAK2_Gauss, added lorentzian functions
    // Param [0] XES_Fit1   Y0 (offset)
    // Param [1] XES_Fit1   A (amplitude)
    // Param [2] XES_Fit1   x0 (peak position)
    // Param [3] XES_Fit1   half width (width)
    // Param [4] XES_Fit1 mixing (Gaussian : Lorentzian)
    // Param [5] XES_Fit1   Skew
    // Param [6] XES_Fit1   LW-ratio
    wave param; variable xx
    variable tempA
        tempA = param[0]+param[1]*(param[4]*ga1(param,xx)+(1-param[4])*lre1(param,xx))  
    return tempA  
end

    Function ga1(param,xx)
        Wave param; Variable xx
        Variable tempB
            tempB= exp((-ln(2)*(xx-param[2])^2)/(((param[3]+(xx-param[2])*param[5])*param[6])^2))
        return TempB
    end
   
    Function lre1(param,xx)
        Wave param;Variable xx
        Variable TempC
            TempC=(param[3]+(xx-param[2])*param[5])^2/((param[3]+(xx-param[2])*param[5])^2+(xx-param[2])^2)
        return tempC
    end







Here are the parameters for the following feature:
0
-8.17597
7.29073e-06
2.26073e-08
0.207583
0.271577
1








Intensity Time (mSec)
0 -1.4993e-06
0 -1.4983e-06
0 -1.4973e-06
0 -1.4963e-06
0 -1.4953e-06
[the rest of the data was removed by aclight because it was too long]
@PiStosHiO17:
I think the extreme size of the data you included in your original post was preventing the post from being displayed. I have edited your original post to remove all but the first few lines of the data. If the data is important, please attach it as a file.

January 2, 2013 at 07:36 am - Permalink

If you can't find an analytic solution for your fit function, you can use the FindRoots operation to find where the function has certain Y values. Use FindRoots/Z=(0.1*W_coef[1]) to find the X value where Y is at 10 percent of the amplitude.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 2, 2013 at 11:07 am - Permalink

You mention that your fit "combines Gaussian and Lorentzian distributions". Depending on how this combination arises, you might also want to consider use of the Voigt function, which is a convolution of a Gaussian and Lorentzian. One simple way to access this fitting function is via the menus Analysis->Packages->Multipeak Fitting->Multipeak Fitting 2.

Read the associated help file Multi-peak Fitting 2 Help.ihf, and try the demo experiment for further details.

January 2, 2013 at 12:34 pm - Permalink

And the Multipeak Fit 2 package includes a batch-fit function that you might be able to use as the basis of what you're doing. Steve's right- the sum of Gaussian and Lorenzian is one method that's sometimes used to approximate a Voigt in order to avoid the convolution. Multipeak Fit uses a different, more accurate approximation.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 2, 2013 at 04:55 pm - Permalink

Thanks for the help

When I use the FindRoots option it gives me the solutions in the command window. For an assymetric function is there a good way to store one of the solutions as a point in a wave?

PiSTosHiO
Stosh Kozimor

January 3, 2013 at 04:40 am - Permalink

Use /Q to suppress the history report.

Access the results via the variables that FindRoots creates. Since your function is univariate, use the variables for a single 1D funtion. From the FindRoots documentation:

The results of finding roots of a single 1D function are put into several variables:
V_numRoots Number of roots found. Either 1 or 2.
V_Root The root.
V_YatRoot Y value of the function at the root. *Always* check this; some discontinuous functions may give an indication of success, but the Y value at the found root isn't even close to zero.
V_Root2 The second root if FindRoots found two roots.
V_YatRoot2 Y value at the second root.

Wrap FindRoots in a user-defined function; you can access the roots from the output variables and assign them to a wave after FindRoots returns. Be sure your code checks V_flag to make sure the root finder didn't encounter obvious problems. Be sure to check V_numRoots to decide if you need to run a second FindRoots for your asymmetric peaks. Be sure to check V_YatRoot (and V_YatRoot2) to make sure that the found root is really a root. As long as the /H and /L brackets are set reasonably, your function is well-behaved and you shouldn't have problems with finding roots that aren't roots.

Read the help file about how FindRoots works- copy this command, paste it into the command line and press Enter:

DisplayHelpTopic "Finding Function Roots"

You can skip the section on polynomial roots :)

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 3, 2013 at 09:12 am - Permalink