function [outGPS,outEcho] = readHypackraw(filenameRAW,varargin)
%READHYPACKRAW  read Hypack *.raw file and store separate echo.txt and gps.txt
%
%   Syntax:
%   readHypackraw('00X_XXXX.RAW')    % this will create the Echo.txt and GPS.txt for further postprocessing of the jetski survey
%   or 
%   [GPS,Echo] = readHypackraw('00X_XXXX.RAW')  % In addition to the Text files, the data will also be stored in the workspace. 
%
%   Input:
%   filenameRAW  of form ('001_1045.RAW') as created by HYPACK hydrographic surveying package
%
%   Output:
%   varargout = Echo and GPS structures containing the data extracted from the .RAW file
%
%   Example Post processing Jetski Survey.
%   1. readHypackraw('00X_XXXX.RAW')
%   2. Transform 'GPS.txt' to RD/NAP coordinate system and store as 'GPS_RDNAP.txt' using the RWS
%   Coordinate Calculator (free download @ http://www.rdnap.nl/)  Set units to d.dd!
%   3. XYZ=Create_Bathy
%   4. ShowBathy(XYZ)
%
%   See also 
%   Create_Bathy.m 
%
%    DISCLAIMER
%    Hydrographic equipment and software used for this data is purchased for educational use.
%    RTK-GPS corrections are provided by LNRglobal.com for non-commercial use.
%    Please contact TUDelft (M.A. de Schipper & Sierd de Vries) when using this data
%
%   --------------------------------------------------------------------
%   Copyright (C) 2009,2010 <TUDELFT>
%       M.A. de Schipper & Sierd de Vries
%
%       m.a.deschipper@tudelft.nl & Sierd.DeVries@TUdelft.nl
%       
%   This library is free software: you can redistribute it and/or
%   modify it under the terms of the GNU Lesser General Public
%   License as published by the Free Software Foundation, either
%   version 2.1 of the License, or (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library. If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% Created: 16 Jun 2009
% Created with Matlab version: 7.7.0.471 (R2008b)

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $

%% Delete non RTK solution points
%   Addition 18-9-2009 M.A de Schipper
%   Sometimes the GPS solution will drop from accurate RTK solution to
%   inaccurate GPS Standalone. These values will be erased automatically.
%   Expert users ;) might want to see also non RTK points. 
%   Then call syntax: readHypackraw('00X_XXXX.RAW','all') 
if nargin == 2
    if varargin{1}=='all' % then user has specified that he wants also wants the inaccurate pointsif vargin
        OnlyRTK=0;
        warning('GPS.txt also contains non RTK solution points!')
    else
        OnlyRTK=1;
    end
else
    OnlyRTK=1;
end

%%
tic

fid = fopen(filenameRAW);

if fid <=0
    error('Input .RAW file not found, check folder')
end

h = waitbar(0,'Determining size of input file');

%% Locate end of file , necessary to determine the initial length of matrices
gi=0; % set counter
go = 1;
while go          
    tline = fgetl(fid);   
    gi=gi+1;
    if tline <= 2  % end of file
        go = 0;
        i_max=gi-1;
    end
end    
% disp([num2str(i_max) ' lines found, start reading']) 
frewind(fid);

%% Allocate matrices, this speeds up the computation enormous!

Echo.PCtime_s=zeros(i_max,1);
Echo.depth1=zeros(i_max,1);
Echo.depth2=zeros(i_max,1);

GPS.GPStime_hhmmss=zeros(i_max,1); GPS.lat=zeros(i_max,1); GPS.lon=zeros(i_max,1); GPS.Mode=zeros(i_max,1); 
GPS.n_sat=zeros(i_max,1); GPS.HDOP=zeros(i_max,1); GPS.altitude=zeros(i_max,1); GPS.geoid=zeros(i_max,1);GPS.age=zeros(i_max,1); 

VTG.COG_true=zeros(i_max,1); VTG.COG_mag=zeros(i_max,1); VTG.speed_knots=zeros(i_max,1); VTG.speed_kmh=zeros(i_max,1);

%reset counters
ei=0;   %echo sounder strings
gi=0;   %GGA strings
vi=0;   %VTG strings
n=1;    %counter for progress

%% Split .raw file
waitbar(0,h,[num2str(i_max) ' lines found, please wait while reading inputfile.']);
set(h, 'name','Dude... it''s reading, man! ' ) %;-)
for i=1:i_max
    tline = fgetl(fid);
            if strcmp(tline(1:3),'TND')
                %% Format of the date/time line is: TND 10:42:29 07/16/2009 -120 ; 4th field is the timezone (minutes deviation PCtime-UTC) 
                TND=strread(tline ,'%s','delimiter', ' ');
                GPS.Surveydate=datenum(TND{3,1},'mm/dd/yyyy'); % write the survey date in matlab format e.g. 07/16/2009 => 733970 
            end
            if strcmp(tline(1:3),'EC2')  % lines with Echosounder info 
                ei=ei+1;
                %%  Format of GGA line is : EC2 1 38548.426 0.000 3.480 (
                %% 3rd field = HYPACK Time in s after midnight (PC clock time), 4th field= dummy (depth frequency 1), 5th field = depth frequency 2(210kHz)
                 [dum dum Echo.PCtime_s(ei) Echo.depth1(ei) Echo.depth2(ei) ] = ...
                 strread(tline ,'%s %f %f %f %f','delimiter', ' ');
            end
            
            if strcmp(tline(1:3),'MSG')   %% lines starting with MSG indicate GPS output
                if strcmp(tline(20:22),'GGA') && length(tline) >= 87  % GGA lines give GPS position output  
                    gi=gi+1;
                    %tline  %% Uncomment this line to see an output of the lines 
                    %%  Format of GGA line is : MSG 0 66406.052 $GPGGA,162646.00,4444.07068,N,00115.03866,W,2,12,0.8,-1.52,M,49.91,M,3.3,0126*75
                    %% lat and lon are given in DegreeMinutes.Decimalminutes
                    %% GPSmode values: 4 means RTK, 2 means no RTK (differential using SBAS) 
                    [GPS.stamp(gi) GPS.GPStime_hhmmss(gi) GPS.lat(gi) GPS.latsign(gi) GPS.lon(gi) GPS.lonsign(gi)...
                        GPS.Mode(gi) GPS.n_sat(gi) GPS.HDOP(gi) GPS.altitude(gi)...
                        GPS.unit(gi) GPS.geoid(gi) GPS.unit2(gi) GPS.age(gi) GPS.stuff(gi)] = ...
                        strread(tline...
                        ,'%s %f %f %s %f %s %f %f %f %f %s %f %s %f %s','delimiter', ',');
                    if OnlyRTK == 1 ; if GPS.Mode(gi) ~= 4 ; gi=gi-1; end; end %% If wanted, Non RTK values are overwritten, becase the counter is setback a position
                end

                if strcmp(tline(20:22) ,'VTG') && length(tline) >= 20 % VTG lines give velocity and heading output
                    vi=vi+1;
                    %%  Format of VTG line is : MSG 0 66406.095 $GPVTG,0.2,T,356.9,M,4.8,N,8.9,K,D*20
                    %%  COG means Course over ground
                    [VTG.stamp(vi) VTG.COG_true(vi) dum VTG.COG_mag(vi) dum VTG.speed_knots(vi) dum VTG.speed_kmh(vi) dum dum] = ...
                    strread(tline ,'%s %f %s %f %s %f %s %f %s %s','delimiter', ',');                    
                end                 
            end
    if i==round(i_max/100)*n
       waitbar(i/i_max)
       n=n+1;
   end
end

close(h);
fclose(fid);
toc

%% see if the depth values are in the Eho.depth1 or Echo.depth2 position Update M.A. de Schipper nov 2010
if max(Echo.depth1)>0
Echo.depth1=Echo.depth1(1:ei,:);
Echo.depth=Echo.depth1;    
end
if max(Echo.depth2)>0
Echo.depth2=Echo.depth2(1:ei,:);
Echo.depth=Echo.depth2;   
end
if ( max(Echo.depth1)>0 && max(Echo.depth2)>0 ) || ( max(Echo.depth1)==0 && max(Echo.depth2)==0 ) % if values are in both depth1 ann depth2 OR no values in any of the two
disp ('Error... No or too many depth values found!')
Echo.depth=Echo.depth2(1:ei,:);  % just use the zeros so the script continues (handy for car surveys in Hypack)
end
Echo = rmfield(Echo, ['depth1';'depth2'])

%% Recut allocated matrices to right length
Echo.PCtime_s=Echo.PCtime_s(1:ei,:); 
GPS.GPStime_hhmmss = GPS.GPStime_hhmmss(1:gi,:);
GPS.lat = GPS.lat(1:gi,:);
GPS.lon = GPS.lon(1:gi,:);
GPS.Mode = GPS.Mode(1:gi,:);
GPS.n_sat = GPS.n_sat(1:gi,:);
GPS.HDOP = GPS.HDOP(1:gi,:);
GPS.altitude = GPS.altitude(1:gi,:);
GPS.geoid = GPS.geoid(1:gi,:);
GPS.age = GPS.age(1:gi,:);

VTG.COG_true = VTG.COG_true(1:vi,:);
VTG.COG_mag = VTG.COG_mag(1:vi,:);
VTG.speed_knots = VTG.speed_knots(1:vi,:);
VTG.speed_kmh = VTG.speed_kmh(1:vi,:);

%% GPS Structure: Add extra fields 

%%% elevation antenna equals antenna altitude above geoid + Geiodal seperation wrt WGS84 
GPS.elevation=GPS.altitude+GPS.geoid;

%%% Add a minus sign to all longitudes West of Greenwich median
GPS.lon(cell2mat(GPS.lonsign') == 'W')=GPS.lon(cell2mat(GPS.lonsign') == 'W')*-1;

%%% Include time in matlab format 
GPS.Surveydate(1:gi,:)=GPS.Surveydate;  %% Turn single value for the surveydate into a vector to match the length of the GPS.time variable
GPS.Matlabtime = datenum([datestr(GPS.Surveydate(:),'yyyymmdd ') num2str(GPS.GPStime_hhmmss(:),'%06.0f\n')],'yyyymmdd HHMMSS'); %% Matlabtime

%%% Include GPS time in seconds after midnight. Needed to compare the
%%% echosounder time with GPS time
GPS.GPStime_s_midnight = (datenum(num2str(GPS.GPStime_hhmmss),'HHMMSS')-datenum(num2str(000000),'HHMMSS'))*24*60*60;

%%% Hypack time(Pc clocktime)  Needed to link the echosounder time to GPS time 
temp = cell2mat(GPS.stamp'); %Rewrite cell structure 2 string
GPS.PCtime_s_midnight=str2num(temp(:,6:15)) ; % allocate time part of GPS.stampstring to a variable GPS.time_s_midnight

%%% Rewrite 2 format deg.decimaldegrees (nessecary 4 coordinate transformation 2 NAP and deg2utm)
GPS.lat_deg=floor(GPS.lat/100) + (GPS.lat-floor(GPS.lat/100)*100)/60;
GPS.lon_deg=floor(GPS.lon/100) + (GPS.lon-floor(GPS.lon/100)*100)/60;

%%% include RD Coordinates
[GPS.rdnap] = plh_etrs2rdnap([GPS.lat_deg,GPS.lon_deg,GPS.elevation]);

%% ECHO Structure: Tranform Echo sounder time (PC Clock) to GPS time
PCtime2GPStime=GPS.GPStime_s_midnight(1)  -  (GPS.PCtime_s_midnight(1)+str2double(TND{4,1})*60) ; % timetag GPS in GGA - timetag arrival in PC (incl timezone) ; latency will be corrected for later
Echo.GPStime_s=Echo.PCtime_s+str2double(TND{4,1})*60-PCtime2GPStime; %% rewrite the Echosounder time to GPStime based on the timetags of the GPS
Echo.GPSMatlabtime=Echo.GPStime_s/(60*60*24)+GPS.Surveydate(1); %% Rewrite Echosounder time in [UTC] Matlab time 

%%% N.B. Pc Clock deviation over the survey and latency can be evaluated
%%% using the Latency.m script

%% Write output files
fid = fopen('Echo.txt','wt');
fprintf(fid,'Survey of %s \n',datestr(GPS.Surveydate(1)));  %% This is used in Create bathy to have the survey date.. don't change position of string
fprintf(fid,'Hypack Echosounder values extracted using readHypackraw.m from %s \n',filenameRAW);
fprintf(fid,'Processed on %s \n',datestr(now));
fprintf(fid,'%10s \n',['Time UTC [HHMMSS.decSS] ; '  'depth [m]']  );  %% Alignment is poor, can be improved later
fprintf(fid,'%010.3f\t %10.2f \n',[str2num(datestr(Echo.GPSMatlabtime,'HHMMSS.FFF'))' ;  Echo.depth'] );  %%% Echosounderfile is written in this format to keep the time information in decimal sec  
fclose(fid);

fid = fopen('GPS.txt','wt');
fprintf(fid,'Survey of %s \n',datestr(GPS.Surveydate(1))); % This is used in Create bathy to have the survey date.. don't change position of string
fprintf(fid,'Hypack GPS values extracted using readHypackraw.m from %s \n',filenameRAW);
fprintf(fid,'Processed on %s \n',datestr(now));
fprintf(fid,'%10s \n',[ 'GPS Time UTC [HHMMSS.decSS]; ' , 'X        ; ' 'Y        ; ' ' Z']  );  %% Alignment is poor, can be improved later
fprintf(fid, 'Hypackf_%08.0f\t %02.9f\t %02.9f\t %0.5f\n', [GPS.GPStime_hhmmss'*100; GPS.lat_deg' ; GPS.lon_deg' ; GPS.elevation' ]); %% necessary to start the first field (time) with a character so the Coordinate Calculator doesn't pick it as a coordinate field
fclose(fid);

fid = fopen('GPS_RDNAP.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 \n',datestr(GPS.Surveydate(1))); % This is used in Create bathy to have the survey date.. don't change position of string
fprintf(fid,'Hypack GPS values extracted using readHypackraw.m from %s \n',filenameRAW);
fprintf(fid,'Processed on %s \n',datestr(now));
fprintf(fid,'%10s \n',[ 'GPS Time UTC [HHMMSS.decSS]; ' , 'X        ; ' 'Y        ; ' ' Z']  );  %% Alignment is poor, can be improved later
fprintf(fid, 'Hypackf_%08.0f\t %02.3f\t %02.3f\t %0.3f\n', [GPS.GPStime_hhmmss'*100; GPS.rdnap(:,1)' ; GPS.rdnap(:,2)' ; GPS.rdnap(:,3)' ]); %% necessary to start the first field (time) with a character so the Coordinate Calculator doesn't pick it as a coordinate field
fclose(fid);



%% FINISH SCRIPT BELOW 
%% make utm output
%[Northing,Easting,utmzone] = deg2utm(GPSraw(:,3),GPSraw(:,4)); % deg2utm should be included

%% Make function output 
outGPS=GPS;
% 
% outGPS.units.time='s  since midnight in UTC';
% outGPS.units.lat='dd.decdegrees';
% outGPS.units.lon='ddd.decdegrees positive East';
% outGPS.units.elev='m underside antenna relative to geoid';
% outGPS.units.UTM_Northing='m';
% outGPS.units.UTM_Easting='m';
% outGPS.units.GPSmode='1 = standalone, 2=, 3=, 4=RTK';

outEcho=Echo;

% outEcho.units.time='s  since midnight in UTC';
% outEcho.units.sounding='m below transducer';
% %%%




