% This prog plots k_r plot of WKB approximation 
function ECE_635_WKB
close all;clear;clc;clf reset
global n1 q delta k0 L a beta
n1 = 1.5; n2 = 1.485;a = 12e-6;q = 2;rn = 0:0.01:1;rnc = 1.01:0.01:2;  
delta = (n1^2 - n2^2)/(2*n1^2);%n1 = 1.50;n2 = 1.47;
nr = n1*sqrt(1 - 2*rn.^q*delta);L = 0;m = 2;%%%%% Graded index
%nr = n1;L = 0;m = 1; %%%%% Step index 
lamda = 1.55e-6;k0 = 2*pi/lamda;k1 = 2*pi*n1/lamda;k2 = 2*pi*n2/lamda;
beta = sqrt(a*k0^2*n1^2 - 4*k0*n1*sqrt(2*delta)*(m + L/2))/sqrt(a)
[k1 k2]
%beta = 0.95*k0*n2; %%%% Selecting beta for refracting ray
kr2cor = nr.^2*k0^2 - L^2./(a^2*rn.^2) - beta^2;
kr2clad = n2.^2*k0^2 - L^2./(a^2*rnc.^2) - beta^2;
kr2 = [kr2cor kr2clad];rnp = [rn rnc]; 
plot(rnp, kr2, '-k  ', 'LineWidth',3);hold on 
plot(rnp, zeros(1,length(rnp)), '--k ', 'LineWidth',2);
hold off; 
set(gcf,'Color',[1 1 1]);set(gca,'FontSize',18)
%title('Plot of \itk_r over core and cladding regions','FontSize',12,'FontWeight','bold');
xlabel('\itr\rm\bf_n = \itr / a','FontSize',22,'FontWeight','bold');
ylabel('\itk_r\rm\bf^2 ( \itr\rm\bf )','FontSize',22,'FontWeight','bold');
axis ([min(rnp)*1 max(rnp)*1 -1e12 max(kr2)*1.2]);
%%%% Test %%%%
%[a*k0*n1*sqrt(delta)/(2*sqrt(2)) m + L/2]
%%%%% Finding turning points
%root_r_L = fzero(@funkr,[0.01 0.5])
root_R_L = fzero(@funkr,[0.5 1])
function fB = funkr(rx)
global n1 q delta k0 L a beta
fB = (n1*sqrt(1 - 2*rx.^q*delta)).^2*k0^2 - L^2./(a^2*rx.^2) - beta^2;