# this script shows how to create a custom time signal

#####################################################
# specify signal

# Use this code to define a time signal in this script.
if (1) {
  t = linspace(0,600e-15, 10000);# the max. time of this vector will affect when does the auto shutoff kicks in. Specify the smallest t_max possible.
  frequency= 300e12;
  w_center = frequency*2*pi;   # 300 THz
  pulselength = 50e-15;
  delta_t = pulselength/(2*sqrt(log(2)));
  offset1   = 150e-15;
  offset2   = 300e-15;

  # sinusoid with double gaussian envelope 
  pulse = sin( -w_center*(t-offset1)) * exp( -(t-offset1)^2/2/delta_t^2 ) + 
          sin( -w_center*(t-offset2)) * exp( -(t-offset2)^2/2/delta_t^2 );
  
  # use this code to write the time signal to a text file
  if (0) { write("time.txt",num2str([t, pulse])); }
} 

# Use this code to read the time signal from a text file
if (0) {
  data  = readdata("time.txt");
  t     = pinch(data,2,1);
  pulse = pinch(data,2,2);
} 

# plot specified time signal
plot(t*1e15,pulse,"time (fs)","signal","Specified pulse");



#######################################################
# convert into complex valued format required by source
# make signal complex valued

# fft time signal
pulse_w=fft(pulse);
w=fftw(t);

# set high frequency components to zero
# and compensate amplitude
n=length(pulse_w);
pulse_w((n/2+1):n)=matrix(n/2); 
pulse_w=pulse_w*2;

# inverse fft the signal
pulse_final=invfft(pulse_w);
t3=fftw(w);

# apply a required phase shift
pulse_final=pulse_final*1i;  


#######################################################
# load signal into source

switchtolayout;
setsourcesignal("custom_source", t3, abs(pulse_final),unwrap(angle(pulse_final)));


#######################################################
# run simulation and compare simulation results with 
# the original signal 
run;


#######################################################
# get simulation data and compare results

# get monitor data.  This signal should be delayed by 
# ~1/4 wavelength, due to monitor position. 
t2=getdata("time","t");
Ez_monitor=getdata("time","Ez");

# get source signal
# Note: this data is sampled at a lower resolution than 
# the original data (sample rate based on nyquist)
t3=getdata("custom_source","time");
source_signal=getdata("custom_source","time_signal");

# plot results
plotxy(t*1e15,pulse,t3*1e15,source_signal,t2*1e15,Ez_monitor,
           "time (fs)","signal","Comparison");
legend("Specified","Source signal","Monitor field (1/4 lambda gap)");