function GPS_RDNAP = Read_Cache(varargin)

%%% This script merges two datafiles and eventually generates XYZ and
%%% accuracies from the lacation of the GPS antenna.

%%% We should have two datafiles for input
%%% 1) GPS_corrected.txt holds the corrected GPS values generated by LGO
%%%     the coordinate system should be ETRS-Cartesian. This file also
%%%     holds data on accuracies
%%% 2) GPS_corrected_NAP.txt holds the corrected GPS values in the RD
%%%     coordinate system this file is generated using file (1) and the
%%%     rijkswaterstaat coordinate calculator.

% initiate data structure
GPS_RDNAP = struct('GPStime',[],'X_RD',[],'Y_RD',[],'Z_RD',[],'PosQlt',[],...
    'HgtQlt',[],'time',[],'Matlabtime',[]);
filesize = [];

%%% first check if there are GPS files, otherwise abort
if ~exist('GPS_Corrected.txt') &&  ~exist('GPS_RDNAP.txt') &&  ~exist('Leica_GPS_RDNAP.txt')
    error('No GPS files found in this folder, check folder and naming convention (names should be Leica_GPS_RDNAP.txt,GPS_RDNAP.txt or GPS_Corrected.txt)')
end

%%% there are four possibilities..
%1. there are no postprocessed GPS files then use the GPS_RDNAP.txt
%2. there are postprocessed GPS files, including RD NAP then use GPS_corrected_NAP.txt (and GPS_corrected.txt for the accuracy values)
%3. There is a postprocessed GPS file but no RD NAP file.
%4. There is a RTK Leica file exported using standard procedure

%%% Sierd, Aug 2010, include RDNAP procedure on GPS_Corrected. Therefore no
%%% external program is needed to convert to RDNAP coordinates is needed.
%%% Sierd, Nov 2010, include reading GPS data from new Leica receiver

% 3
if ~exist('GPS_Corrected_NAP.txt') && exist('GPS_Corrected.txt')
    disp('Post processed GPS file available, reading file')
    [GPS_RDNAP,filesize] = read_GPSCorrected(GPS_RDNAP,filesize);
    %%% Convert and add RD coordinates
    disp('Converting to RD coordinates')
    for i=1:length(GPS_RDNAP.X_cart)
        rdnap(i,1:3)=etrs2rdnap([GPS_RDNAP.X_cart(i)...
            GPS_RDNAP.Y_cart(i) GPS_RDNAP.Z_cart(i)])';
    end
    GPS_RDNAP.X_RD=rdnap(:,1)';
    GPS_RDNAP.Y_RD=rdnap(:,2)';
    GPS_RDNAP.Z_RD=rdnap(:,3)';

    %%% Write to GPS_Corrected_NAP.txt
    disp('Writing GPS_Corrected_NAP.txt')
    fid = fopen('GPS_Corrected_NAP.txt','wt');
    % fprintf(fid,'GPS coordinates converted using new RDNAPTRANS 2008 procedure \n');
    % fprintf(fid,'GPS coordinates converted using new RDNAPTRANS 2008 procedure \n');
    fprintf(fid,'GPS coordinates converted using new RDNAPTRANS 2008 procedure \n');
    fprintf(fid,'Survey of %s ',datestr(GPS_RDNAP.Surveydate(1))); % This is used in Create bathy to have the survey date.. don't change position of string
    fprintf(fid,'Processed on %s \n',datestr(now));
    % fprintf(fid,'GPS values extracted using GPS_Corrected.txt');
    fprintf(fid,'%10s \n',[ 'GPS Time UTC [HHMMSS.decSS]; ' , 'X        ; ' 'Y        ; ' ' Z']  );  %% Alignment is poor, can be improved later
    fprintf(fid, 'septent_%08.0f\t %02.3f\t %02.3f\t %0.3f\n', [GPS_RDNAP.time*100 ; GPS_RDNAP.X_RD(1,:) ;  GPS_RDNAP.Y_RD(1,:) ; GPS_RDNAP.Z_RD(1,:) ]); %% necessary to start the first field (time) with a character so the Coordinate Calculator doesn't pick it as a coordinate field
    fclose(fid);

%2
elseif exist('GPS_Corrected_NAP.txt') && exist('GPS_Corrected.txt')
    disp('Post processed GPS files including RDNAP available, reading files')
    [GPS_RDNAP,filesize] = read_GPSCorrected(GPS_RDNAP,filesize);
    [GPS_RDNAP,filesize] = read_GPSCorrected_NAP(GPS_RDNAP,filesize);

