Search code examples
gnuplotdata-fitting

Calculate area under curve in gnuplot


I have data (data can be downloaded here: gauss_data) and need to find the area of a particular peak. From my data set, the peak seems to have some contribution from another peak. I made the fit on my data with 3 Gaussians using this code:

# Gaussian fit                                              

reset 
set terminal wxt enhanced

# Set fitting function      

f(x)  = g1(x)+g2(x)+g3(x)
g1(x) = p1*exp(-(x-m1)**2/(2*s**2))
g2(x) = p2*exp(-(x-m2)**2/(2*s2**2))
g3(x) = p3*exp(-(x-m3)**2/(2*s3**2))


# Estimation of each parameter      

p1 = 100000
p2 = 2840
p3 = 28000
m1 = 70
m2 = 150
m3 = 350
s = 25
s2 = 100
s3 = 90

# Fitting & Plotting data                           

fit [0:480] f(x) 'spectrum_spl9.txt' via p1, m1, s, p2, m2, s2, p3, m3, s3
plot [0:550] 'spectrum_spl9.txt' lc rgb 'blue', f(x) ls 1, g1(x) lc rgb 'black', g2(x) lc rgb 'green' , g3(x) lc rgb 'orange' 

and the result is shown in fig below

enter image description here

I need to calculate the area under the peak i.e. area f(x) - area g3(x). Is there any way to find/calculate the area of each function in Gnuplot?


Solution

  • Your data is equidistant in x-units with a step width of 1. So, you can simply sum up the intensity values multiplied by the width (which is 1). If you have irregular data then this would be a bit more complicated.

    Code:

    ### determination of area below curve
    reset session
    FILE = "SO/spectrum_spl9.txt"
    
    # fitting function      
    f(x)  = g1(x)+g2(x)+g3(x)
    g1(x) = p1*exp(-(x-m1)**2/(2*s1**2))
    g2(x) = p2*exp(-(x-m2)**2/(2*s2**2))
    g3(x) = p3*exp(-(x-m3)**2/(2*s3**2))
    
    # Estimation of each parameter      
    p1 = 100000
    p2 = 2840
    p3 = 28000
    m1 = 70
    m2 = 150
    m3 = 350
    s1 = 25
    s2 = 100
    s3 = 90
    
    set fit quiet nolog
    fit [0:480] f(x) FILE via p1, m1, s1, p2, m2, s2, p3, m3, s3
    
    set table $Difference
        plot myIntegral=0 FILE u 1:(myIntegral=myIntegral+f($1)-g3($1),f($1)-g3($1)) w table
    unset table
    
    set samples 500    # set samples to plot the functions
    plot [0:550] FILE u 1:2 w p lc 'blue' ti FILE noenhanced, \
          f(x) ls 1, \
          g1(x) lc rgb 'black', \
          g2(x) lc rgb 'green', \
          g3(x) lc rgb 'orange', \
          $Difference u 1:2 w filledcurves lc rgb 0xddff0000 ti sprintf("Area: %.3g",myIntegral)
    ### end of code
    

    Result:

    enter image description here