function [downcast, upcast] = loadMVPraw(path2file)
%
% =========================================================================
% loadMVPraw  Version 1.1 25-Jun-2013
%
% Usage:
%   [downcast, upcast] = loadMVPraw(path2file)
%
% Description:
%   This MATLAB function opens and reads data from an MVP .raw data file.
%
% Input:
%   An MVP .raw processed data filename
%
% Output:
%   downcast - a structure containing all downcast MVP data
%   upcast - a structure containing all upcast MVP data
%
% Authors:
%   Blair Greenan
%   Bedford Institute of Oceanography
%   28-Dec-2011
%   Matlab 7.13
%   blair.greenan@dfo-mpo.gc.ca
%
%   Adapted by Susanne Craig
%   Bedford Institute of Oceanography
%   18-Jul-2012
%
%   Modified by Dave Hebert
%   dave.hebert@dfo-mpo.gc.ca
%   25 June 2013
%
%   Modified by Susanne Craig
%   susanne.craig@dfo-mpo.gc.ca
%   16 Oct 2013
%
% Revision History:
%
%   18 Jul 2012 -  Adapted from Blair Greenan's load_MVP_m1_convert2vpr.m
%                  function. The last part of this original file was
%                  unnecessary as it interfaced with a VPR visualisation
%                  tool, and has simply been omitted.
%
%   25 June 2013 - Change the arbitary 6 minutes of elapsed time used to
%                  calculate upcast location and time to time based on data
%                  sampled at 25 Hz.  Time is start of upcast but could
%                  easily be modified for middle of upcast. Added
%                  extraction of analog signal, calculation of chl and
%                  correction of LOPC to be counts/sec
%
%   16 Oct 2013 -  Changed the field name for fluorescence from 'anal'
%                  (?!!) to 'analogV'
%                  
%                  Corrected the calculation of Chl from voltage.
%                  Previously, this was set as simple multiplication by 6.
%                  The correct formual is:
%                  Chl = 6 * (analogV - 0.065)
%                  (Source: Blair Greenan's load_MVP_m1_convert2vpr' file)
%
% =========================================================================

%% Parse out latitude

% Open the MVP .m1 file which is a processed data file generated by the
% MVP data logging software
fid = fopen(path2file);

% Bottom Depth from ship's sounder. Units - meters
downcast.bottomDepth = cell2mat(textscan(fid, 'Bottom Depth (m): %f', 'HeaderLines', 6));
upcast.bottomDepth = downcast.bottomDepth;

% Ship's speed. Units - knots
ShipsLog = cell2mat(textscan(fid, 'Ship''s Log (kt): %f'));

% Latitude at start of MVP cast - read as a string and then convert
LatString = textscan(fid, 'LAT ( ddmm.mmmmmmm,N):  %12c,%c');

% Convert first two characters to degrees latitude
LatDeg_downcast = str2double(LatString{1}(1:2));

% Convert characters 3-12 to decimal minutes
LatMin_downcast = str2double(LatString{1}(3:12));

% Convert to decimal degrees
downcast.lat = LatDeg_downcast + (LatMin_downcast/60);

% Record whether latitude is North or South
LatNS = LatString{2};
if strcmpi(LatNS,'S')
    downcast.lat = -1.*downcast.lat;
end

%% Parse out longitude

% Longitude at start of MVP cast - read as a string and then convert
LonString = textscan(fid, 'ON (dddmm.mmmmmmm,E): %13c,%c');

% Convert first three characters to degrees longitude
LonDeg_downcast = str2double(LonString{1}(1:3));

% Convert characters 4-13 to decimal minutes
LonMin_downcast = str2double(LonString{1}(4:13));

% Convert to decimal degrees
downcast.lon = LonDeg_downcast + (LonMin_downcast/60);

% Record whether longitude is East or West
LonEW = LonString{2};
if strcmpi(LonEW,'W')
    downcast.lon = -1.*downcast.lon;
end

%% Parse out time

% Time at start of MVP cast - UTC
TimeArray = textscan(fid, 'Time (hh|mm|ss.s): %2d:%2d:%f');
THour = double(TimeArray{1});
TMinute = double(TimeArray{2});
TSecond = double(TimeArray{3});
clear TimeArray

% Date at start of MVP cast - UTC
DateArray = textscan(fid, 'Date (dd/mm/yyyy): %2d/%2d/%4d');
DDay = double(DateArray{1});
DMonth = double(DateArray{2});
DYear = double(DateArray{3});
clear DateArray

