Problem on integral and ploting...

Hello there!

I do have a problem with symbolic integration. The original code looks like this:

function f = fmbloss(phi)
n1 = 1.491;
n2 = 1.405;
ns1 = n1*n1;
ns2 = n2*n2;
lambda = 0.58e-6;
ra1 = 0.5e-3;
kw = 2*pi/lambda;
bend_r = 10e-3;
thetac = acos(n2/n1);
rad = @(theta) (bend_r+ra1)*n1./(n2*cos(theta));
n1_phi = @(theta) (sqrt(ns1-ns1*cos(theta).*cos(theta)));
n2_phi = @(theta) (sqrt(ns1*cos(theta).*cos(theta)-ns2));
p = @(theta) (rad(theta)/(bend_r+ra1));
ps = @(theta) (p(theta).*p(theta));
ps1 = @(theta) (sqrt(ps(theta)-1));
log_p = @(theta) (log(p(theta)+ps1(theta)));
abs_t = @(theta) (4*n1_phi(theta).*n2_phi(theta)/(ns1-ns2));
expf = @(theta) (-2*kw*ra1*n1*cos(theta).*(log_p(theta)-ps1(theta)./p(theta)));
bt = @(theta) (abs_t(theta).*expf(theta));
sub_sin = @(theta) (sqrt(sin(theta).*sin(theta)-sin(thetac)*sin(thetac)));
sub_sins = @(theta) (sin(theta)+sub_sin(theta).*sub_sin(theta));
bt1 = @(theta) (4*sin(theta).*sub_sin(theta)./(sub_sins(theta)));
exp_f1 = @(theta) exp(-bt(theta)*phi./(2*theta));
exp_f2 = @(theta) exp(-bt(theta)*bend_r*phi.*theta/(4*ra1));
exp_f3 = @(theta) exp(-bt1(theta)*bend_r*phi/(4*ra1));
fun1 = @(theta) (exp_f1(theta));
fun2 = @(theta) (exp_f2(theta));
fun3 = @(theta) (exp_f3(theta));
theta_max = sqrt(2*ra1/bend_r);
int1 = quadl(fun1,0,theta_max);
int2 = quadl(fun2,theta_max,thetac);
int3 = quadl(fun3,thetac,pi/2);
f = pi*(int1+int2+int3)/2;

and input fmbloss(pi/3) in the command window firstly,

and error message：

>> fmbloss(pi/3)
Warning: Infinite or Not-a-Number function value encountered. 
> In quadl at 105
In fmbloss at 40 
Warning: Infinite or Not-a-Number function value encountered. 
> In quadl at 105
In fmbloss at 41 
ans =

   Inf

input fmbloss(0) in the command window secondly,

>> fmbloss(0) ans = 2.467401100272340

It is not working at all.I hope someone can give me a good advise, how to deal with that. Thank you!

[Walter Roberson]

Your routine contains no symbolic integration ?

You are hitting infinity in your integration. You should learn how to use the debugger to figure out why. In particular,

dbstop if naninf