MATLAB: Plotting error for Bessel Function for “large” input values.

besselgaussianplot

Sample of code in question:
%% Inputs
InputAngle = 4;
Wavelength = 0.65;
z = 10100;
FiberCoreRadius = 100;
%% Window
r = 0:10:10000;
%% Dependents
k = 2*pi/Wavelength;
AngleRad= InputAngle/180*pi;
w0 = FiberCoreRadius*0.67; % Assuming fiber is long enough to generate
Zr = k*w0^2/2;
wz = w0*sqrt(1+z^2/Zr^2);
Rz = (Zr^2+z^2)/(z);
alpha = k*sin(AngleRad);
l = alpha*FiberCoreRadius;
%% BG Code
Gaussian = exp(-(1/wz^2-1i*k/(2*Rz)).*(r.^2+alpha^2*z^2/k^2));
Bessel = besselj(l,(alpha.*r)./(1+1i*z/Zr));
FullBeamProfile = abs(Gaussian.*Bessel).^2;
%% Plotting
hold on;
plot(r,FullBeamProfile/max(FullBeamProfile),'b');
For larger Z values, the plot is returning zero for all values and not even attempting to plot larger values of r which is preventing me from modeling what I need to for my research. For smaller Z values there is no issue, but I need larger values to confirm what I'm working on.
The posted example gives the gaussian curve but does not go the full 10000 'r' distance I requested which becomes problematic when I need to see beyond that (for larger z values). It looks as though my Bessel function is going beyond the limits of a 64bit float and returning INF, and my gaussian function is returning below the limit, thus returning 0. I'm not 100% sure if this is causing this issue though.
Does anyone have any idea on how to fix this kind of error or have any suggested work arounds? I would love to have the full Gaussian curve generated, but if it helps, having the distance to peak value would work as well.

Best Answer

Hi Greg,
Time to look at logs. The bessel function has a scaled version that keeps it from blowing up, and you can look at the log of the whole works as is done below. The plot shows the (complex) log of (Gaussian*Bessel)^2 and is free of overflow or underflow problems. Exponentiating that quantity and taking abs will get rid of the imaginary part at which point only the real part matters. The real part has a maximum at about r = 700, but it drops down to around -3e4 at r = 10000. exp(-3e4) is so ridiculously small that I don't think this a double precision problem, really. It's more an issue of what is even physically meaningful.
%% Inputs
InputAngle = 4;
Wavelength = 0.65;
z = 10100;
FiberCoreRadius = 100;
%% Window
r = 0:10:10000;
%% Dependents
k = 2*pi/Wavelength;
AngleRad= InputAngle/180*pi;
w0 = FiberCoreRadius*0.67; % Assuming fiber is long enough to generate
Zr = k*w0^2/2;
wz = w0*sqrt(1+z^2/Zr^2);
Rz = (Zr^2+z^2)/(z);
alpha = k*sin(AngleRad);
l = alpha*FiberCoreRadius;
% -------------same up to this point ----------------
%% BG Code
G = exp(-(1/wz^2-1i*k/(2*Rz)).*(r.^2+alpha^2*z^2/k^2));
logG = (-(1/wz^2-1i*k/(2*Rz)).*(r.^2+alpha^2*z^2/k^2));
% Bessel = besselj(l,(alpha.*r)./(1+1i*z/Zr));
argB = (alpha.*r)./(1+1i*z/Zr);
logsfac = abs(imag(argB));     % scaled bessel is bessel divided by exp(logsfac) 
logBscaled = log(besselj(l,(alpha.*r)./(1+1i*z/Zr),1));       % note third input
%FullBeamProfile = abs(Gaussian.*Bessel).^2;
logFBP = 2*(logG +logsfac +logBscaled);
figure(1)
plot(r,real(logFBP),r,imag(logFBP))
return
Related Question