% This program calculates the onset from temperature data. 
%The program is developed using the wavelet program that is freely
%available from "http://paos.colorado.edu/research/wavelets/" which was written January 1998 by C. Torrence

% This program has to be in the folder wave_matlab and program saves the onset marker in a file with a file name starting with Onset



close all;
clear all;
clc;



[FileName,PathName]=uigetfile('*.txt','Open temperature data file (sample file is SampleTemp.txt) ');

data=load([PathName FileName]);
 
sst2=data;


ofile=['Onset_' FileName];

variance = std(sst2)^2;
m_sst=mean(sst2);
sst = (sst2 - m_sst)/sqrt(variance) ;


n = length(sst);
% dt is the sampling time, here it is 15 minutes
dt = 15/60;
time = [0:length(sst)-1]*(dt/24);  
xlim = [0 Inf];  
pad = 1;      
dj=0.25/64;
s0 = 2*dt;
j1 = 7/dj;    
lag1 = 0.72;  



mother='DOG';

% Wavelet transform:
[wave,period,scale,coi] = wavelet(sst,dt,pad,dj,s0,j1,mother);
power = abs(wave).^2 ;        % compute wavelet power spectrum
Cdelta=3.476;
AP_real=real(wave);

as=find((period>=0)&(period<=100));
red_power=power(as,:);
red_period=period(as);
red_AP_real=AP_real(as,:);
red_scale=scale(as);

for i=1:length(time)
    
    I_pksA=zeros(1,5);
    [pks,loc]=findpeaks(red_power(:,i));
    [s_pks,I_pks]=sort(pks,'descend');
    
    for ki=1:length(I_pks);
        I_pksA(ki)=I_pks(ki);
    end;
    
    
     if(I_pksA(1)==0)
         
         [maxpower posfreq1]=max(red_power(:,i));
              
     else
                  
   posfreq1=loc(I_pks(1));
   maxpower=pks(I_pks(1));
     end;
     
     cf11=(variance*dt)/(Cdelta*red_scale(posfreq1));
     cf1=sqrt(cf11);
     insamp(i)=cf1*red_AP_real(posfreq1,i);
     
     
      
end;

kp=1;

time_minp(1)=0;
dtmin=1;

for i=2:length(insamp)-1
  
    
    
    if(((insamp(i-1)<=0)&&(insamp(i)> 0)))
        
       minp(kp)=0; %#ok<AGROW>
       timeOnset(kp)=time(i);
      
               kp=kp+1; 
       
    end;
    
    
        %#ok<AGROW>
    
end;


kpp=2;
timeOnsetmod(1)=timeOnset(1); 
minpmod(1)=0;

dtime=timeOnset(2)-timeOnset(1);
for j=2:length(timeOnset)-1
    
   
    
    if(dtime>0.9)
        
        timeOnsetmod(kpp)=timeOnset(j);
        minpmod(kpp)=0;
        kpp=kpp+1;
        
    end;
    
    dtime=timeOnset(j+1)-timeOnsetmod(kpp-1);
    
end;

figure(1)
plot(time,insamp,'k');
hold on;
plot(timeOnsetmod,minpmod,'r*');
axis([0 Inf -1.2 1.2]);

fid=fopen(ofile,'w');
fprintf(fid,'%f \n',timeOnsetmod);
status=fclose(fid);