%4
elseif exist('Leica_GPS_RDNAP.txt')
    disp('Leica Receiver GPS output available, reading file')
    [GPS_RDNAP,filesize] = read_Leica(GPS_RDNAP,filesize);
% 1
else
    disp('No post processed GPS files available, reading Hypack output files')
    [GPS_RDNAP,filesize] = read_GPS_RDNAP(GPS_RDNAP,filesize);
end

%% Include GPS time in seconds after midnight.
GPS_RDNAP.GPStime(1,1:filesize) = (datenum(num2str(GPS_RDNAP.time(:)/100),'HHMMSS')-datenum(num2str(000000),'HHMMSS'))*24*60*60;

%%% Include time in matlab format, this mambo jambo was a puzzle getting right
date = datenum(GPS_RDNAP.Surveydate,'mm/dd/yyyy');
GPS_RDNAP.Matlabtime(1,1:filesize) = datenum([datestr(date(:),'yyyymmdd ') num2str(GPS_RDNAP.time(:)/100,'%06.0f\n')],'yyyymmdd HHMMSS');

GPS_RDNAP = rmfield(GPS_RDNAP,'Surveydate');


%% END OF FILE... Functions Below



%% Read GPS_corrected (File generated by Leico Geo Office)
function [GPS_RDNAP,filesize] = read_GPSCorrected(GPS_RDNAP,filesize)
Filenames.GPS_Corrected='GPS_Corrected.txt';
fid = fopen(Filenames.GPS_Corrected);
gi=0; % set counter
go = 1;
while go
    tline = fgetl(fid);
    if length(tline) >= 50
        gi=gi+1;

        %%% replace annoying '_' in string by a tab
        tline = strrep(tline,'_','	');

        %%% Read tline to structure
        %%% tline info: Point Id,Epoch,Stored Status,Ambiguity Status,GNSS Type,Type,Solution type,Frequency,X,Y,Z,Posn. Qlty,Hgt. Qlty,Posn. + Hgt. Qlty
        [dum GPS_RDNAP.time(gi)... %
            GPS_RDNAP.Surveydate(gi) dum dum...
            dum dum dum...
            GPS_RDNAP.X_cart(gi) GPS_RDNAP.Y_cart(gi) GPS_RDNAP.Z_cart(gi) ...
            GPS_RDNAP.PosQlt(gi) GPS_RDNAP.HgtQlt(gi) dum] =...
            strread(tline, '%s %f %s %s  %s %s %s %s  %f %f %f %f  %f %f', 'delimiter', '	');

        %%% Include GPS time in seconds after midnight.
        %         GPS_Corrected.GPStime(gi) = datenum([0 0 0 floor(GPS_Corrected.time(gi)/1e6)...
        %                 floor((GPS_Corrected.time(gi)-floor(GPS_Corrected.time(gi)/1e6)*1e6)/1e4)...
        %                 floor((GPS_Corrected.time(gi)-floor(GPS_Corrected.time(gi)/1e4)*1e4)/1e2)])*24*60*60;
    else
        go = 0;
        %check if files are of equal length
        if isempty(filesize)
            filesize=gi;
        elseif filesize ~= gi
            error('Files are not of equal length and therefore not compatible')
        end
    end

end
fclose(fid);
% GPS_RDNAP.Surveydate(1,1:gi)=datenum(Surveydate,'mm/dd/yyyy');

%% Read GPS_corrected_NAP.txt (file generated by RWS Coordinate Calculator from LGO file)
function [GPS_RDNAP,filesize] = read_GPSCorrected_NAP(GPS_RDNAP,filesize)
Filenames.GPS_Corrected_RDNAP = 'GPS_Corrected_NAP.txt';
fid = fopen(Filenames.GPS_Corrected_RDNAP);
for i=1:3;tline = fgetl(fid);end %(Get rid of header lines)

gi=0; % set counter
go = 1;

