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
回答（1 个）
Titus Edelhofer
2011-8-27
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
2 个评论
Mikheil
2011-8-27
Titus Edelhofer
2011-8-27
You mean, z should be defined? Yes, something like z=zeros(size(f));
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
