clear; closeall;

############################################################
# Measure impulse response of system over desired bandwidth
# at the location x,y=0,z=25um using the 'Ex' field component
# over the wavelength range 500-1500nm
# This simulation uses the standard 'short' pulse.

  switchtolayout;
  setnamed("impulse","enabled",1);
  setnamed("custom","enabled",0);
  run;

  E = getresult("impulseResponse","E");
  f = E.f;
  nf=length(f);
  Ex= pinch(E.Ex);
  impulse = pinch(E.Ex);  # use near field Ex data as impulse response


############################################################
# Next, we want to calculate the response of the system to a 
# 50fs gaussian pulse with a carrier of 1um wavelength.  The 
# maximum amplitude occurs at 200fs. The source named 'custom' is 
# already configured in this manner.
# Notice that spectrum of this pulse is completly contained
# within the range of the impulse response data collected
# in the previous section (0.5-1.5um)
# This will be done in two ways:


# Method 1: Directly measure response with an FDTD simulation
# This data will be used to confirm the results obtained 
# from the impulse response calculation
  switchtolayout;
  setnamed("impulse","enabled",0);
  setnamed("custom","enabled",1);
  run;

  # get measured data from the time monitor
  E_t_FDTD = getresult("time","E");
  response_FDTD = E_t_FDTD.Ex;
  response_FDTD_t = E_t_FDTD.t;



# Method 2: Obtain response to an arbitrary input using 
# the impulse response of the system.  

# The first step is to specify the desired time signal.
# Generally this would be directly specified by the user.
# However, in this example, we get the signal from the 'custom' 
# source so the result can be compared with the simulated result
  custom_signal_t = getdata("custom","time");  
  custom_signal   = getdata("custom","time_signal");

# Next, calculate custom signal spectrum.  Use the same frequency 
# vector as the impulse response data
dt = custom_signal_t(2)-custom_signal_t(1); 
custom_signal_f = czt(custom_signal,custom_signal_t,f*2*pi)*dt;


# calculate the time domain response of system to the custom 
# time signal by multiplying the impulse response by the 
# custom spectrum, then inverse fourier transforming.

response_impulse_f = impulse*pinch(custom_signal_f);
complex_factor = 2;                                 # frequency vector does not contain negative frequencies
dw = (f(2)-f(1))*2*pi;
response_impulse_t=linspace(0,1e-12,10000);
response_impulse = complex_factor/(2*pi)*czt(response_impulse_f,-f*2*pi,response_impulse_t)*dw;  


# plot results
plot(c/f*1e9,abs(impulse), 
             "wavelength (nm)","spectrum","Impulse response");
plot(c/f*1e9,abs(custom_signal_f), 
             "wavelength (nm)","spectrum","Input spectrum");


plotxy(custom_signal_t*1e15,real(custom_signal)+4,
       response_FDTD_t*1e15,real(response_FDTD)+2,
       response_impulse_t*1e15,response_impulse+0,
       "time (fs)","field");
legend("input signal","response FDTD","response Impulse");