# get the Ex field above the structure
t = getdata("above","t");
E = pinch(getdata("above","Ex"));

# fourier transform
Ef = fft(E);
w = fftw(t);

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

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

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

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

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

# get the level populations
N0 = getdata("medium","storage_Ex_4");
N1 = getdata("medium","storage_Ex_5");
N2 = getdata("medium","storage_Ex_6");
N3 = getdata("medium","storage_Ex_7");

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