# get the Ez field at the output
t = getdata("output","t");
Ez = pinch(getdata("output","Ez"));

# fourier transform
Ezf = fft(Ez);
w = fftw(t);

# apply filters in frequency domain
w1 = 2*pi*c/750e-9;
w2 = 2*pi*c/1500e-9;

filter1 = 2*exp(-(w-w1)^2/(150e12)^2);
filter2 = 2*exp(-(w-w2)^2/(150e12)^2);

# plot the spectrum and filters
plot(w/(2*pi)*1e-12,abs(Ezf)/max(abs(Ezf)),filter1,filter2,"f (THz)","spectra (a.u.)");
legend("spectrum","filter1","filter2");

# calculate the pump and lasing field in the time domain
Ezpump = invfft(Ezf*filter1);
Ezlase = invfft(Ezf*filter2);
# get new time vector, different due to zero padding
t2 = fftw(w);

# plot the pump and lasing output
plot(t2*1e12,abs(Ezpump),abs(Ezlase),"t (ps)", "|E|");
legend("pump","laser output");

# get the level populations
N0 = getdata("medium","storage_Ez_4");
N1 = getdata("medium","storage_Ez_5");
N2 = getdata("medium","storage_Ez_6");
N3 = getdata("medium","storage_Ez_7");

# plot the level populations vs time
plot(t*1e12,N0,N1,N2,N3,"t (ps)","population");
legend("N0","N1","N2","N3");