% Modified Matlab routine for IEEE_802.3aq examples of "kink" in DMD.
% code which generates index profile & perturbation are included
% code which calculates mode delays and estimated DMD measurement has been taken out.
%----------------------------------------------------------
% Note that this is an "example" of a profile causing a "kink" in the DMD and is
% not claimed to be a paradigm, archtype, or best example.
%--------------------------------------------------------------

%====================================================================
% Plot index profile & index error for kink.
% Share with IEEE_802.3aq  channel modeling ad hoc
%  John Abbott Corning 10/27/2004
%=============================================================


Delta =  zeros(111,2); %    % Call it delta (2 columns
r=0.:.01:1.1;
Delta2=r';
Delta(:,2)=Delta2;
Delta1=1.946*(1.-r.^1.97)/(1-r(10)^1.97);
for i=1:length(r)
   if(r(i)>1.00)
      Delta1(i)=0.;
   end
end
Delta1=Delta1';
Delta(:,1)=Delta1;
%-------------------------------------------------------
figure(1)

x1 = Delta(:,2) ;
y1 = Delta(:,1);
y0=-y1;
whos
subplot(2,2,1),...
   plot(x1,y0,'b-'),,...
xlabel(' (r/a)'),ylabel('Delta'),...
title (' Profile')
% fiber dmd axis ([0 1.2 -1. 0.2])
axis ([0 1.1 0. 2.]);
% Now generate  the index error Nr
%  Nr =  (Delta(r) - Delta_opt(r)) / Delta_0  * 100.
%
N_r = Delta(:,1);
N_0 = N_r(10,1);
N_opt = N_0.*(1.-x1.^1.97)/(1.-x1(10)^1.97);  % typical 1300nm alpha, assumes data at 1300nm

N_opt(101:111,1) = 0.;

%-----------kink----------------
a_k = .02 ; 
b_k = 20*5/31.25;
D_N = 0.*x1;
r_k = .56;
for i_rad = 1:length(x1)
   D_N(i_rad) = a_k*(exp(-b_k*abs((x1(i_rad)-r_k)))-exp(-b_k*abs((x1(10)-r_k))));
      
   
end
D_N(101:111,1)=0.;
N_r1 = N_opt + D_N;
N_r = N_r1;
%----------end kink----------------------------

 N_r = (N_r - N_opt)./N_0 .*100.;
 hold on
 plot(x1,N_r1,'b')
plot(x1,N_opt,'r-');
hold off
%-----------------------------
subplot(2,2,2),...
  
plot(x1,N_r,'r-');
xlabel(' (r/a)'),ylabel('Percent Delta Error'),...
   title (' Profile Error = N_r'),axis ([0 1.1 -3. 3. ]);
%-----------------------------------------------

r33 = 0:1/31.25:32/31.25;
N_r33 = interp1(x1,N_r,r33,'linear');
hold on
plot(r33,N_r33,'o',...
   'MarkerFaceColor','g',...
   'MarkerSize',4);
hold off
%====================================
   

figure(5)
   plot(x1,y0,'b-'),,...
xlabel(' (r/a)'),ylabel('Delta'),...
title (' Profile')
% fiber dmd axis ([0 1.2 -1. 0.2])
axis ([0 1.1 0. 2.]);
hold on
 plot(x1,N_r1,'b')
plot(x1,N_opt,'r-');
hold off