closeall;
clear;
w0 = 2*pi*c/1580e-9; # center frequency
Imin = 0.01; Imax = 0.5;

load("forward.fsp");
Ey = getdata('monitor_2','Ey');
Esy = getdata('monitor_2_1','Ey');
t = getdata('monitor_2','t');
I = linspace(Imin, Imax, length(t));
zero_pad = 2^16;
w = fftw(t,1,zero_pad);
filter1 = 2*exp(-(w-w0)^2/(150e12)^2);

# remove high frequency components
Ey_through_w = 2*((1:length(w))<=(length(w)/2+0.1))*fft(pinch(Ey),1,zero_pad);
Ey_through_t = invfft(pinch(Ey_through_w)*filter1);
Ey_through_t = Ey_through_t(1:length(t));
Esy_through_w = 2*((1:length(w))<=(length(w)/2+0.1))*fft(pinch(Esy),1,zero_pad);
Esy_through_t = invfft(pinch(Esy_through_w)*filter1);
Esy_through_t = Esy_through_t(1:length(t));

load("reverse.fsp");
Ey_r = getdata('monitor_2','Ey');
Esy_r = getdata('monitor_2_1','Ey');
t = getdata('monitor_2','t');
I_r = linspace(Imax, Imin, length(t));
zero_pad = 2^16;
w = fftw(t,1,zero_pad);
filter1 = 2*exp(-(w-w0)^2/(150e12)^2);

# remove high frequency components
Ey_through_w_r = 2*((1:length(w))<=(length(w)/2+0.1))*fft(pinch(Ey_r),1,zero_pad);
Ey_through_t_r = invfft(pinch(Ey_through_w_r)*filter1);
Ey_through_t_r = Ey_through_t_r(1:length(t));
Esy_through_w_r = 2*((1:length(w))<=(length(w)/2+0.1))*fft(pinch(Esy_r),1,zero_pad);
Esy_through_t_r = invfft(pinch(Esy_through_w_r)*filter1);
Esy_through_t_r = Esy_through_t_r(1:length(t));

# plot results
plotxy(I, abs(Ey_through_t)^2/abs(Esy_through_t)^2,
I_r, abs(Ey_through_t_r)^2/abs(Esy_through_t_r)^2,
'I (MW/cm^2)','(E_out/E_source)^2');
legend("Increasing Source Intensity", "Decreasing Source Intensity");