dblquad integration problem

I am trying to do double integral using dblquad function, it seems that everything works fine the only problem is I am setting an "if" condition in function window and the program does not take it into consideration it still goes beyond the condition I have set. here is the function code: lcs.m

function z = lcs(E_gamma1, phi,E_gamma) 

gamma = 85;   % Lorentz Factor
El = 0.00466106; % keV
Em = 4*(gamma^2)*El; % Maximum X-ray energy
%E_gamma = Em/(1+(gamma^2)*(theta^2));
theta = (sqrt(-1+Em./E_gamma1))*(1/gamma);    % scattering angle
%theta = theta/10;
theta_d = 0;
phi_d = 0;
theta_xd = theta_d*cos(phi_d);
theta_yd = theta_d*sin(phi_d);
theta_x = theta*cos(phi);
theta_y = theta*sin(phi);
sigma_e = 100; %5.12E-3;
sigma_gam = (2*Em*sigma_e)/(gamma*(1+(theta_d^2)*(gamma^2))^2);
sigma_x = 5; %0.002;
sigma_y = 5; %0.003;
f = (1+cos(2*phi)).*((E_gamma1)/Em).*(((E_gamma1)/Em)-1)+0.5;

if (E_gamma1 > Em) | (E_gamma1 < El)

    f = 0;
else

    z = f.*(exp(-((E_gamma1-E_gamma).^2)/(2*sigma_gam^2))).*(exp(-((theta_xd-theta_x).^2)/(2*sigma_x^2))).*(exp(-((theta_yd-theta_y).^2)/(2*sigma_y^2)));
end

And this is the dblquad code: integr.m

clc
clear all
for E_gamma = 1:250
Q(E_gamma) = dblquad(@lcs,10,140,0,2*pi,[],[],E_gamma);
end
plot(Q, '.')

I would appreciate your help

   -M

Answer 1

[Titus Edelhofer]

Hi,

first: I guess in the "if" it should be z=0? Second: your function should allow for vector inputs. Therefore the if should be something like

idx = (E_gamma1 > Em) | (E_gamma1 < El);
z(idx) = 0;
z(~idx) = f(~idx).* ...

Titus

[Mikheil]

Thanks for your reply,

when I used your suggestion about idx, it should be defined otherwise it gives me errors.

any suggestion?

[Titus Edelhofer]

You mean, z should be defined? Yes, something like z=zeros(size(f));