###############################################
# script file: nanobeams_grating.lsf
#
#     This script computes the diffraction efficiency and deflection angles for the nanobeam grating example.
#     The script sweeps the angle of incidence of the plane wave source and collects the relative power into
#     each propagating grating order together with the deflection angles for each source angle. The diffraction
#     efficiency is taken to be the relative power into the first grating order.
#
# Copyright 2017 Lumerical Solutions
###############################################

clear;
closeall;

# -------
# Inputs:
# -------

incidenceAngle = linspace(-20.0,20.0,21);
fileName = "nanobeam_grating";
load(fileName);

# ---------------
# Grating Orders:
# ---------------

# retrieve substrate refractive index
centerFreq = getnamed("DGTD::source","center frequency");
backIndexT = getdgtdindex("substrate",centerFreq,centerFreq,centerFreq);
# retrieve periodicity vector (meters)
periodicityVector = getperiodicity("DGTD");
# compute k-vector corresponding to each source angle
sourceKUnitVector = matrix(3,length(incidenceAngle));
switchtolayout;
for(jjj = [1:length(incidenceAngle)]){
   setnamed("DGTD::source","angle theta",incidenceAngle(jjj));
   sourceKUnitVector(:,jjj) = getsourcedirection("DGTD::source");
}
clear(jjj);
# compute the number of T/R grating orders
goT = gratingorders(periodicityVector,sourceKUnitVector(:,1),centerFreq,backIndexT);
goR = gratingorders(periodicityVector,sourceKUnitVector(:,1),centerFreq);

# -------------------
# Grating Projection:
# -------------------

normProjectionT = deflectionAngleT = matrix(length(incidenceAngle),length(goT.N));
normProjectionR = deflectionAngleR = matrix(length(incidenceAngle),length(goR.N));
for(jjj = [1:length(incidenceAngle)])
{ 
    # simulation
    switchtolayout;
    setnamed("DGTD::source","angle theta",incidenceAngle(jjj));
    run;
    # T projection
    fieldsT = getresult("DGTD::T","fields");
    fluxT = getresult("DGTD::T","flux");
    gp = gratingprojection(fieldsT,periodicityVector,sourceKUnitVector(:,jjj),backIndexT);
    normProjectionT(jjj,:) = gp.projection*fluxT.T_front;
    deflectionAngleT(jjj,:) = gp.theta;
    # R projection
    fieldsR = getresult("DGTD::R","fields");
    fluxR = getresult("DGTD::R","flux");
    gp = gratingprojection(fieldsR,periodicityVector,sourceKUnitVector(:,jjj));
    normProjectionR(jjj,:) = gp.projection*fluxR.T_back;
    deflectionAngleR(jjj,:) = gp.theta;
}
clear(jjj,fieldsT,fieldsR,fluxT,fluxR,gp);

# ------
# Plots:
# ------

# relative power for all transmitted and reflected grating orders
plot(-incidenceAngle,normProjectionT,"incidence angle (deg)","relative power");
legStr1 = cell(length(goT.N));
for(iii = 1:length(goT.N)){
    legStr1{iii} = "T - N = "+num2str(goT.N(iii));
}
holdon;
plot(-incidenceAngle,normProjectionR,"incidence angle (deg)","relative power");
legStr2 = cell(length(goR.N));
for(iii = 1:length(goR.N)){
    legStr2{iii} = "R - N = "+num2str(goR.N(iii));
}
holdoff;
legend(legStr1,legStr2);

# plot deflection angle for all transmitted orders
plot(-incidenceAngle,deflectionAngleT,"incidence angle (deg)","deflection angle (deg)");
holdon;
plot(-incidenceAngle,deflectionAngleR,"incidence angle (deg)","deflection angle (deg)");
holdoff;
legend(legStr1,legStr2);
clear(iii,legStr1,legStr2);

# load the reference reflection angle
deflectionAngleRefData = readdata("nanobeam_grating_ref_angle.txt");
incidenceAngleRef = deflectionAngleRefData(:,1);
deflectionAngleRef = deflectionAngleRefData(:,2);
deflectionAngleRef = interp(deflectionAngleRef,incidenceAngleRef,-incidenceAngle);
clear(deflectionAngleRefData,incidenceAngleRef);

# plot the deflection angle for the first order aganinst the reference result
plot(-incidenceAngle,deflectionAngleRef,deflectionAngleT(:,3),"incidence angle (deg)","deflection angle (deg)");
legend("reference","DGTD");

# load the reference diffraction efficiency
efficiencyRefData = readdata("nanobeam_grating_ref_efficiency.txt");
incidenceAngleRef = efficiencyRefData(:,1);
efficiencyRef = efficiencyRefData(:,2);
efficiencyRef = interp(efficiencyRef,incidenceAngleRef,-incidenceAngle);
clear(efficiencyRefData,incidenceAngleRef);

# plot the diffraction efficiency for the first order against the reference result
plot(-incidenceAngle,efficiencyRef*100.0,normProjectionT(:,3)*100.0,"incidence angle (deg)","diffraction efficiency (%)");
legend("reference","DGTD");