#graded index fiber
#this script calculates the error for 2 modes of a graded index fiber as a function of grid accuracy
#The script also calculates the dispersion error as a function of wavelength

use_matlab=0; # 1 to launch Matlab calculations and 0 to run MODE calculations only.
conformal_mesh_on=1; # 1 for "conformal variant 1" and 0 for "staircase"; conformal mesh provides better accuracy

load("graded_index_fiber.lms");
cleardcard;
loaddata("graded_index_fiber.ldf");

lam = 1.0e-6;
n0 = sqrt(4);
k0 = 2*pi*n0/lam;
H = 40e-6;
m = 0;
n = 0;
alpha = 2+2*(m+n);
beta = k0*sqrt(1-alpha/k0/H);
None = n0*sqrt(1-alpha/k0/H);

m = 0;
n = 1;
alpha = 2+2*(m+n);
Ntwo = n0*sqrt(1-alpha/k0/H);

m = 1;
n = 1;
alpha = 2+2*(m+n);
Nthree = n0*sqrt(1-alpha/k0/H);

mesh_cell = matrix(1,10);
mc=20;
for(i=1:length(mesh_cell)){
   mesh_cell(i)=round(mc);
   mc=mc*sqrt(2);
}
N1=matrix(1,length(mesh_cell));
N2=matrix(1,length(mesh_cell));
N3=matrix(1,length(mesh_cell));

for(ii = 1:length(mesh_cell)){
    switchtolayout;
    if(conformal_mesh_on){
      setnamed("MODE","mesh refinement","conformal variant 1");
    } else {
      setnamed("MODE","mesh refinement","staircase");    
    }
    setnamed("MODE","mesh cells x",mesh_cell(ii));
    setnamed("MODE","mesh cells y",mesh_cell(ii));
    findmodes;
    N1(ii) = getdata(bestoverlap("TE00"),"neff");
    N2(ii) = getdata(bestoverlap("TE01"),"neff");
    N3(ii) = getdata(bestoverlap("TE11"),"neff");

    #do dispersion calculation only for 80 grid cells
    if(mesh_cell(ii)==80){
    #if(mesh_cell(ii)==453){
        selectmode(3);  # TE01 mode
        setanalysis("detailed dispersion calculation", 0); #disable detailed dispersion
        frequencysweep;
        neff = getdata("frequencysweep", "neff");
        D2=getdata("frequencysweep","D");
        f_neff=getdata("frequencysweep", "f"); 
        f_D=getdata("frequencysweep","f_D");

        w_D=2*pi*f_D;
        m = 0;
        n = 1;
        alpha = 2+2*(m+n);
        D2th=alpha^2/(8*pi*H^2*n0)*(w_D^2)*((w_D^2-alpha*c*w_D/H/n0)^(-3/2));        
    }
}

p1 = 100*abs((N1-None)/None);
p2 = 100*abs((N2-Ntwo)/Ntwo);
p3 = 100*abs((N3-Nthree)/Nthree);
D2err=abs(D2th-D2)*100/D2th;
lam2=c/f_D*1e6;

#error specifications for 160x160 grid points
if(p2(7)>2e-4){
   ?"ERROR: TE01 mode exceeded error tolerance in graded_index_fiber";
   #break;
}
if(p3(7)>1e-3){
   ?"ERROR: TE11 mode exceeded error tolerance in graded_index_fiber";
   #break;
}
if(max(D2err)>1){
   ?"ERROR: TE01 dispersion exceeded error tolerance in graded_index_fiber";
   #break;
}

if(use_matlab){
    matlabput(mesh_cell,p2,p3,D2err,lam2,D2,D2th);
    matlab("
       close all;
       figure(1);
       loglog(mesh_cell,p2,'-ro',mesh_cell,p3,'-bo','LineWidth',4,'MarkerSize',12);
       set(gca,'FontSize',20,'LineWidth',2);
       h=legend('TE_{01}','TE_{11}');
       set(h,'LineWidth',2);
       axis([10 1000 1e-6 1]);
       h=xlabel('points per side','FontSize',20);
       h=ylabel('error amplitude (%)','FontSize',20);
       grid on;
       box on;
       print -djpeg90 graded_index_fiber.jpg;
       figure(2);
       plot(lam2,D2th*1e6,'-r',lam2,D2*1e6,'ro','LineWidth',4,'MarkerSize',12);
       set(gca,'FontSize',20,'LineWidth',2);
       h=legend('analytic value','calculated value');
       set(h,'LineWidth',2);
       axis([1 1.5 0.1 0.2]);
       h=xlabel('wavelength (microns)','FontSize',20);
       h=ylabel('dispersion (ps/nm/km)','FontSize',20);
       grid on;
       box on;
       print -djpeg90 graded_index_fiber2.jpg;
       figure(3);
       semilogy(lam2,D2err,'-ro','LineWidth',4,'MarkerSize',12);
       set(gca,'FontSize',20,'LineWidth',2);
       h=legend('80 points per side','Location','SouthWest');
       set(h,'LineWidth',2);
       axis([1 1.5 1e-3 1e0]);
       h=xlabel('wavelength (microns)','FontSize',20);
       h=ylabel('dispersion error (%)','FontSize',20);
       grid on;
       box on;
       print -djpeg90 graded_index_fiber3.jpg;
    ");
}else{
    plot(mesh_cell,p2,p3,"points per side","error amplitude (%)","graded_index_fiber");
    legend("TE_01","TE_11");
    setplot("log10y",true);
    setplot("log10x",true);
    exportfigure("lum_graded_index_fiber");
    
    plot(lam2,D2th*1e6,D2*1e6,"wavelength (microns)","dispersion (ps/nm/km)","graded_index_fiber_2");
    legend('analytic value','calculated value');
    exportfigure("lum_graded_index_fiber_2");
    
    plot(lam2,D2err,"wavelength (microns)","dispersion error (%)","graded_index_fiber_3");
    legend("80 points per side");
    setplot("log10y",true);
    exportfigure("lum_graded_index_fiber_3");  
}

plot((c/f_neff)/1e-6, real(neff), "wavelength (microns)","effective index","neff for TE01 from numerical calculation");
legend("n80 points perside");