while go
    tline = fgetl(fid);
    if length(tline) >= 47
        gi=gi+1;

        %%% Replace annoying comma's and '_'
        tline = strrep(tline,',','.');
        tline = strrep(tline,'_',' ');

        %%% read tline to structure
        %%% tline info: Point Id, Epoch, X_RD, Y_RD, NAP
        [dum dum GPS_RDNAP.X_RD(gi)...
            GPS_RDNAP.Y_RD(gi) GPS_RDNAP.Z_RD(gi)] = strread(tline, '%s %f %f %f %f', 'delimiter', ' ');

        %         GPS_Corrected_RDNAP.GPStime(gi) = datenum([0 0 0 floor(GPS_Corrected_RDNAP.time(gi)/1e6)...
        %                         floor((GPS_Corrected_RDNAP.time(gi)-floor(GPS_Corrected_RDNAP.time(gi)/1e6)*1e6)/1e4)...
        %                         floor((GPS_Corrected_RDNAP.time(gi)-floor(GPS_Corrected_RDNAP.time(gi)/1e4)*1e4)/1e2)])*24*60*60;
    else
        go = 0;
        %check if files are of equal length
        if isempty(filesize)
            filesize=gi;
        elseif filesize ~= gi
            error('Files are not of equal length and therefore not compatible')
        end
    end

end
fclose(fid);

%% Read GPS_RDNAP.txt (file generated by RWS Coordinate Calculator from Hypack output)
function [GPS_RDNAP,filesize] = read_GPS_RDNAP(GPS_RDNAP,filesize)
Filenames.GPS_unCorrected_RDNAP = 'GPS_RDNAP.txt';
fid = fopen(Filenames.GPS_unCorrected_RDNAP);
for i=1:3;tline = fgetl(fid);end %(Get rid of header lines)
tline = fgetl(fid);  Surveydate=tline(11:21); % read surveydate see end of function
for i=1:3;tline = fgetl(fid);end %(Get rid of header lines)
gi=0; % set counter
go = 1;

while go
    tline = fgetl(fid);
    if length(tline) >= 40
        gi=gi+1;

        %%% Replace annoying comma's and '_'
        tline = strrep(tline,',','.');
        tline = strrep(tline,'_',' ');

        %%% read tline to structure
        %%% tline info: Point Id, Epoch, X_RD, Y_RD, NAP
        [dum GPS_RDNAP.time(gi) GPS_RDNAP.X_RD(gi)...
            GPS_RDNAP.Y_RD(gi) GPS_RDNAP.Z_RD(gi)] = strread(tline, '%s %f %f %f %f', 'delimiter', ' ');

        %         GPS_Corrected_RDNAP.GPStime(gi) = datenum([0 0 0 floor(GPS_Corrected_RDNAP.time(gi)/1e6)...
        %                         floor((GPS_Corrected_RDNAP.time(gi)-floor(GPS_Corrected_RDNAP.time(gi)/1e6)*1e6)/1e4)...
        %                         floor((GPS_Corrected_RDNAP.time(gi)-floor(GPS_Corrected_RDNAP.time(gi)/1e4)*1e4)/1e2)])*24*60*60;
    else
        go = 0;
        %check if files are of equal length
        if isempty(filesize)
            filesize=gi;
        end
    end

end
fclose(fid);

%this takes a lot of time and should be better.
for i = 1:gi
    GPS_RDNAP.Surveydate{1,i}=datestr(Surveydate,'mm/dd/yyyy');
end

function  [GPS_RDNAP,filesize] = read_Leica(GPS_RDNAP,filesize)
Filenames.Leica_GPS_RDNAP = 'Leica_GPS_RDNAP.txt';
fid = fopen(Filenames.Leica_GPS_RDNAP);
gi=0; % set counter
go = 1;
while go
    tline = fgetl(fid);

    if length(tline) >= 40
        %     replace ':'
    tline = strrep(tline,':','');
        gi=gi+1;
        
%         if gi == 1
%                 tline = strrep(tline,'/','');
%                     [dum Surveydate(gi) dum ...
%             dum dum dum ...
%             dum dum ] = ...
%             strread(tline, '%s %f %f %f %f %f %f %f','delimiter','	');
%         end
                tline = strrep(tline,'/','');
        [dum Surveydate(gi) GPS_RDNAP.time(gi) ...
            GPS_RDNAP.X_RD(gi) GPS_RDNAP.Y_RD(gi) GPS_RDNAP.Z_RD(gi) ...
            GPS_RDNAP.PosQlt(gi) GPS_RDNAP.HgtQlt(gi) ] = ...
            strread(tline, '%s %f %f %f %f %f %f %f','delimiter','	');
    else
        go = 0;
        %check if files are of equal length
        if isempty(filesize)
            filesize=gi;
        end
    end
end

GPS_RDNAP.time=GPS_RDNAP.time*100;

%this takes a lot of time and should be better see also previous
%for similar coommand
for i = 1:gi
    GPS_RDNAP.Surveydate{1,i}=datestr(datenum(num2str(Surveydate(i)),'ddmmyyyy'),'mm/dd/yyyy');
end
        