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
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?
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.
Script:
### 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 script
Result:
Addition: (irregular spacing in x)
In case your x-data is not spaced by 1 and not even equidistant you would have to adjust the procedure to calculate the area under a curve.
The area under the curve between two datapoints (x0,y0)
and (x1,y1
is simply (x1-x0)*(y0+y1)*0.5
.
So, you can calculate the area under the curve during plotting the original data.
In the example below, the plotting line with impulses
is just there to better illustrate the different spacings of the datapoints.
Script:
### calculate area under curve for irregular x-spacings
reset session
# create some random test data
set table $Data
N = 15
set samples N
plot x0=0 '+' u (x0=x0+rand(0)*10):(100*exp(-((x0-N*2)/N)**2)) w table
unset table
set style fill solid 0.3 border
set key noautotitle
plot auc=0 $Data u ($0==0?(x1=$1,y1=$2):0, y0=y1, y1=$2, x0=x1, x1=$1): \
(auc=auc+(x1-x0)*(y0+y1)*0.5,y1) w lp pt 7 ps 0.5 lw 2, \
'' u 1:2 w impulses lc "red", \
'' u 1:2 w filledcurves y=0 lc rgb 0xccff0000 ti sprintf("Area: %g",auc)
### end of script
Result: