design={"A","B","C"};
#design={"A"};

#d2=0.4e-6;
d2=linspace(0.27e-6,0.75e-6,20); # uncomment this line for a matrix of thickness
neff1 = neff2 = ng1 = ng2 = matrix(length(d2),length(design));
NRPS = 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));
        
        # forward simulations
        s=1;
        MO_NRPS_Bahlmann_material;
        
        setanalysis("use max index",true);
        
        n = findmodes;
        ?mname = findfirstTMmode(n);
        neff1(i,j) = getdata(mname,"neff");
        ng1(i,j) = getdata(mname,"ng");
        
        # backward simulations
        switchtolayout;
        s=-1;
        MO_NRPS_Bahlmann_material;
        
        n = findmodes;
        ?mname = findfirstTMmode(n);
        neff2(i,j) = getdata(mname,"neff");
        ng2(i,j) = getdata(mname,"ng");
        ?NRPS(i,j) = real(2*pi/1.3e-6*(neff1(i,j)-neff2(i,j)));
        
    }
}
    
plot(d2*1e6,NRPS*1e-2,"total film thickness d2 (microns)","NRPS (1/cm)","Figure 2 Bahlmann et al.","linewidth=3");
legend(design);