Numerical integration of tabulated data

5 views (last 30 days)
Jonathan on 7 Sep 2017
Commented: Jonathan on 8 Sep 2017
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 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.

Answers (1)

Steven Lord
Steven Lord on 7 Sep 2017
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.
  1 Comment
Jonathan on 8 Sep 2017
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?

Sign in to comment.


Find more on Particle & Nuclear Physics in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!