# In this example we demonstrate how to set chirp coefficient or chirp parameter to achive desired change in bragg wavelength.
# This example is instructive to help understand how the chirp is defined.
#
# Steps to run:
# 1. Open Chirped_TWLM_test.icp project file.
# 2. Run this script file to set chirp.
# 3. Run the simulation.
# 4. Visualize ONA transmission to confirm bragg wavelength change.
#
# To set TWLM for passive mode do the following (already set):
# 1. DC source amplitude = 0
# 2. gain coefficient = 0
# 3. spontaneous emission factor = 0
# 4. linewidth enhancement factor = 0 (to prevent index variation as a function of carrier density)
# 5. Add ONA
# 6. Set the same frequency (units and value) in ONA as in TWLM. This is important since small discrepancy may 
#    fail the simulation since both of these elements are sources, so there may be a clash. If the simulation 
#    complains about the simulation bandwidth do explicit copy/paste from TWLM to ONA and back (make units equal prior to this).
# 7. set frequency range in ONA to %sample rate% (in expression column) to make it the same as for TWLM.
# 8. Set ONA to impulse mode (since TWLM is a time-domain-only model)
# 9. Set other desired properties of ONA (e.g. number of points to 10000)
#
#
# NOTES:
# - beta stands for propagation constant (chirp has the same units and it is additive to beta):
#   beta = beta0 + dbeta/df*delta_f = 2*pi*f0/c*neff + 1/vg*2*pi*(f-f0)
#   if neff = ng beta is just the zeroth order term 2*pi*f/c*neff
# - Grating model equation (s-matrix) is parametrized at center frequency
# - The grating transmission shows a dip at the Bragg wavelength due to quarter wave shift
#   (pi/2 phase slip in the middle)



#USER INPUT
chirp_function = "chirp coefficient";
#chirp_function = "chirp parameter";
desired_bragg_shift = -5e-9; #m (1320 nm -> 1315 nm)



neff = getnamed("DFB Laser","effective index 1");
ng = getnamed("DFB Laser","group index 1");
grating_period = getnamed("DFB Laser","grating period");
sample_rate = getnamed("::Root Element","sample rate");

lambda_bragg = 2*neff*grating_period;
vg = c/ng; #group velocity

#First subtract old grating propagation constant and add the desired propagation constant 
#(old propagation constant perturbed by chirp)
beta_bragg_unperturbed = pi/grating_period;
grating_period_shifted = (lambda_bragg + desired_bragg_shift)/2/neff;
beta_bragg_perturbed = pi/grating_period_shifted;
chirp = beta_bragg_perturbed - beta_bragg_unperturbed; #the signs are such because beta_new = beta - (beta_bragg + chirp)

#Then account for beta dispersion at the new bragg frequency if neff is different than ng
fc = getnamed("DFB Laser","frequency"); #grating model equation is parametrized at center frequency
fbragg = c/lambda_bragg; #old bragg frequency
beta_unperturbed = 2*pi*fbragg/c*neff + 1/vg*2*pi*(fc - fbragg); #old propagation constant calculated at fc
fbragg_perturbed = c/(lambda_bragg + desired_bragg_shift); #desired bragg frequency
beta_perturbed = 2*pi*fbragg_perturbed/c*neff + 1/vg*2*pi*(fc - fbragg_perturbed); #new propagation constant calculated at fc
chirp = chirp + beta_unperturbed - beta_perturbed; #the signs are such because beta_new = beta - (beta_bragg + chirp)

#Coordinates for user defined chirp function
L = getnamed("DFB Laser","length");
z = linspace(0.001,1,5000)*L; #avoid singularity at zero due to 1/z in formulas below

#Calculate and set either chirp coefficent or chirp parameter depending on user choice
setnamed("DFB Laser","chirp function","user defined");
setnamed("DFB Laser","load chirp from file",false);
setnamed("DFB Laser","user defined chirp function",chirp_function);
if(chirp_function == "chirp coefficient"){
    chirp_coeff = -chirp*lambda_bragg^2/4/pi/neff/z;
    plot(z,chirp_coeff);
    setnamed("DFB Laser","chirp table",[z/L,chirp_coeff]);
}else if(chirp_function == "chirp parameter"){
    chirp_param = chirp*L^2/z;
    plot(z,chirp_param);
    setnamed("DFB Laser","chirp table",[z/L,chirp_param]);
}else{
    assert("chirp function not supported",false);
}