tStart = 100e-15;
tStop = 500e-15; #short simulation
#tStop = 6000e-15; #long simulation

# ---------------------------------------

# Initial frequency range:
fmin = c/450e-9;
fmax = c/400e-9;
# ---------------------------------------

# Frequency range to get 418nm resonance: 
#fmin = 715e12; 
#fmax = 719e12;
# ---------------------------------------  

# Frequency range to get 428nm resonance: 
#fmin = 695e12;
#fmax = 705e12;
# ---------------------------------------

logfilename = "multiple_resonances.log";
numSteps = 51;
monitorName = "timemonitor1";

# ---------------------------------------

if(fileexists(logfilename)){rm(logfilename);}
logmessage = "t_max (fs), resonance freq. (THz), error estimate" + endl;
write(logfilename, logmessage);


t = getresult(monitorName,"t");
curComp = "Hz";
rawSignal = pinch(getresult(monitorName,curComp));
 
tFilter = find( (t>tStart) & (t<tStop) );
tFiltered = t(tFilter);
rawSignal = rawSignal(tFilter);
ptsPerInterval = floor(length(tFiltered)/numSteps);
 
for(stopIdx = 1:numSteps) {
    curT = tFiltered(stopIdx*ptsPerInterval);
    tFilter2 = find( (tFiltered>tStart) & (tFiltered<curT));
    tFiltered2=t(tFilter2);
    sFiltered=rawSignal(tFilter2);
    
    result = findresonances(tFiltered2,sFiltered,[fmin,fmax]);    
    
    if( length(result)>0 ) {
        THz_frequ = result(:,1)/1e12;
        fs_time = curT*ones(length(THz_frequ),1)/1e-15;
        error_estimate = result(:,6);      
        
        size_result = size(result);
        num_res = size_result(1); # number of resonances found
              
        if(num_res>1){
            logmessage = num2str([fs_time, THz_frequ, error_estimate]) + endl;                        
            write(logfilename,logmessage);
        }        
        
        plot(fs_time,THz_frequ,"t_max (fs)","f (THz)","Resonance Frequencies "+monitorName,"plot points");
        setplot("show legend", 0);
        holdon;
    }
}
holdoff;
 
for(stopIdx = 1:numSteps) {
    curT = tFiltered(stopIdx*ptsPerInterval);
    tFilter2 = find( (tFiltered>tStart) & (tFiltered<curT));
    tFiltered2=t(tFilter2);
    sFiltered=rawSignal(tFilter2);
    
    result = findresonances(tFiltered2,sFiltered,[fmin,fmax]);
    
    if( length(result)>0 ) {
        Q = result(:,3);
        fs_time = curT*ones(length(Q),1)/1e-15;
        
        plot(fs_time,Q,"t_max (fs)","Q","Q-factors "+monitorName,"plot points"); #!!! REMOVE FACTOR OF 2 ONCE CODE IS FIXED!
        setplot("show legend", 0);
        setplot("log10y",1);
        holdon;
    }
}
holdoff;