v = (0.01*2.0e-12) ;u = 0.036;
K = 8.617e-5;
T = 300; hBar = 1;
no = 10e16; dg = 1;
trig = (sin(4))^2;
t = 2.0e-12; e = 1.0;
deltag = 0.0156; deltaxi =0.0204;
w = [0.2 0.5 5.0];
Betaxi = linspace(0,10, 201);
p = (0.9*pi);
delta1 = deltag./(K.*T); delta2 = deltaxi./(K.*T);
I0 = besseli(0, delta1); I1 = besseli(1, delta1);
I3 = besseli(0, delta2); I4 = besseli(1, delta2);
I2 = (I0./I1); I5 = (I3./I4);
joxi = ((no*e*deltag*dg)/hBar).*I5;
[mu] = meshgrid(-10000:100:10000);
legendStrings = cell(length(w), 1);
for k1 = 1:length(w)
thisN = w(k1);
for i = 1:length(Betaxi)
B1 = Betaxi(i);
J1 = besselj(mu,B1).^2;
R = (v + (mu.*thisN));
S = ((v + (mu.*thisN)).^2);
Z = (1+u.*((I2.*exp(1i.*mu.*p))+1)).*trig;
Z1 = (1+u.*((I5.*exp(1i.*mu.*p))+1));
tmp = ((J1.*R)./(1+(2.*u)+(u.^2)+S)).*Z;
tmp1 = ((J1.*R)./(1+(2.*u)+(u.^2)+S)).*Z1;
J(i) = joxi.*sum(tmp(:));
J1(i) = joxi.*sum(tmp1(:));
JJ = J(i) +J1(i);
end
legendStrings{k1} = sprintf('w = %.3f', thisN);
plot(Betaxi, real(JJ), '-', 'LineWidth', 2, 'MarkerSize', 15);
hold on;
drawnow;
end
grid on;
fontSize = 20;
xlabel('\beta_\xi', 'FontSize', fontSize)
ylabel('\it j_\xi /Am^{-2}', 'FontSize', fontSize)
title('\it j_\xi /Am^{-2} vs. \beta_\xi [Axial axis]', 'FontSize', fontSize)
legend(legendStrings, 'Location', 'best');
hFig.WindowState = 'maximized';
