#arrow waveguide
#this script calculates the error for modes of an arrow  waveguide as a function of grid accuracy

use_matlab=0;
conformal_mesh_on=1;

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

#The analytic answer
lam=632.8e-9;
N=1.457941265;
L=5.4189e-8*40*pi/lam/log(10)*1e-2; #(in dB/cm)

mesh_cell = matrix(1,7);
mesh_cell2 = matrix(1,7);
mesh_cell(1) = 200;
for(i=2:length(mesh_cell)){
   mesh_cell(i)=mesh_cell(i-1)*2;
}
N2=matrix(1,length(mesh_cell));
L2=matrix(1,length(mesh_cell));

for(ii = 1:length(mesh_cell)){
    switchtolayout;
    if(conformal_mesh_on){setnamed("MODE","mesh refinement","conformal variant 1");}
    setnamed("MODE","mesh cells y",mesh_cell(ii));
    findmodes;
    N2(ii) = real(getdata(bestoverlap("TE1"),"neff"));
    L2(ii) = getdata(bestoverlap("TE1"),"loss")/100;
    simulation;
    mesh_cell2(ii)=getnamed("MODE","actual mesh cells y");
}   

p3 = 100*abs((N2-N)/N);
p4 = 100*abs((L2-L)/L);

if(p3(4)>1e-5){
   ?"ERROR: mode index exceeded error tolerance in arrow_waveguide";
   #break;
}
if(p4(4)>1){
   ?"ERROR: mode loss exceeded error tolerance in arrow_waveguide";
   #break;
}

N=matrix(1,length(mesh_cell2))+N;
L=matrix(1,length(mesh_cell2))+L;

if(use_matlab){
    matlabput(mesh_cell2,p3,p4,N,N2,L,L2);
    matlab("
       close all;
       figure(1);
       subplot('position',[0.2143 0.1738 0.6 0.7512]);
       [ax,h1,h2]=plotyy(mesh_cell2,N2,mesh_cell2,L2,@semilogx);
       set(h1,'LineWidth',4,'LineStyle','none','Marker','o','Color','r','MarkerSize',12);
       set(h2,'LineWidth',4,'LineStyle','none','Marker','o','Color','b','MarkerSize',12);
       set(ax(1),'YColor','r');
       set(ax(2),'YColor','b');
       set(ax(1),'FontSize',20,'LineWidth',2);
       set(ax(2),'FontSize',20,'LineWidth',2);
       axes(ax(1));
       line([1e2 1e5],[N(1) N(1)],'LineWidth',4,'Color','r');
       ylabel('effective index','FontSize',20);
       xlabel('grid points','FontSize',20);
       set(ax(1),'XGrid','on');
       axes(ax(2));
       line([1e2 1e5],[L(1) L(1)],'LineWidth',4,'Color','b');
       ylabel('propagation loss (dB/cm)','FontSize',20);
       box on;
       print -djpeg90 arrow_waveguide.jpg;
       figure(2);
       loglog(mesh_cell2,p3,'-ro',mesh_cell2,p4,'-bo','LineWidth',4,'MarkerSize',12);
       set(gca,'FontSize',20,'LineWidth',2);
       h=legend('effective index','propagation loss','Location','NorthEast');
       set(h,'LineWidth',2);
       axis([1e2 1e5 1e-8 1e1]);
       h=xlabel('grid points','FontSize',20);
       h=ylabel('error amplitude (%)','FontSize',20);
       grid on;
       box on;
       print -djpeg90 arrow_waveguide2.jpg;
    ");
}else{
    plot(log10(mesh_cell2),N2,N,"grid points","effective index","arrow_waveguide","");
    legend("eigenmode solver","comparison results");
    setplot("x min", 2);
    setplot("x max", 5);
    #setplot("y min", 1.45794);
    #setplot("y max", 1.45794);
    exportfigure("lum_arrow_waveguide");
    plot(log10(mesh_cell2),L2,L,"grid points","propagation loss (dB/cm)","arrow_waveguide_1","");
    legend("eigenmode solver","comparison results");
    setplot("x min", 2);
    setplot("x max", 5);
    setplot("y min", 0.04);
    setplot("y max", 0.06);
    exportfigure("lum_arrow_waveguide_1");
    plot(log10(mesh_cell2),p3,p4,"grid points","error amplitude (%)","arrow_waveguide_2","logplot");
    legend("effective index","propagation loss");
    setplot("x min", 2);
    setplot("x max", 5);
    setplot("y min", -8);
    setplot("y max", 1);
    exportfigure("lum_arrow_waveguide_2");
}