#hollow metal wg
#This script calculates the dispersive properties of a hollow metal waveguide
#This script also examines the error in the calculated wavenumber as a function
#of number of grid points.

#if not set flags
use_matlab=0;
conformal_mesh_on=1;

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

#stage 1: look at wavenumber for 3 modes from 5GHz to 25GHz

#calculate anlaytic freq response
a = 2e4 * 1e-6; #length of wg
b = 1e4 * 1e-6; #width of wg
fth = linspace(5e9,25e9,1401);
wth = 2*pi*fth;
k=matrix(3,1401);

m=1; n=0;
k(1,1:1401) = sqrt( (wth/c)^2 - pi^2*( (m/a)^2 + (n/b)^2 ) );
m=0; n=1;
k(2,1:1401) = sqrt( (wth/c)^2 - pi^2*( (m/a)^2 + (n/b)^2 ) );
m=1; n=1;
k(3,1:1401) = sqrt( (wth/c)^2 - pi^2*( (m/a)^2 + (n/b)^2 ) );
k=real(k);

switchtolayout;
setanalysis("number of trial modes",5);
if(conformal_mesh_on){setnamed("MODE","mesh refinement","conformal variant 1");}
setnamed("MODE","mesh cells x",40);
setnamed("MODE","mesh cells y",80);
findmodes;
selectmode(bestoverlap("TE10"));
frequencysweep;
f=getdata("frequencysweep","f");
TE10=2*pi*f/c*real(getdata("frequencysweep","neff"));
selectmode(bestoverlap("TE01"));
frequencysweep;
TE01=2*pi*f/c*real(getdata("frequencysweep","neff"));
selectmode(bestoverlap("TE11"));
frequencysweep;
TE11=2*pi*f/c*real(getdata("frequencysweep","neff"));

f=f/1e9;
fth=fth/1e9;

if(use_matlab){
    matlabput(fth,k,f,TE10,TE01,TE11);
    matlab("
       close all;
       f10=f(TE10>1e-6);
       f01=f(TE01>1e-6);
       f11=f(TE11>1e-6);
       TE10=TE10(TE10>1e-6);
       TE01=TE01(TE01>1e-6);
       TE11=TE11(TE11>1e-6);
       TE10t=k(1,:);
       TE01t=k(2,:);
       TE11t=k(3,:);
       f10t=fth(TE10t>1e-6);
       f01t=fth(TE01t>1e-6);
       f11t=fth(TE11t>1e-6);
       TE10t=TE10t(TE10t>1e-6);
       TE01t=TE01t(TE01t>1e-6);
       TE11t=TE11t(TE11t>1e-6);
       plot(f10t,TE10t,'r-',f01t,TE01t,'b-',f11t,TE11t,'g-','LineWidth',4,'MarkerSize',12);
       hold on;
       plot(f10,TE10,'ro',f01,TE01,'bo',f11,TE11,'go','LineWidth',4,'MarkerSize',12);
       set(gca,'FontSize',20,'LineWidth',2);
       h=legend('TE_{01}','TE_{01}','TE_{11}','Location','NorthWest');
       set(h,'LineWidth',2);
       axis([5 25 0 500]);
       h=xlabel('frequency (GHz)','FontSize',20);
       h=ylabel('wavenumber (m^{-1})','FontSize',20);
       grid on;
       box on;
       print -djpeg90 hollow_metal_wg.jpg;
    ");
}else{
    TE10t=interp(k(1,1:length(fth)),fth,f);
    TE01t=interp(k(2,1:length(fth)),fth,f);
    TE11t=interp(k(3,1:length(fth)),fth,f);

    plot(f,TE10,TE10t,TE01,TE01t,TE11,TE11t,"frequency (GHz)","wavenumber (1/m)","hollow_metal_wg","");
    legend("TE10","TE10t","TE01","TE01t","TE11","TE11t");
    setplot("x min", 5);
    setplot("x max", 25);
    setplot("y min", 0);
    setplot("y max", 500);
    exportfigure("lum_hollow_metal_wg");
}

#stage 2: look at error in index at 20GHz
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;
    setanalysis("number of trial modes",5);
    if(conformal_mesh_on){setnamed("MODE","mesh refinement","conformal variant 1");}
    setnamed("MODE","mesh cells x",mesh_cell(ii)/2);
    setnamed("MODE","mesh cells y",mesh_cell(ii));
    setanalysis("frequency",20e9);
    findmodes;
    N1(ii) = 2*pi*20e9/c*real(getdata(bestoverlap("TE10"),"neff"));
    N2(ii) = 2*pi*20e9/c*real(getdata(bestoverlap("TE01"),"neff"));
    N3(ii) = 2*pi*20e9/c*real(getdata(bestoverlap("TE11"),"neff"));
}

a = 2e4 * 1e-6; #length of wg
b = 1e4 * 1e-6; #width of wg
wth = 2*pi*20e9;

m=1; n=0;
N1th = sqrt( (wth/c)^2 - pi^2*( (m/a)^2 + (n/b)^2 ) );
m=0; n=1;
N2th = sqrt( (wth/c)^2 - pi^2*( (m/a)^2 + (n/b)^2 ) );
m=1; n=1;
N3th = sqrt( (wth/c)^2 - pi^2*( (m/a)^2 + (n/b)^2 ) );
p1 = 100*abs((N1-N1th)/N1th);
p2 = 100*abs((N2-N2th)/N2th);
p3 = 100*abs((N3-N3th)/N3th);

if(p1(7)>5e-4){
   ?"ERROR: TE10 mode exceeded error tolerance in hollow_metal_wg";
   #break;
}
if(p2(7)>2e-2){
   ?"ERROR:  TE01 mode exceeded error tolerance in hollow_metal_wg";
   #break;
}
if(p3(7)>2.5e-2){
   ?"ERROR:  TE11 mode exceeded error tolerance in hollow_metal_wg";
   #break;
}

if(use_matlab){
    matlabput(mesh_cell,p1,p2,p3);
    matlab("
       figure(2);
       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('TE_{10}','TE_{01}','TE_{11}','Location','SouthWest');
       set(h,'LineWidth',2);
       axis([10 1000 1e-6 1]);
       h=xlabel('points on long side','FontSize',20);
       h=ylabel('error amplitude (%)','FontSize',20);
       grid on;
       box on;
       print -djpeg90 hollow_metal_wg2.jpg;
    ");
}else{
    plot(log10(mesh_cell),p1,p2,p3,"points on long side","error amplitude(%)","hollow_metal_wg_2","logplot");
    legend("TE10","TE01","TE11");
    setplot("x min", 1);
    setplot("x max", 3);
    setplot("y min", -6);
    setplot("y max", 0);
    exportfigure("lum_hollow_metal_wg_2");
}