% ------------------------------------------------------------
%
% read_tb3 (filename, interval, type)
%
% Function to Read and Process TB3 Rainfall Format
%
% Args:     filename should be with path and enclosed by ' '
%           interval is an integer: 1 = Daily Totals (mm)
%                                   2 = Hourly Totals (mm)
%                                   3 = 30-min Totals (mm)
%           type is the tip size:   1 = 0.2 mm
%                                   2 = 0.254 mm
% 
% Script by S. Yatheendradas and E.R. Vivoni (ASU)
%   Contact E.R. Vivoni (vivoni@asu.edu) for questions
%
% ------------------------------------------------------------
 
function read_tb3(filename, interval, type)

%Open File
fid = fopen(filename, 'rt');

%Read First Line Header (DUR, SITENAME, OTHER INFO)
tline = fgetl(fid);

%Initialize Variables
pulse = [];
index_i = 1;
index_f = 0;
index = 0;

%Read Each Line of File
while feof(fid) == 0
   tline = fgetl(fid);
   
   %Dissect the line for tip time stamp (HH:MM:SS)
   if(length(tline)<15)
       index = index + 1;
       
       %Find the location of the :
       sep = strfind(tline,':');
       
       %Convert date into numerical value for hour/min/sec
       if(length(sep)~=0)
        d = datenum(0,0,0,str2num(tline(1:(sep(1)-1))), str2num(tline((sep(1)+1):(sep(2)-1))), str2num(tline((sep(2)+1):(sep(2)+2))));
        pulse = [pulse d];
       end
       

   %Dissect the line for daily accumulation (HH:MM:SS MM:DD:YY Total)
   else
       index_f = index;
       
       %Find the location of the /
       sep = strfind(tline,'/');
       
       %Assign day and month
       if type == 1
            day = str2num(tline((sep(1)-2):(sep(1)-1)));
            month = str2num(tline((sep(1)+1):(sep(2)-1)));
       end
       if type == 2
            month = str2num(tline((sep(1)-2):(sep(1)-1)));
            day = str2num(tline((sep(1)+1):(sep(2)-1)));
       end
       
       %Determine year
       year = str2num(tline((sep(2)+1):(sep(2)+2))) + 2000;
       
       %Convet date into numerical value
       add = datenum(year,month,day);
       
       %Add year/month/day to hour/min/sec
       pulse(index_i:index_f) = pulse(index_i:index_f) + add;
       
       index_i = index + 1;
   end
end

%Close File
fclose(fid);

%Convert back to date
pulse_2 = datevec(pulse);
lengthpulse = length(pulse_2);

%Read First and Last Dates  
firstDate = datenum(pulse_2(1,:));
lastDate = datenum(pulse_2(lengthpulse,:));

%Determine number of days and hours in record
NumDays = round((lastDate-firstDate)+1);
NumHours = NumDays*24;

%Create file with zeros for different interval

count = 0;

%Daily Interval
if interval == 1
    OutDataMat = zeros(NumDays,4);

    for ii = 1:NumDays
        
        %Convert back to date
        DateVecDay = datevec(firstDate+ii-1);
        
        %Assign Year, Month, Day
        OutDataMat( ii ,1:3) = DateVecDay(1,1:3);
        
        %Find daily number of tips
        for iv = 1:lengthpulse
            if(pulse_2(iv,1:3) == DateVecDay(1,1:3))
                count = count+1; 
            end
        end
        
        %Count tips
        dailycount(ii,1) = count;
        
        %Assign daily rainfall values
        if type == 1
            OutDataMat(ii,4) = dailycount(ii,1)*0.2;
        end
        if type == 2
            OutDataMat(ii,4) = dailycount(ii,1)*0.254;
        end
        
        count = 0;
        
    end

end

