Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to integrate a spectrum I obtained from a fourier transform infrared spectroscopy experiment, in order to determine the global warming potential of green house gases.

In order to do this, I would like to integrate the spectrum that I have in chunks of 10 wave numbers (x axis is wave numbers). IE, I would like to generate a list that has the integral from 525 - 535 in the first slot, 535-545 in the next slot

Here is a link to my data: https://drive.google.com/file/d/0B5gGZ4-SD7-ZZDRoNE5XMzFnYnM/edit?usp=sharing

data = Import["/home/marco/Documents/Mathematica/GlobalWarming/nitrogenforMathematicaTest.csv"]
spectraInterp = Interpolation[data, InterpolationOrder->3] (* Here we  create an interpolation function to represent our data *)


ListPlot[data]
Plot[spectraInterp[x], {x,525,1400}, PlotRange->{{525, 1400}, {0, 0.3}}]
NIntegrate[spectraInterp[x], {x,600,1400}]
share|improve this question
    
Your data runs in {{0, 7}, {16000, 22000}} ... Where are those wavelenghts? –  belisarius Nov 21 '13 at 19:37
    
@belisarius, for an FTIR dataset I would expect the X-axis to be in units of wavenumbers, which are proportional to frequency rather than wavelength. But I don't see why the first column of data runs between 0 and 7 either. What units are they suppose to be @olliepower? The main problem I see is that your interpolation function is valid for values between 0 and 7, but you want to integrate between 600 and 1400. –  Jason B Nov 21 '13 at 19:41
    
@JasonB What you say is correct. I wrote this bit of code using a different data set, and neglected to change my range of integration. I want to integrate between 525 and 1400 wave numbers. –  olliepower Nov 22 '13 at 6:00
    
Are you sharing the right data? The filename in the google-drive-linked notebook references gas chromatography. –  bobthechemist Dec 3 '13 at 2:48
    
@bobthechemist Here is the correct spectrum. –  olliepower Dec 3 '13 at 23:00
add comment

1 Answer

up vote 1 down vote accepted

Since there is some concern about what the interval of integration should be, why not use a general procedure like this to create all the integrals you need?

bIntegrate[f_, {x_, a_, b_}, dx_] := Block[{x0, n},
    n = Floor[(b - a)/dx];
    x0 = Ceiling[a/dx] dx;
    Table[Integrate[f, {x, x0 + k dx, x0 + (k + 1) dx}], {k, 0, n - 1}]
    ]

It will split the integral into subintegrals whose extrema are subsequent multiples of dx. You will lose something at the beginning and at the end of the interval if a and b are not multiples of dx themselves. This is a necessary measure to keep the interval length constant over all integrals.

Toy example

bIntegrate[ f[x], {x, 427, 446}, 5 ]

will give you a list with integrals on the intervals {430-435}, {435-440} and {440,445}.

In your case, with the data ranging from 525 to 1400 you will invoke it as

bIntegrate[ spectralInterp[x], {x, 525, 1400}, 10 ]

The integrals in the returned list will start from 530 (the closest multiple of 10 contained in the interval of integration) and will end in 1400 (likewise) with intervals of integration of length 10.

EDIT: I see that you used NIntegrate. Just change Integrate in NIntegrate in the above function, then. A better name would then be biNIntegrate. :-) If you want to start from 525, modify the code in this way (just remove x0 and use a instead - even turn the procedure into a one-liner)

biNIntegrate[f_, {x_, a_, b_}, dx_] := 
    Table[NIntegrate[f, {x, a + k dx, a + (k + 1) dx}], {k, 0, Floor[(b - a)/dx]- 1}]

The numerical integrals will be over {525,535}, {535,545}, {545,555}, ..., {1385,1395}.

Further improvement: you can add some checking on the arguments to make sure that f is indeed a function of x, the extrema are indeed numerical and so on...

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.