function [timeseries,years,seriesnames] = rwl2mat(filename)

% RWL2MAT  Read Tucson formatted raw ring width measurements into Matlab
%       
%    [timeseries,years,seriesnames] = rwl2mat(filename)
%
% Input:
%  filename = string with filename and optionally the path
%
% Outputs:
%  timeseries  = an [n x p] matrix containing the raw series in columns
%  years       = vector of years with length [n]
%  seriesnames = cell array containing the code for each core
%

% Change History
%  14 September 2010 - [KJA] modified how the last year was determined so extra empty
%                            rows are no longer being created
%  20 May 2011       - [KJA] modified so that '999' not recognized as stop marker; 
%                            replaced with -999, which might break some 0.01 mm precision rwl files
%                            Thanks to Mandy Freund for spotting this problem
 
% Open existing file in text mode
fid=fopen(filename,'rt');
headerlines = 0;

% read in data
[data]=textscan(fid,'%8s%f%f%f%f%f%f%f%f%f%f%f%*s','TreatAsEmpty',' ','EmptyValue',NaN); 

% what if there is a header? 
while isempty(data{2}) || isempty(data{3})
 headerlines = headerlines + 1;
 [data]=textscan(fid,'%8s%f%f%f%f%f%f%f%f%f%f%f%*s','Headerlines',headerlines,'EmptyValue',NaN);  
end

% Construct matrix with measurement data
datamatrix = [data{3} data{4} data{5} data{6} data{7} data{8} data{9} data{10} data{11} data{12}];

% Warn the user if the value of 999 is found in the data matrix
if sum(sum(any(datamatrix==999))) > 0
    disp('Warning! If value 999 is used as stop marker, this function might fail to extract your file properly')
    pause(2)
end

% find stop markers
decades = data{2}; [cbR,cbC] = find(datamatrix == -999 | datamatrix == -9999);
cbRx = find(decades(cbR)==max(decades));  % pointers to lines with the latest (oldest) decade
minimumYear = min(decades); 
latestCore = max(sum(~isnan(datamatrix(cbR(cbRx),:)),2)); % the longest string of non-missing values in the last decade row
maximumYear = max(decades) + latestCore - 1; % ... since '999' and '-9999' will count as a non-missing value

if ~all(all([datamatrix(cbR(cbRx),:) ~= -999 & datamatrix(cbR(cbRx),:) ~= -9999])) % ... since if all the latest decades simply have the missing value flag, they are actually empty lines as far as data is concerned
 maximumYear = maximumYear - 1;
end

years = [minimumYear:maximumYear]';
timeseries = NaN(length(years),length(cbR));
cbR = sort([cbR]); seriesnames = data{1}(cbR);

% append a zero so we can loop over all the series
cbR = sort([0; cbR]);

% now, read each individual into a timeseries matrix
for i=1:length(cbR)-1;
  thisStartYear = find(years==decades(cbR(i)+1));
  thisCore = datamatrix(cbR(i)+1:cbR(i+1),:)'; thisCore=thisCore(:); thisCore = thisCore(~isnan(thisCore) & (thisCore~=999 & thisCore ~= -9999));
  timeseries(thisStartYear:thisStartYear+length(thisCore)-1,i) = thisCore;
  if years(thisStartYear)+length(thisCore)>maximumYear+1; 
        % [i years(thisStartYear) years(thisStartYear)+length(thisCore)]
        error('Core is longer than maximum time period. Something has gone wrong!'); 
  end
end

fclose(fid);