##############################################################
# Scriptfile: binary_phase_mask_design.lsf
#
# Description: This file is used to plot the unpolarized
#    0th order transmission calculated by the parameter
#    sweep object in the fsp file binar_phase_mask_duty_cycle_sweep.fsp.
#    It then calculates the best duty cycle profile to target a 
#    parabolic intensity pattern in the 0th order. It saves
#    the target duty cycle as a function of device period to
#    the matlab file target_profile.mat.
#
# Copyright 2014 Lumerical Solutions Inc.
##############################################################

load("binary_phase_mask_duty_cycle_sweep");

# run the sweep if it does not already have the data
if(!havesweepresult("polarization")) { 
  runsweep("polarization");
  save;
}

# collect the 0th order transmisison
T0 = getsweepresult("polarization","T0");
T0_0 = pinch(T0.T0,2,1);
T0_90 = pinch(T0.T0,2,2);
T0_unpolarized = (T0_0+T0_90)/2;
plot(T0.duty_cycle,T0_0,T0_90,T0_unpolarized,"duty cycle","T0");
legend("polarization 0","polarization 90","unpolarized");

# chose the target profile
a = 3e-6;
x_mask = (-15:15)*a; #-45 to 45 micrions
Itarget = x_mask^2/max(x_mask)^2;
plot(x_mask*1e6,Itarget,"x (microns)","I","target intensity profile");
setplot("show legend",false);

#normalize the calculated I0 profile to a max value of 1
I0 = T0_unpolarized/max(T0_unpolarized);
duty_cycle0 = T0.duty_cycle;
#choose upper half or lower half
pos_upper = find(I0,min(I0)):length(I0);
pos_lower = 1:find(I0,min(I0));
# keep only the lower half
I0 = I0(pos_lower);
duty_cycle0 = duty_cycle0(pos_lower);
# find the duty cycle that gives the target intensity
duty_cycle = interp(duty_cycle0,I0,Itarget);
plot(x_mask*1e6,duty_cycle,"x (microns)","duty cycle","duty cycle vs position");
setplot("show legend",false);

# save the target profile to a file
matlabsave("target_profile",duty_cycle,Itarget,x_mask,a);