#########################################################################
# scriptfile: kerr1.lsf
#
# Description: This script file analyses and compares the results
#	of kerr1_linear.fsp and kerr1_nonlinear.fsp.
#
# Copyright 2007 Lumerical Solutions Inc.
########################################################################

############################################
# COLLECT NONLINEAR RESULTS
# load the linear simulation file
load("kerr1_nonlinear");

# make sure that CW normalization is used for the analysis
cwnorm;

# get the frequency if the fourier monitors
f = getdata("transmission","f");

# calculate the power vs time at the beam waist
power_vs_t_nonlinear = getdata("time_waist","power");
t_nonlinear = getdata("time_waist","t");

# calculate the average power over the simulation times
power_avg_nonlinear = integrate(power_vs_t_nonlinear,2,t_nonlinear)/(max(t_nonlinear)-min(t_nonlinear));

# calculate the average power assuming an 80 MHz pulse rate
power_avg_nonlinear_80MHz = integrate(power_vs_t_nonlinear,2,t_nonlinear)*80e6;

# display results on source power
?"note: <> represents time average over an optical cycle";
?"nonlinear:";
?"    maximum <power> measured from P(t) at waist = " + num2str(0.5*max(power_vs_t_nonlinear)) + " Watts";
?"    maximum <power> measured from P(w) at waist = " + num2str(transmission("waist")*sourcepower(f)) + " Watts";
?"    maximum <power> calculated at from source = " + num2str(sourcepower(f)) + " Watts";
?"    Average power over simulation time = " + num2str(power_avg_nonlinear ) + " Watts";
?"    Average power at 80MHz pulse rate = " + num2str(power_avg_nonlinear_80MHz) + " Watts";

# collect the data for the x-normal monitor
x_x_normal=getdata("x-normal","x");
y_x_normal=getdata("x-normal","y");
z_x_normal=getdata("x-normal","z");
E2_x_normal_nonlinear=getelectric("x-normal");

# collect data for the monitor at the beam waist
x_waist=getdata("waist","x");
y_waist=getdata("waist","y");
z_waist=getdata("waist","z");
E2_waist_nonlinear=getelectric("waist");

# collect data for the monitor at the beam waist
x_y_normal=getdata("y-normal","x");
y_y_normal=getdata("y-normal","y");
z_y_normal=getdata("y-normal","z");
E2_y_normal_nonlinear=getelectric("y-normal");

############################################
# COLLECT LINEAR RESULTS
# load the linear file
load("kerr1_linear");

# make sure that CW normalization is used for the analysis
cwnorm;

# calculate the power vs time at the beam waist
power_vs_t_linear = getdata("time_waist","power");
t_linear = getdata("time_waist","t");

# calculate the average power over the simulation times
power_avg_linear = integrate(power_vs_t_linear,2,t_linear)/(max(t_linear)-min(t_linear));

# calculate the average power assuming an 80 MHz pulse rate
power_avg_linear_80MHz = integrate(power_vs_t_linear,2,t_linear)*80e6;

# display results on source power
?"linear:";
?"    maximum <power> measured from P(t) at waist = " + num2str(0.5*max(power_vs_t_linear)) + " Watts";
?"    maximum <power> measured from P(w) at waist = " + num2str(transmission("waist")*sourcepower(f)) + " Watts";
?"    maximum <power> calculated at from source = " + num2str(sourcepower(f)) + " Watts";
?"    Average power over simulation time = " + num2str(power_avg_linear ) + " Watts";
?"    Average power at 80MHz pulse rate = " + num2str(power_avg_linear_80MHz) + " Watts";

# collect the data for the x-normal monitor
E2_x_normal_linear=getelectric("x-normal");

# collect data for the monitor at the beam waist
E2_waist_linear=getelectric("waist");

# collect data for the monitor at the beam waist
x_y_normal=getdata("y-normal","x");
y_y_normal=getdata("y-normal","y");
z_y_normal=getdata("y-normal","z");
E2_y_normal_linear=getelectric("y-normal");

########################################
# SHOW IMAGES OF THE RESULTS

# plot power(t) at the waist in both cases
plot(t_linear*1e15,power_vs_t_linear,power_vs_t_nonlinear,"time (fs)","power(t)","power(t) vs t at waist");
legend("linear","nonlinear");

