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

d2=0.4e-6;
#d2=linspace(0.27e-6,0.5e-6,10); # uncomment this line for a matrix of thickness
NRPS=matrix(length(d2),length(design));
matlabload("data_fde"); # loading the group index data from FDE simulation
if (length(ng(:,1))!=length(d2)){
    ?"Warning! Check FDE simulation parameters"; break;
}
for(j=1:length(design)) {
    switchtolayout;
    setnamed("::model","design", design{j});
    
    for (i=1:length(d2)){
        
        switchtolayout;
        setnamed("::model","total thickness",d2(i));
        k0 = 2*pi/1.3e-6;
        att1 = getnamed("F2_rib","grid attribute name");
        att2 = getnamed("F2_film","grid attribute name");
        att3 = getnamed("F1_film","grid attribute name");
        setnamed("F2_rib","grid attribute name","");
        setnamed("F2_film","grid attribute name","");
        setnamed("F1_film","grid attribute name","");
        select("source");
        updatesourcemode;
        neff=getresult("source","neff");
        beta=neff.neff*k0;
        setnamed("FDTD","kz",real(beta));
        setnamed("F2_rib","grid attribute name",att1);
        setnamed("F2_film","grid attribute name",att2);
        setnamed("F1_film","grid attribute name",att3);
        
        # forward simulations
        s=1;
        MO_NRPS_Bahlmann_material;
        save(pwd+"\job"+num2str(i)+","+num2str(j)+"_for");
        addjob(pwd+"\job"+num2str(i)+","+num2str(j)+"_for");
        
        # backward simulations
        switchtolayout;
        s=-1;
        MO_NRPS_Bahlmann_material;
        save(pwd+"\job"+num2str(i)+","+num2str(j)+"_back");
        addjob(pwd+"\job"+num2str(i)+","+num2str(j)+"_back");
        
    } runjobs;
    
    for (i=1:length(d2)){
        
        load(pwd+"\job"+num2str(i)+","+num2str(j)+"_for");
        spectrum_for=getresult("bandstructure","spectrum");
        load(pwd+"\job"+num2str(i)+","+num2str(j)+"_back");
        spectrum_back=getresult("bandstructure","spectrum");
        
        f_for=spectrum_for.f;
        f_back=spectrum_back.f;
        
        spectrum_for=spectrum_for.fs;
        spectrum_back=spectrum_back.fs;
        
        f_for_peak=f_for(findpeaks(spectrum_for));
        f_back_peak=f_back(findpeaks(spectrum_back));
        
        # negative sign is intentional and correct
        # it is due to the translation of a frequency difference at fixed beta
        # to a difference in beta at a fixed frequency
        # The waveguide with the higher frequency resonance for a given beta
        # actually has smaller beta at fixed frequency.
        ?NRPS(i,j) = -2*pi*(f_for_peak-f_back_peak)*ng(i)/c;
        
    }
}

plot(d2*1e6,NRPS*1e-2,"total film thickness d2 (microns)","NRPS (1/cm)","Figure 2 Bahlmann et al.","linewidth=3");
legend(design);

load("MO_NRPS_Bahlmann.fsp");