# This script file is for the FDE to calculate the group index
clear;
load("MO_NRPS_Bahlmann_group_index.lms");

switchtolayout;
#design={"A","B","C"}; # make sure the same parameter in the FDTD simulation
design = {"A"};

d2=0.4e-6; # make sure the same parameter in the FDTD simulation
#d2=linspace(0.27e-6,0.5e-6,10); # uncommment this line for a matrix of d2
ng=matrix(length(d2),length(design));

function findfirstTMmode(n) {
    for(i=1:n) {
        mname = "mode"+num2str(i);
        if(getdata(mname,"TE polarization fraction") < 0.5) {
            return mname;
        }
    }
    return "";
}

for(j=1:length(design)) {
    switchtolayout;
    setnamed("::model","design", design{j});
    for (i=1:length(d2)){
        
        switchtolayout;
        setnamed("::model","total thickness",d2(i));
        run;
        setanalysis("wavelength",1.3e-6);
        setanalysis("number of trial modes",5);
        setanalysis("calculate group index",true);
        n = findmodes;
        ?TMMode = findfirstTMmode(n);
        ?ng(i,j) = real(getdata(TMMode,"ng"));   
    }
    matlabsave("data_fde",ng);
}