# image 3 different cross sections
image(x_waist*1e9,y_waist*1e9,E2_waist_linear,"x (nm)","y (nm)","LINEAR: |E|^2 at beam waist");
image(x_waist*1e9,y_waist*1e9,E2_waist_nonlinear,"x (nm)","y (nm)","NONLINEAR: |E|^2 at beam waist");

image(y_x_normal*1e9,z_x_normal*1e9,E2_x_normal_linear,"y (nm)","z (nm)","LINEAR: |E|^2 at x=0");
image(y_x_normal*1e9,z_x_normal*1e9,E2_x_normal_nonlinear,"y (nm)","z (nm)","NONLINEAR: |E|^2 at x=0");

image(x_y_normal*1e9,z_y_normal*1e9,E2_y_normal_linear,"x (nm)","z (nm)","LINEAR: |E|^2 at y=0");
image(x_y_normal*1e9,z_y_normal*1e9,E2_y_normal_nonlinear,"x (nm)","z (nm)","NONLINEAR: |E|^2 at y=0");

# plot |E|^2 vs z at x=y=0
x = x_x_normal;
y = y_x_normal;
z = z_x_normal;
plot(z*1e9,0.5*E2_x_normal_linear(1,find(y,0),1:length(z),1),0.5*E2_x_normal_nonlinear(1,find(y,0),1:length(z),1),"z (nm)","max(<|E|^2>)","max(<|E|^2>) vs z");
legend("linear","nonlinear");

# plot the maximum refractive index change averaged over an optical cycle
chi3 = 2e-21;
eps_linear = 1.455^2 + 0*0.5*E2_x_normal_linear(1,find(y,0),1:length(z),1);
eps_nonlinear = 1.455^2 + chi3*0.5*E2_x_normal_nonlinear(1,find(y,0),1:length(z),1);
plot(z*1e9,eps_linear,eps_nonlinear,"z (nm)","max(<epsilon>)","max(<epsilon>) vs z at x=y=0");
legend("linear","nonlinear");

# attempt to find the beam radius as a function of z in x-normal plane
x = x_x_normal;
y = y_x_normal;
z = z_x_normal;
beam_radius_linear_x = matrix(length(z));
beam_radius_nonlinear_x = matrix(length(z));
y_new = linspace(min(y),max(y),20000);
for(i=1:length(z)) {
    temp = pinch(E2_x_normal_linear,3,i) / max(pinch(E2_x_normal_linear,3,i));
    temp = spline(temp,y,y_new);
    temp = temp > 0.5;
    beam_radius_linear_x(i) = 0.5*integrate(temp,1,y_new);
    temp = pinch(E2_x_normal_nonlinear,3,i) / max(pinch(E2_x_normal_nonlinear,3,i));
    temp = spline(temp,y,y_new);
    temp = temp > 0.5;
    beam_radius_nonlinear_x(i) = 0.5*integrate(temp,1,y_new);
}

# attempt to find the beam radius as a function of z in y-normal plane
x = x_y_normal;
y = y_y_normal;
z = z_y_normal;
beam_radius_linear_y = matrix(length(z));
beam_radius_nonlinear_y = matrix(length(z));
x_new = linspace(min(x),max(x),20000);
for(i=1:length(z)) {
    temp = pinch(E2_y_normal_linear,3,i) / max(pinch(E2_y_normal_linear,3,i));
    temp = spline(temp,x,x_new);
    temp = temp > 0.5;
    beam_radius_linear_y(i) = 0.5*integrate(temp,1,x_new);
    temp = pinch(E2_y_normal_nonlinear,3,i) / max(pinch(E2_y_normal_nonlinear,3,i));
    temp = spline(temp,x,x_new);
    temp = temp > 0.5;
    beam_radius_nonlinear_y(i) = 0.5*integrate(temp,1,x_new);
}

# plot the beam radius
plot(z*1e9,beam_radius_linear_x*1e9,beam_radius_nonlinear_x*1e9,
           beam_radius_linear_y*1e9,beam_radius_nonlinear_y*1e9,
           "z (nm)","beam radius (nm)","radius vs z in different planes");
legend("linear (x=0)","nonlinear (x=0)", "linear (y=0)","nonlinear (y=0)");