function completed = load_MVP_m1_convert2vpr(fname,logfile)
%
% =========================================================================
% load_MVP_m1_convert2vpr  Version 1.0 28-Dec-2011
%
% Usage:
%   load_MVP_m1_convert2vpr(filename,vprlogfile)
%
% Description:
%   This MATLAB function opens and reads data from an MVP .m1 data file
%   and generates the output in the form of a .mat file and simulates the
%   output of the VPR to generate a logfile that can be read into the Video
%   Plankton Recorder vizualization software
%
% Input:
%   An MVP .m1 processed data filename
%   Note: Set VPRindex as a global variable in the Matlab workspace prior
%   to executing this function for the first time.
%
% Output:
%   - MATLAB .mat file
%   - Simulated VPR logfile
%
% Author:
%   Blair Greenan
%   Bedford Institute of Oceanography
%   28-Dec-2011
%   Matlab 7.13
%   blair.greenan@dfo-mpo.gc.ca
%
% Revision History:
%   11 Jan 2011 - modified program to also retrieve the ship's heading as
%                 provided in the NMEA string in the header
%   11 Jan 2011 - previously just reading in the downcast, but now reading
%                 both the downcast and upcast
%
% =========================================================================
%
% 
% Generate an index for each data line (simply to meet the needs of the 
% VPR logfile)
global VPRindex
%
% Change to the working directory
cd 'S:\mvp\nbp1201\rois\vprlog\'
%
% Open the MVP .m1 file which is a processed data file generated by the
% MVP data logging software
fid = fopen(fname);
% Bottom Depth from ship's sounder. Units - meters
BottomDepth = cell2mat(textscan(fid, 'Bottom Depth (m): %f', 'HeaderLines', 6));
% 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');
if strcmpi(LatString{2},'S')
    MVPLatString_downcast = ['-' LatString{1}(1:8)];
else
    MVPLatString_downcast = LatString{1}(1:8);
end
% 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));
% Record whether latitude is North or South
LatNS = LatString{2};
if strcmpi(LatNS,'S')
    LatDeg_downcast = -1.*LatDeg_downcast;
end
clear LatString
% Longitude at start of MVP cast - read as a string and then convert
LonString = textscan(fid, 'ON (dddmm.mmmmmmm,E): %13c,%c');
if strcmpi(LonString{2},'W')
    MVPLonString_downcast = ['-' LonString{1}(1:9)];
else
    MVPLonString_downcast = LonString{1}(1:9);
end
% 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));
% Record whether longitude is East or West
LonEW = LonString{2};
if strcmpi(LonEW,'W')
    LonDeg_downcast = -1.*LonDeg_downcast;
