#step index fiber
#this script calculates the error for 3 TM modes of a step index fiber as a function of grid accuracy

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("step_index_fiber.lms");
cleardcard;
loaddata("step_index_fiber.ldf");

if(use_matlab) {
    #get the analytic answer from MATLAB script
    matlab("
       clear all;
       Nth=step_index_fiber;
    ");
    matlabget(Nth);
    rm("step_index_fiber.txt");
    format long;
    write("step_index_fiber.txt", num2str(Nth));
    format short;
}else{
    Nth=readdata("step_index_fiber.txt");
}

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) = real(getdata(bestoverlap("TM01"),"neff"));
    N2(ii) = real(getdata(bestoverlap("TM02"),"neff"));
    N3(ii) = real(getdata(bestoverlap("TM03"),"neff"));

}    

p1 = 100*abs((N1-Nth(1))/Nth(1));
p2 = 100*abs((N2-Nth(2))/Nth(2));
p3 = 100*abs((N3-Nth(3))/Nth(3));

#error specifications for 160x160 grid points
if(p1(7)>6e-5){
   ?"ERROR: TM01 mode exceeded error tolerance in step_index_fiber";
   #break;
}
if(p2(7)>1e-3){
   ?"ERROR: TM02 mode exceeded error tolerance in step_index_fiber";
   #break;
}
if(p3(7)>5e-3){
   ?"ERROR: TM03 mode exceeded error tolerance in step_index_fiber";
   #break;
}

if(use_matlab){
    matlabput(mesh_cell,p1,p2,p3);
    matlab("
       close all;
       figure(1);
       loglog(mesh_cell,p1,'-ro',mesh_cell,p2,'-bo',mesh_cell,p3,'-go','LineWidth',4,'MarkerSize',12);
       set(gca,'FontSize',20,'LineWidth',2);
       h=legend('TM_{01}','TM_{02}','TM_{03}','Location','SouthWest');
       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 step_index_fiber.jpg;
    ");
}else{
    plot(mesh_cell,p1,p2,p3,"points per side","error amplitude (%)","step_index_fiber");
    setplot("log10y",true);
    setplot("log10x",true);
    legend("TM01","TM02","TM03");
    exportfigure("lum_step_index_fiber");
}