%Hourly Interval
if interval == 2
    OutDataMat = zeros(NumHours,5);     
        
    for ii = 1:NumDays
        
        %Convert back to date
        DateVecDay = datevec(firstDate+ii-1);
        
        %Assign Year, Month, Day
        OutDataMat( ((ii-1)*24+1):((ii-1)*24+24) ,1:3) = ones(24,1)*DateVecDay(1,1:3);
        
        %Assign Hour
        OutDataMat((ii-1)*24+1:(ii-1)*24+24,4)=[0:23]';
        
        %Find hourly number of tips
        for ix = 1:24
            for iv = 1:lengthpulse
                
                %Compare tip times and hours
                if(OutDataMat((ii-1)*24+ix,1:4) == pulse_2(iv,1:4))
                    count = count + 1; 
                end 
            end
            
            %Assign total tips
            hourlycount((ii-1)*24+ix,1) = count;
                
            count = 0;
        end
    end

    %Assign hourly rainfall values based on type 
    if type == 1  
        OutDataMat(:,5) = hourlycount(:,1)*0.2;
    end
    if type == 2  
        OutDataMat(:,5) = hourlycount(:,1)*0.254;    
    end
    
end

%30-minute Interval
count30 = 0;
count60 = 0;

if interval == 3
    OutDataMat = zeros(2*NumHours,5);   
         
    for ii = 1:NumDays
        
        %Convert back to date
        DateVecDay = datevec(firstDate+ii-1);
        
        %Assign Year, Month, Day
        OutDataMat( ((ii-1)*48+1):((ii-1)*48+48) ,1:3) = ones(48,1)*DateVecDay(1,1:3);
        
        %Assign Hour
        halfhour = [0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23];
        OutDataMat((ii-1)*48+1:(ii-1)*48+48,4)=halfhour';
        
        %Find half-hourly number of tips
        for ix = 1:24
            for iv = 1:lengthpulse
                
                %Compare tip times and hours
                if(OutDataMat((ii-1)*48+(ix-1)*2+1,1:4) == pulse_2(iv,1:4))
                    if(pulse_2(iv,5) < 30)
                        count30 = count30 + 1;
                    else
                        count60 = count60 + 1;
                    end 
                end 
            end
            
            %Assign total tips
            halfhourlycount((ii-1)*48+(ix-1)*2+1,1) = count30;
            halfhourlycount((ii-1)*48+(ix-1)*2+2,1) = count60;
                
            count30 = 0;
            count60 = 0;
        end   
    end
    
    %Assign hourly rainfall values based on type 
    if type == 1  
        OutDataMat(:,5) = halfhourlycount(:,1)*0.2;
    end
    if type == 2  
        OutDataMat(:,5) = halfhourlycount(:,1)*0.254;    
    end
    
end


%Compose output filename
loc = strfind(filename,'.');
name = filename(1:loc(1)-1);

%Different names for each interval
if interval == 1
    fileout = strcat(name,'_daily.txt');
end
if interval == 2
    fileout = strcat(name, '_hourly.txt');
end
if interval == 3
    fileout = strcat(name, '_halfhourly.txt'); 
end

%Open file
fid = fopen(fileout, 'w');

%Write header
if interval == 1
    strdate = strcat('Year','\t','Month','\t','Day','\t','Rain','\n');
else
    strdate = strcat('Year','\t','Month','\t','Day','\t','Hour','\t','Rain','\n');
end
fprintf(fid,strdate);

%Write data
if interval == 1
    for kk=1:NumDays
        fprintf(fid,'%d\t%d\t%d\t%7.3f\n',OutDataMat(kk,:));
    end
end
if interval == 2
    for kk=1:NumHours
        fprintf(fid,'%d\t%d\t%d\t%d\t%7.3f\n',OutDataMat(kk,:));
    end
end
if interval == 3
    for kk=1:NumHours*2
        fprintf(fid,'%d\t%d\t%d\t%d\t%7.3f\n',OutDataMat(kk,:));
    end
end

%Close file
fclose(fid); 

return