###############################################
# file: usr_spectrum_custom.lsf
#
# This script file will create a real time signal 
# with a power spectrum specified by the user


###############################################
# START OF SETUP SECTION

# Specify length of time signal
# Longer times will allow the spectrum to be more 
# accurately reproduced, but require a longer
# simulation time.
max_time = 30e-15;

# Specify power spectrum
f=linspace(0,1000e12,1000);
spectrum=matrix(1000);

# sloped top hat example spectrum
f0=find(f,250e12);
f1=find(f,300e12);
f2=find(f,600e12);
f3=find(f,650e12);
spectrum(f0:f1)=linspace(0,1,f1-f0+1);
spectrum(f1:f2)=1+matrix(f2-f1+1);
spectrum(f2:f3)=linspace(1,0,f3-f2+1);

## sin function shape example spectrum
# f1=find(f,300e12);
# f2=find(f,600e12);
# spectrum(f1:f2)=sin(linspace(0,pi,f2-f1+1));

# plot the spectrum
plot(f/1e12,spectrum,"f (THz)","spectrum");
legend("Specified spectrum");

# END OF SETUP SECTION
###############################################


###############################################
# calculate time signal
fmax = max(f*(spectrum>1e-8*max(spectrum)));
dt = 1/fmax/200;
t = 0:dt:max_time;
w = fftw(t);
f2=w/2/pi;  # frequency vector based on fft results
temp = interp(sqrt(spectrum),f,f2)*exp(1i*w*max_time/2);
t = fftw(w);

signal=invfft(temp);
signal=signal/max(abs(signal));

p1 = find(t,max_time);
t = t(1:p1);
signal = signal(1:p1);
# end of calculate time signal


###############################################
# plot results
plot(t*1e15,signal,"t (fs)","signal");
legend("Calculated signal");

spectrum_calc=real(signal);
spectrum_calc=abs(fft(spectrum_calc))^2;
spectrum_calc=spectrum_calc/max(spectrum_calc);
plot(f2/1e12,interp(spectrum,f,f2),spectrum_calc,"f (THz)","spectrum");
legend("Specified spectrum","Calculated spectrum");


###############################################
# load time signal into source
switchtolayout;
setsourcesignal("source1", t, abs(signal),unwrap(angle(signal)));


###############################################
# run simulation
run;


###############################################
# compare with fdtd results 


# use source  data and source power function
source_t=getdata("source1","time");
source_t_signal=getdata("source1","time_signal");

nonorm;  # temporarily disable CW norm
spectrum_sourcepower=sourcepower(f);
cwnorm;  # re-enable CW norm
spectrum_sourcepower=spectrum_sourcepower/max(spectrum_sourcepower);



# check with time monitor.  
monitor_t = getdata("time","t");
monitor_t_signal =     getdata("time","Ez");
monitor_f=1/2/pi*fftw(monitor_t);
spectrum_monitor= abs(fft(monitor_t_signal))^2;
spectrum_monitor= spectrum_monitor/max(spectrum_monitor);
spectrum_monitor= interp(spectrum_monitor,monitor_f,f);



plot(source_t*1e15,source_t_signal,"t (fs)","Signal");
legend("FDTD signal");

plot(f/1e12,spectrum,
            interp(spectrum_calc,f2,f),
            spectrum_sourcepower,
            spectrum_monitor,
            "f (THz)","spectrum");
legend("Specified spectrum","Calculated spectrum","Actual source spectrum","Monitor measurement");