function [transectStruct] = jarkus_transect(filename, originselection)
%JARKUS_TRANSECT   reads Jarkus raw data into Jarkus struct
%
%   [transect] = jarkus_transect(filename, originselection)
%
% See web : <a href="http://www.watermarkt.nl/kustenzeebodem/">www.watermarkt.nl/kustenzeebodem/</a>
% See also: JARKUS_TRANSECT2GRID  , JARKUS_NETCDF2GRID, JARKUS_UPDATEGRID,
%           JARKUS_TRANSECT2NETCDF, JARKUS_GRID2NETCDF

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$

% Read a jarkus transect file

% The file looks like this:

%         15  1965     0     0     0  1210    46
%    150  -2805      160  -3005      170  -4505      180  -9605      190  -9705
%    200 -11605      210 -13305      220 -14705      230 -16205      240 -17005
%  ...
%    550 -34005      560 -33405      570 -33205      580 -33005      590 -34105
%    600 -35205    999999999999    999999999999    999999999999    999999999999

%     The first line is a header line for the following section.
%     Multiple header, section combinations are stored in each file
%     Each header, section combination defines a year of jarkus
%     transect data

%     The header line is
%     areacode, year, metre, dummy, dateTopo, dateBathy, n
%     the numbers in the section are
%     metrering (distance along the shoreline), altitude (surface elevation relative to NAP) with a code as the last digit.,
%     The code is interpreted as
%     Origin of the data:
%     id=1 non-overlap     beach data
%     id=2     overlap     beach data
%     id=3     interpolation     data (between beach and off shore)
%     id=4     overlap off shore data
%     id=5 non-overlap off shore data

%     With originselection control is given which data is transformed
%     to the transectStruct. By default origin 1,3 and 5 are combined

if nargin < 2
    originselection = [1 3 5];
end

fid                         =     fopen(filename);

transectStruct = jarkus_createtransectstruct();
transectStruct = repmat(transectStruct, 1, 500);

i = 0; % count starting by 0
while ~feof(fid)    % Statement is valid until the end of each area data file is reached
    
    %% read index line
    
    line    = fscanf(fid,'%i',7);
    numbers = num2cell(line);
    if (length(numbers) ~= 7)
        error(['invalid line detected:' line]);
    end
    
    %% create a new structure

    transect          = jarkus_createtransectstruct();
    [transect.areaCode, transect.year, transect.alongshoreCoordinate, dummy, dateTopo, dateBathy, transect.n] = deal(numbers{:});
    transect.id       = transect.areaCode*1e6 + transect.alongshoreCoordinate; %TODO: we need unique ids! Maybe find another way
    if transect.id < 0
        fprintf('WARNING: %s\n%s\n', filename, line)
    end
    transect.areaName = rws_kustvak(transect.areaCode);
    
    %% read data
    
    dataPoints        = fscanf(fid,'%f',[2 transect.n])';
    dummy             = fgetl (fid);
    
    % we have to do a bit of cleaning up
    % - dates converting to datenums
    % - altitude -> strip off last digit,
    % - combine different origins and remove non-used data
    
    transect.time          = datenum([transect.year  7  1]) - datenum([1970 1 1]); % assign data to center of year
    transect.timelims(1,:) = datenum([transect.year  1  1]) - datenum([1970 1 1]); % start
    transect.timelims(2,:) = datenum([transect.year 12 31]) - datenum([1970 1 1]); % end
    
    % split up (format is day*100+month)
    transect.dayTopo    = div(dateTopo ,100);
    transect.dayBathy   = div(dateBathy,100);
    transect.monthTopo  = mod(dateTopo ,100);
    transect.monthBathy = mod(dateBathy,100);
    
    %% compute to days since 1970 ...
    
    transect.timeTopo  = datenum([transect.year transect.monthTopo,  transect.dayTopo ])- datenum([1970 1 1]);
    transect.timeBathy = datenum([transect.year transect.monthBathy, transect.dayBathy])- datenum([1970 1 1]);
    
    % the last digit in the 2nd number is a code. See jarkus_createtransectstruct
    % creation for definition
    % we want everything in row vectors because that makes concatting
    % easier and is the default in matlab (think 1:3)
    
    origin                        =   mod(abs(dataPoints(:,2)), 10)';     % use last digit
    altitude                      = floor(abs(dataPoints(:,2))/ 10).*sign(dataPoints(:,2))/100 ; % bathymetry, set last digit to 0, convert from cm+NAP to m+NAP <This works correct for both + and - altitudes.
%   altitude                      =     floor(dataPoints(:,2) / 10)/100 ; % bathymetry, set last digit to 0, convert to m <This causes negative altitudes to be 1cm too low!
    crossShoreCoordinate          =           dataPoints(:,1)';           % cross-shore coordinate
    
    %% strip out data where origin is not 1,3,5
    
    keepIndex                     = ismember(origin, originselection);
    if ~any(keepIndex); continue; end % Do not store empty transects
    transect.origin               =               origin(keepIndex);
    transect.crossShoreCoordinate = crossShoreCoordinate(keepIndex);
    transect.altitude             =             altitude(keepIndex);
    
    %% store the results in the big array
    transectStruct(i+1)           = transect;
    i = i + 1;
    
end


transectStruct(i+1:end) = [];

end % end of function jarkus2transect
