#####################################
# gain_pc.lsf
# The analysis script file that goes with gain_pc.fsp
#  Sweep Lorentz permittivity  and calculates the slope of the
# log plot of the time signal from the time monitor.
#
######################################
# sweep lorentz gain parameter
clear;
closeall;

plot_timesignal=0;

#number of permittivity values to sweep over
p=11;
permittivity=linspace(0,-0.1,p);
q = matrix(p,1);
n = matrix(p,1);

for(i=1:p){

switchtolayout;
?"Simulation "+num2str(i)+" of "+num2str(p);
setmaterial("Gain_Material", "Lorentz Permittivity", permittivity(i));
run;

    # get center frequency of resonance
    hz = pinch(getmagnetic("freq"));
    f = pinch(getdata("freq","f"));
    shz = size(hz);
    for (j=1:shz(1)) {
        if (hz(j)==max(hz)) {fc = f(j);}
    }
	
	#
    # find slope of time exponential
	# first get data from d-card
	
    t = pinch(getdata("time","t"));
    hz = pinch(getdata("time","Hz"));
	
	#plot the time signal / logplot if plot_signal is set to 1
	if(plot_timesignal){
	plot(t, hz,"t","Hz","lorentz permittivity="+num2str(permittivity(i)));
	plot(t, hz,"t","Hz","lorentz permittivity="+num2str(permittivity(i)),"logplot");
	}
	
	#estimate the max value that occurs near the start and end of the time signal, use them to calculate the slope of the log plot
    totalpoints = size(hz);
    dt = t(2) - t(1);
    T = 1/fc;
    points = round(T/dt);
    hz_start = matrix(points,1);
    hz_end = matrix(points,1);
    
    for (k=1:points) {
        hz_start(k) = hz(k);
        hz_end(k) = hz(totalpoints(1) - points + k);
    }    

    for (l=1:points) {
        if (hz_start(l)==max(hz_start)) {ts = t(l);}
        if (hz_end(l)==max(hz_end)) {te = t(totalpoints(1) - points + l);}
    }

    slope = (log10(max(hz_start)) - log10(max(hz_end))) / (ts - te);
    q(i) = slope;
	
	#get and store the material index for plotting analysis
	n(i)=getindex("Gain_Material",fc);
}

#plot slope vs k
k=imag(n);
plot(k,q,"k","slope of log10(Hz(t))");