Numerical integration of tabulated data

Hello,

I'm trying to follow the tutorial here so that I can generate a method of drawing a random number (between 0 and 1) to determine the scattering angle of a photon in a Monte Carlo simulation (of light propagating through water).

My phase function (beta tilde), scattering per Steradian vs. angle, is stored in a Matlab vector, betaFF. From the tutorial, a probability density function (PDF) is required to generate scattering angles from the random number (Tau)..." to determine the polar scattering angle, we draw a U[0,1] random number (Tau) and solve":

Furthermore: "In general (this equation) must be solved numerically because of the complicated shape of most phase functions, or when the phase function is defined by tabulated data at a finite number of scattering angles and is fit with a spline (or other) function to generate a continuous function of scattering angle."

A typical output, which is what I'm aiming to obtain, is shown at the bottom of that tutorial page:

I can fit the tabulated betaFF values using Matlab's "fit" function and 'linearinterp', as described in the tutorial (tabulated data at a finite number of scattering angles...is fit with a spline (or other) function):

func = fit(psi,betaFF,'linearinterp');

where psi is a vector containing angles from 0 to 180 degrees.

From there, I get pretty lost. I don't know how I'm supposed to use this linear fit (func) to numerically evaluate the equation above. Is it possible to create a function in Matlab containing both "func" and the sin(psi) term and integrate that? Any pointers would be very welcome.

Answer 1

[Steven Lord]

Since you've fitted a function you can use the integrate function on the func object. If you wanted to numerically integrate a vector or matrix of data, use trapz or cumtrapz as shown in this example.

[Jonathan]

Hello,

Thanks for your response. Am I correct in saying that the "integrate" function only works for a cfit object?

My issue is that I need to integrate the product of a cfit object (found using 'linearinterp' and my tabulated data for Beta tilde) and the function sin(psi). I don't think I can do this using the "integrate" function?

I have tried to define a function:

func = @(Theta) sin(Theta).*fit(Theta);

Where "fit" is my 'linearinterp' fit.

And then integrate that using the integral command:

q = integral(fun,0,pi)

But I get the error message:

Error using  .* 
Matrix dimensions must agree.

Error in @(Theta) sin(Theta).*fit(Theta)

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

I don't understand why that doesn't work. The dimensions of Theta, sin(Theta), fit(Theta) and betaFF match (all are 1800x1 doubles), and if I simply calculate the product the result is also a 1800x1 double. Why do I get the "matrix dimensions must agree" error?