# This script calculates the NRPS for the design described in Zhou et al

t=0.15e-6; # comment out this line to simulate multiple thickness
#t=linspace(0.05e-6, 0.5e-6, 10); # Si thickness for MO(+)/Si/MO(-)
#t=linspace(0.18e-6, 0.35e-6, 10); # Si thickness for MO/Si/air

NRPS=matrix(length(t));

for (i=1:length(t)) {
    
    # forward
    switchtolayout;
    setnamed("Si", "y span", t(i));
    setnamed("source","x",-0.1e-6);
    setnamed("source","direction","forward");
    run;
    Hzf = getdata("line","Hz");
    x = getdata("line","x");
    
    # backward
    switchtolayout;
    setnamed("source","x",9.5e-6);
    setnamed("source","direction","backward");
    run;
    Hzb = getdata("line","Hz");
    
    phase_diff = -(unwrap(angle(Hzf))+unwrap(angle(Hzb))); # negative sign to make phase difference positive
    p1 = find(x,1e-6); p2 = find(x,9e-6);
    p = p1:p2;
    fit = polyfit(x(p),phase_diff(p),1);
    ?NRPS(i)=fit(2);
    if(i==1) {
        plot(x(p)*1e6,phase_diff(p),fit(2)*x(p)+fit(1),"x (microns)","phase difference (rad)"); # enable to plot the phase_diff
        legend("phase difference","fit");
    }
}

plot(t*1e6,NRPS*1e-3,"thickness (microns)","NRPS (rad/mm)");