end
clear LonString
% Time at start of MVP cast - UTC
TimeArray = textscan(fid, 'Time (hh|mm|ss.s): %2d:%2d:%f');
THour = TimeArray{1};
TMinute = TimeArray{2};
TSecond = TimeArray{3};
clear TimeArray
% Date at start of MVP cast - UTC
DateArray = textscan(fid, 'Date (dd/mm/yyyy): %2d/%2d/%4d');
DDay = DateArray{1};
DMonth = DateArray{2};
DYear = DateArray{3};
clear DateArray
% Extract the ship's heading from the GPS NMEA string
ShipsHeading = textscan(fid,'BOT F10:$GPVTG,%f', 'HeaderLines', 20);
% Based on the header information, estimate a lat/lon location for the upcast
% assuming that the upcast occurs 6 minutes after the downcast
% Estimate the distance traveled (km) in 6 minutes based on the the ship
% speed
% The ShipsLog is divided by 10 since we want 6/60 minutes
% 1.852 is the conversion factor from nautical miles to kilometers
xkm = (ShipsLog/10) * 1.852;
% 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 6 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)));
delta_lon = delta_lon/(cos(Lat_downcast*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
Lat_upcast = (sign(LatDeg_downcast)*(abs(LatDeg_downcast) + (LatMin_downcast/60))) + delta_lat;
LatDeg_upcast = fix(Lat_upcast);
LatDeg_upcast_str = num2str(LatDeg_upcast);
if (abs(LatDeg_upcast) < 10)
    if (LatDeg_upcast < 0)
        LatDeg_upcast_str(1) = '0';
        LatDeg_upcast_str = ['-' LatDeg_upcast_str];
    else
        LatDeg_upcast_str = ['0' LatDeg_upcast_str];
    end
end
LatMin_upcast = abs(Lat_upcast - LatDeg_upcast)*60;
LatMin_upcast_str = num2str(LatMin_upcast);
if (LatMin_upcast < 10)
    LatMin_upcast_str = ['0' LatMin_upcast_str];
end
% if this string is less than 6 digits long then you need to pad it
for ii = 1:(6-length(LatMin_upcast_str))
    LatMin_upcast_str = [LatMin_upcast_str '0'];
end
Lon_upcast = (sign(LonDeg_downcast)*(abs(LonDeg_downcast) + (LonMin_downcast/60))) + delta_lon;
LonDeg_upcast = fix(Lon_upcast);
LonDeg_upcast_str = num2str(LonDeg_upcast);
if (abs(LonDeg_upcast) < 100)
    if (abs(LonDeg_upcast) < 10)
        if (LonDeg_upcast < 0)
            LonDeg_upcast_str(1) = '0';
            LonDeg_upcast_str = ['-0' LatDeg_upcast_str];
        else
            LonDeg_upcast_str = ['00' LatDeg_upcast_str];
        end
    else
        if (LonDeg_upcast < 0)
            LonDeg_upcast_str(1) = '0';
            LonDeg_upcast_str = ['-' LatDeg_upcast_str];
        else
            LonDeg_upcast_str = ['0' LatDeg_upcast_str];
        end
    end
end
LonMin_upcast = abs(Lon_upcast - LonDeg_upcast)*60;
LonMin_upcast_str = num2str(LonMin_upcast);
if (LonMin_upcast < 10)
    LonMin_upcast_str = ['0' LonMin_upcast_str];
end
% if this string is less than 6 digits long then you need to pad it
for ii = 1:(6-length(LonMin_upcast_str))
    LonMin_upcast_str = [LonMin_upcast_str '0'];
end
% Generate a latitude and longitude string for the upcast
MVPLatString_upcast = [LatDeg_upcast_str,LatMin_upcast_str(1:6)];
MVPLonString_upcast = [LonDeg_upcast_str,LonMin_upcast_str(1:6)];
clear LatDeg_upcast_str LatMin_upcast_str LonDeg_upcast_str LonMin_upcast_str
% Skip the rest of the header and load the data
DataMatrix = textscan(fid,'%f,%f,%f,%f,%f,%f,%f,%f,%f', 'HeaderLines', 16);
% Find the maximum depth and use this to divide the array between up cast
% and downcast
findmax = find(DataMatrix{2} == max(DataMatrix{2}));
% Pressure - Units dbar
Pressure_downcast = DataMatrix{1}(1:findmax);
Pressure_upcast = DataMatrix{1}(end:-1:findmax+1);
% Depth - Units meters
Depth_downcast = DataMatrix{2}(1:findmax);
Depth_upcast = DataMatrix{2}(end:-1:findmax+1);
% Sound Velocity - Units meters per second
SoundVelocity_downcast = DataMatrix{3}(1:findmax);
SoundVelocity_upcast = DataMatrix{3}(end:-1:findmax+1);
% Temperature - Units degrees Celsius
Temperature_downcast = DataMatrix{4}(1:findmax);
Temperature_upcast = DataMatrix{4}(end:-1:findmax+1);
% Conductivity - Units mS/cm
Conductivity_downcast = DataMatrix{5}(1:findmax);
Conductivity_upcast = DataMatrix{5}(end:-1:findmax+1);
% Salinity - PSS-78
Salinity_downcast = DataMatrix{6}(1:findmax);
Salinity_upcast = DataMatrix{6}(end:-1:findmax+1);
% Density - Units kg per cubic meter
Density_downcast = DataMatrix{7}(1:findmax);
Density_upcast = DataMatrix{7}(end:-1:findmax+1);
% Analog Voltage - Units Volt (from ECO fluorometer)
AnalogVoltage_downcast = (DataMatrix{8}(1:findmax))./1000;
AnalogVoltage_upcast = (DataMatrix{8}(end:-1:findmax+1))./1000;
% Convert analog voltage to chlorophyll based on WET Labs calibration
% Analog Range 1
Chl_downcast = 6.*(AnalogVoltage_downcast - 0.065);
Chl_upcast = 6.*(AnalogVoltage_upcast - 0.065);
% Analog Range 2
%Chl = 13.*(AnalogVoltage - 0.035);
% Analog Range 4
%Chl = 25.*(AnalogVoltage - 0.021);
% LOPC - Units counts per second
LOPC_downcast = DataMatrix{9}(1:findmax);
LOPC_upcast = DataMatrix{9}(end:-1:findmax+1);
clear DataMatrix findmax
fclose(fid);
completed = 1;
%
% Convert the data set to the Video Plankton Recorder log file format
fid2 = fopen(logfile,'a+');
% Create a new time array
StartTime = datenum(double(DYear),double(DMonth),double(DDay),double(THour),double(TMinute),double(TSecond));
MVPTime_downcast(1) = StartTime;
for ii = 2:length(Depth_downcast)
    % Increment the MVP time for each depth recorded during free fall
    % Increment each sample by 1 second becuase the VPR software cannot
    % handle lines with the same time
    % Convert from seconds to fractional day
    MVPTime_downcast(ii) = MVPTime_downcast(ii-1) + (1/86400); 
end
%
% Convert MVPTime to Matlab date array
MVPTime_downcast_date_array = datevec(MVPTime_downcast);
MVPHour_downcast = MVPTime_downcast_date_array(:,4);
MVPMin_downcast = MVPTime_downcast_date_array(:,5);
MVPSec_downcast = floor(MVPTime_downcast_date_array(:,6));
MVPDay_downcast = MVPTime_downcast_date_array(:,3);
MVPMonth_downcast = MVPTime_downcast_date_array(:,2);
MVPYear_downcast = MVPTime_downcast_date_array(:,1)-2000;

MVPTime_upcast(1) = MVPTime_downcast(end) + 1/86400;
for ii = 2:length(Depth_upcast)
    % Increment the MVP time for each depth recorded during free fall
    % Increment each sample by 1 second becuase the VPR software cannot
    % handle lines with the same time
    % Convert from seconds to fractional day
    MVPTime_upcast(ii) = MVPTime_upcast(ii-1) + (1/86400); 
end
%
% Convert MVPTime to Matlab date array
MVPTime_upcast_date_array = datevec(MVPTime_upcast);
MVPHour_upcast = MVPTime_upcast_date_array(:,4);
MVPMin_upcast = MVPTime_upcast_date_array(:,5);
MVPSec_upcast = floor(MVPTime_upcast_date_array(:,6));
MVPDay_upcast = MVPTime_upcast_date_array(:,3);
MVPMonth_upcast = MVPTime_upcast_date_array(:,2);
MVPYear_upcast = MVPTime_upcast_date_array(:,1)-2000;
% 
% Write the MVP data to a file in a Video Plankton Recorder logfile format
% Note: Some of the columns are left blank because there is no MVP data
% analagous to the VPR (e.g. pitch, roll, etc.)
for jj = 1:length(Depth_downcast)
    VPRindex = VPRindex + 1;
    fprintf(fid2,'%d %02d%02d%02d %02d%02d%02d %s %s %5.1f 0.0 0.0 0.0 0.0 %5.2f 0.0 %8d %6.3f %6.3f %6.2f \r\n', ...
        VPRindex, MVPHour_downcast(jj), MVPMin_downcast(jj), MVPSec_downcast(jj), ...
        MVPDay_downcast(jj), MVPMonth_downcast(jj), MVPYear_downcast(jj), ... 
        MVPLatString_downcast, MVPLonString_downcast, -1.*BottomDepth, ...
        Chl_downcast(jj), LOPC_downcast(jj), Salinity_downcast(jj), ...
        Temperature_downcast(jj), Depth_downcast(jj));
end

for jj = 1:length(Depth_upcast)
    VPRindex = VPRindex + 1;
    fprintf(fid2,'%d %02d%02d%02d %02d%02d%02d %s %s %5.1f 0.0 0.0 0.0 0.0 %5.2f 0.0 %8d %6.3f %6.3f %6.2f \r\n', ...
        VPRindex, MVPHour_upcast(jj), MVPMin_upcast(jj), MVPSec_upcast(jj), ...
        MVPDay_upcast(jj), MVPMonth_upcast(jj), MVPYear_upcast(jj), ... 
        MVPLatString_upcast, MVPLonString_upcast, -1.*BottomDepth, ...
        Chl_upcast(jj), LOPC_upcast(jj), Salinity_upcast(jj), ...
        Temperature_upcast(jj), Depth_upcast(jj));
end

%Save results to a .mat file
[token, remain] = strtok(fname,'.');
matfile = [token '.mat'];
save(matfile)

fclose(fid2);