% Make a datenumber
downcast.daten = datenum([DYear DMonth DDay THour TMinute TSecond]);

% Extract the ship's heading from the GPS NMEA string
ShipsHeading = textscan(fid,'BOT F10:$GPVTG,%f', 'HeaderLines', 20);

%% Load the data

% Use this method instead of the readtext function because the space
% delimiter in the raw files is inconsistent.
fseek(fid,0,'eof');
eof = ftell(fid);
fseek(fid,0,'bof');

file_pointer = 0;
cnt = 0;
a = cell(1e6,1);
while file_pointer < eof
    cnt = cnt+1;
    tline = fgetl(fid);
    a{cnt,1} = tline;
    file_pointer = ftell(fid);
end
fclose(fid);

% Remove excess elements
a = a(1:cnt);

hidx = find(strcmp('**********', a), 10, 'last');
% This occurs several times in the header. Choose the last occurrence to
% mark the end of the header.
hidx = max(hidx);

dat = a(hidx+3:end,:);

% There are sometimes rows with differing number of columns. So, first
% check that there are 5 columns, then convert to a number. If not, enter a
% row of NaNs
DataMatrix = zeros(length(dat), 5);
for n = 1:length(dat)
    rowData = str2num(dat{n,1});
    if length(rowData) == 5
        DataMatrix(n,:) = rowData;
    else
        DataMatrix(n,:) = nan(1,5);
        %DataMatrix = [DataMatrix; rowData];
    end
end

% Now get rid of NaNs
idx = find(all(isnan(DataMatrix), 2));
DataMatrix(idx,:) = [];

% Find the maximum depth and use this to divide the array into up cast and
% downcast
pres = DataMatrix(:,1);             % Units: dbar
cond = DataMatrix(:,2);             % Units: mS cm-1
temp = DataMatrix(:,3);             % Units: degrees Celsius
analogV = DataMatrix(:,4) / 1000;   % Units: V (raw units in mV)
lopc = DataMatrix(:,5);             % Units: counts (0.5 s)-1
findmax = find(pres == max(pres));

% There may be several maximum depths. If so, include them in the downcast.
if length(findmax) > 1
    findmax = max(findmax);
end

% Downcast
downcast.pres = pres(1:findmax);
downcast.cond = cond(1:findmax);
downcast.temp = temp(1:findmax);
downcast.analogV = analogV(1:findmax);
downcast.chl = 6 * (downcast.analogV - 0.065); % Analogue range 1 calibration
downcast.lopc = lopc(1:findmax) * 2;           % Convert to counts s-1
idx = find(downcast.lopc == 0);                % Set LOPC values of zero to NaN 
downcast.lopc(idx) = NaN;

%% Calculate upcast co-ords

% Based on the header information, estimate a lat/lon location for the
% upcast that occurs (findmax/25)/60 minutes after the downcast. Estimate
% the distance traveled (km) in x minutes based on the the ship speed.
% 1.852 is the conversion factor from nautical miles to kilometers
xkm = (ShipsLog * 1.852)*(findmax/25.)/60./60.;

% Convert to radians and then to degrees. 6371 is the Earth's radius in km
xdeg = (xkm/6371)*(180/pi);

% Change in latitude & longitude over the x minutes
delta_lon = xdeg*sin(ShipsHeading{1}*(pi/180));

% Adjust the delta_lon based on the latitude
% Lat_downcast = (sign(LatDeg_downcast)*(abs(LatDeg_downcast) + (LatMin_downcast/60))); % decimalise latDown
delta_lon = delta_lon/(cos(downcast.lat * pi/180));
delta_lat = xdeg*cos(ShipsHeading{1}*(pi/180));

% Compute the latitude and longitude of the upcast based on the ship's
% speed and heading and assuming the upcast happens 6 minutes after the
% downcast
upcast.lat = downcast.lat + delta_lat;
upcast.lon = downcast.lon + delta_lon;

% Calculate a datenumber 6 mins later than downcast
upcast.daten = downcast.daten + (findmax/25.)/60./60./24.;

% Upcast
upcast.pres = pres(end:-1:findmax+1);
upcast.cond = cond(end:-1:findmax+1);
upcast.temp = temp(end:-1:findmax+1);
upcast.analogV = analogV(end:-1:findmax+1);
upcast.chl = 6 * (upcast.analogV - 0.065); % Analogue range 1 calibration
upcast.lopc = lopc(end:-1:findmax+1)*2;    % Convert to counts s-1
idx = find(upcast.lopc == 0);              % Set LOPC values of zero to NaN 
upcast.lopc(idx) = NaN;

