how can I pass the points in which a funtion is evaluated in an integral to a different one?

Hi to everyone, as the title sais my question is:

how can I pass the points in which a funtion is evaluated in an integral to a different one?

I am trying to evaluate the value of the absorption coefficient of a fluorescent slab. The idea is to fit experimental data whit a model and the data correspond to the intensity registered by an optical fibre with the source at different distances.

The problem is that the source is a spot composed by two components: the spot due to the primary light and the ring due to the fluorescence light.

The contribution due to the primary light is quite simple and correspond to the integral on the spot integrated using integral2. The integrand function is:

I_0*eta*nems*e^(-abs*sqrt((dist-x)^2+y^2))*(dist-x))/(2*pi*b*h*((dist-x)^2+y^2))

Where I made the assumption that the zero of my reference frame is at the spot centre and I_0 is the initial intensity, abs is the absorption coefficient, dist the distance of the spot centre from the optical fibre (that here is considered a point-like object for simplicity), nems is the normalized emission spectrum and eta the emission efficiency at each wavelength. x and y correspond to the coordinates of a point I the primary spot, and they range from range [-b,b] and [0,h] respectively.

The contribution due to the fluorescent ring can be treated as before since it can be written as:

I_F*eta*nems*exp(-abs*sqrt((dist-x_a^2+y_a^2))*(dist-x_a))/(2*pi*b*h*((dist-x_a)^2+y_a^2))

where now x_a and y_a correspond to the coordinates of a point I the fluorescence ring, and they range from range [b,-ln(0.1)/abs] and [h,-ln(0.1)/abs] respectively.

The problem is that I_F is the intensity of each point of the fluorescence ring due to the primary spot, it has to be computed as the integral3 of:

(I_0*eta*nems*exp(-abs*sqrt((x-x_p)^2+(y-y_p)^2)))/(h*b)

evaluated on the primary spot so with x = [-b,b] and y = [0,h]. The third dimension is due to the integral on all the wavelength lambda = [lambdamin, lambdamax].

So, exist a way thanks to which these tree integral functions can exchange their domain points (x,y) and (x_p,y_p)? Or I must write an integral function on my own in which I can control the points in which the function is evaluated?

I am sorry I don't add any code, but I have no idea how to do this. Maybe I am just to worked up in the problem and can't see the wood for the trees. Help is greatly appreciated.

Thanks in advance and best regards,

Giulio

Answer 1

[Alan Stevens]

% I’m not sure I fully understand the system you are trying to model,
% but how about using a Monte-Carlo approach to the integrals?
%
% Begin by randomly scattering M points over the spot. 
% Each of these points, (x,y), represents on average a “differential” area
% of dxdy = Aspot/M (where Aspot is the area of the spot).
% Approximate your first integral by summation:
% Ispot = sum(I_0*eta*nems*e^(-abs*sqrt((dist-x)^2+y^2))*(dist-x))/(2*pi*b*h*((dist-x)^2+y^2)))*dxdy;
%
% Now scatter N random points over the fluorescent ring.
% Each of thee points (xp,yp) represents on average a "differential" area
% of dxpdyp = Aring/N. Approximate your third integral by
% I_F(xp,yp) =
% sum((I_0*eta*nems*exp(-abs*sqrt((x-x_p)^2+(y-y_p)^2)))/(h*b))*dxpdyp
% where the summation is over the M values of the (x,y) points in the spot.
% There will be N values of I_F, one for each of the (xp,yp) points.
%
% Now you can approximate your second integral by summing 
% I_Fcontribution  = 
% sum(I_F(xp,yp)*eta*nems*exp(-abs*sqrt((dist-xp^2+yp^2))*(dist-xp))/(2*pi*b*h*((dist-xp)^2+yp^2)))*dxpdyp
%
% The larger you make N and M, the more accurate your integrals will be,
% though the greater the number of calculations, of course!
%
% Just a suggestion!  It might or might not make sense in the context of
% your system.

I meant to add that it isn't a good idea to use abs as a variable name, as it is a built-in function (for absolute value).

[Giulio Mangherini]

Thank you for the answer and for the add on.

A Monte-Carlo approach could make sense in my context, but I have a problem with it.

After that I computed the various integrals,I want to fit them with my experimental data.

So calling: intensity_correction = I_S + I_R

where I_S and I_R are the integral computed with your suggestion for the spot and the ring respectively. I use "lsqcurvefit" to fit them as following:

for i=1:size_red(1)
red_line = red(i,:);    
[coeff,~,R,~,~,~,J] = lsqcurvefit (@intensity_correction,coeff,distances,red_line,[1000 0]);
temp = nlparci(coeff,R,'jacobian',J);
err_coeff (i,:) = [temp(1) temp(3) temp(2) temp(4)]; 
fit_coeff(i,:) = coeff;
end

where:

• "coeff" is a vector 1x2 with coeff(1) = I_0 and coeff(2) = abs_coeff (instead of "abs")
• "red" is the matrix of my experimental data having on the row the different wavelength and on the colum the different distances
• "distances" is the vector with the various distances

The purpose of this for loop is to find the best I_0 and abs_coeff at all distances but for one wavelength, and then cycle all the wavelengths: Where I am using "lsqcurvefit" instead of "nlinfit" abs_coeff must be positive.

The problem is that "lsqcurvefit" does not converge for most of the wavelength.

Maybr is becuse I am generating new sets of point (x,y) and (xp,yp) every time that "lsqcurvefit" calls "intensity_correction", but I am not sure it is just guess.

[Alan Stevens]

Firstly, the kernels of your integrals can be expressed in the form

k*e^(sqrt((dist-x)^2+y^2))*(dist-x))/(2*pi*b*h*((dist-x)^2+y^2)));

where

k = I_0*eta*nems*exp(-abs_coeff); 

(with analogous expressions for the others).

This means you don't need to include k when doing the Monte_Carlo stuff. Just do the Monte-Carlo on the

e^(sqrt((dist-x)^2+y^2))*(dist-x))/(2*pi*b*h*((dist-x)^2+y^2)))

type expressions. This can be done "off-line", as it were, independently of the experimental data. Then when you fit to experimental data, these integrals are just constant multipliers on the terms in k (I_0, abs_coeff, etc.) to which you are trying to find the best fit values.

Secondly, make sure you try a few different sets of values for M and N (by increasing them) to ensure you have integrals that have converged to within acceptable tolerances.