function [theResult]  = read_ndc(varargin)
% read_ndc: Read SEAS XBT files in NODC ascii exchange format (seas2000)
% 
% DESCRIPTION: Paul Chinn generates the so-called nodc ascii format from
% the binary messages received from seas 2000 ships.  This format is also sent
% to AOML PIs (Molly) to archive high density xbt profiles.  This program will read
% the profile and metadata in one ndc file and return a structure.
% 
% INPUT:
%  ndcFile: character array containing path and name of the ndc file.
% 
% OUTPUT:
%  xbt: structure array (length 1) containing the xbt profile information with
%          metadata.  
%
% USAGE:
% >> xbt = read_ndc('mydata/ax10011.ndc')
%       will read the file and return one structure.
% >> read_ndc 
% >> xbt = read_ndc 
%       will print or return an empty data structure for inspection.  Helpful if you don't
%  know the form of the structure that will be returned.
% 
% SEE ALSO:
%  read_seas4x.m
% 
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AUTHOR: Derrick Snowden
%         NOAA/AOML/PhOD
%         Tue Sep 21 13:14:34 EDT 2004
%

% Some defaults
%
% missing data defaults to -9, since that is what is written in a majority of the XBT files
%   saved in "AOML CTD format".  Note that there are some 99.999 default values in some old 
%   STACS ctd files, but this is not a good choice since for example oxygen values may be 
%   within this range (using umol/kg units).
%
error(nargchk(1,2,nargin))
if nargin == 2
	theFileName = varargin{1};
	missingData = varargin{2};
elseif nargin == 1
	theFileName = varargin{1};
	missingData = -9.;
end

% Open the file and test the result.
%    Use read 'r' permissions for text 't' file input
fid = fopen(theFileName,'rt');
if (fid < 0) % make sure file opened correctly
	warning(['Unable to open the file ',theFileName])
	theResult = [];
	return
end

% Parse the data first and then work with the meta data
fullres = get_fullres(fid,theFileName);
twometer = get_twometer(fid,theFileName);
infpts = get_infpts(fid,theFileName);

%
% Rewind and
% Save all header lines in a cell array "UserData"
% The header lines are assumed to be the *first* N nonblank lines.
% 
frewind(fid);
nHdrLines = count_header_lines(fid);
UserData = parse_metadata(fid,max(nHdrLines,12));

% Initialize output structure with
xbt = struct('UserData',[]);
xbt.UserData = UserData;

% Work on the individual lines in UserData.
% Each section is prefaced by a line from  a sample of the header 
% information I am using as a model.
theHdr = UserData;

% SEAS Version:  6.12
idx = strmatch('SEAS Version',theHdr);
if isempty(idx)
	xbt.seas_version = 'Unknown';
else
	xbt.seas_version = sscanf(theHdr{idx},'SEAS Version:%f');
end
% Ship Name: Horizon Hawaii
idx = strmatch('Ship Name',theHdr);
if isempty(idx)
	xbt.ship_name = 'Unknown';
else
	[first,remain] = strtok(theHdr{idx},'Ship Name');
	if isempty(remain)
		xbt.ship_name = 'Unknown';
	else
		xbt.ship_name = fliplr(deblank(fliplr(deblank(remain))));
	end
end

% Call Sign: TEST
idx = strmatch('Call Sign:',theHdr);
if isempty(idx)
	xbt.call_sign = 'Unknown';
else
	[first,remain] = strtok(theHdr{idx},'Call Sign');
	if isempty(remain)
		xbt.call_sign = 'Unknown';
	else
		xbt.call_sign = fliplr(deblank(fliplr(deblank(remain))));
	end
end

% Lloyds Number:  1234567
idx = strmatch('Lloyds Number:',theHdr);
if isempty(idx)
	xbt.lloyds_number = 99999;
else
	[first,remain] = strtok(theHdr{idx},'Lloyds Number');
	if isempty(remain)
		xbt.lloyds_number = 99999;
	else
		xbt.lloyds_number = fliplr(deblank(fliplr(deblank(remain))));
	end
end

% Date/Time(dd/mm/yyyy): 05/06/2004 20:30 GMT
match = 'Date/Time(dd/mm/yyyy):';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.obsDateTime = 'Unknown';
	xbt.obsTimeZone = 'Unknown';
else
	parts = strsplit(match,theHdr{idx});
	parts = sscanf(parts{2},'%d/%d/%d %d:%d %s');
	yr = parts(3); mnth = parts(2); day = parts(1); hr = parts(4); mins = parts(5);
	timezone = char(parts(6:end))';
	xbt.obsDateTime = datestr ([yr mnth day hr mins 0],0);
	xbt.obsTimeZone = timezone;
end

% Latitude(ddd.ddd): 036.888 N
match = 'Latitude(ddd.ddd):';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.lat = nan;
else
	parts = strsplit(match,theHdr{idx});
	parts = sscanf(parts{2},'%f %s');
	hemi = char(parts(2));
	if isequal(upper(hemi),'N')
		sgn = 1;
	else
		% TODO: % Make this logic more robust for non N or S strings.
		sgn = -1;
	end
	xbt.lat = sgn.*parts(1);
end

% Longitude(ddd.ddd): 0072.717 W
match = 'Longitude(ddd.ddd):';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.lon = nan;
else
	parts = strsplit(match,theHdr{idx});
	parts = sscanf(parts{2},'%f %s');
	hemi = char(parts(2));
	if isequal(upper(hemi),'E')
		sgn = 1;
	else
		% TODO: % Make this logic more robust for non N or S strings.
		sgn = -1;
	end
	xbt.lon = sgn.*parts(1);
end

% Probe Type: Sippican Deep Blue
match = 'Probe Type:';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.probe_type = 'Unknown';
else
	parts = strsplit(match,theHdr{idx});
	xbt.probe_type = ddeblank(parts{2});
end

% Probe Code:  52
match = 'Probe Code:';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.probe_code = 'Unknown';
else
	parts = strsplit(match,theHdr{idx});
	xbt.probe_code = str2num(ddeblank(parts{2}));
end

% Probe Serial No: 0000000
match = 'Probe Serial No:';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.probe_serial_num = 'Unknown';
else
	parts = strsplit(match,theHdr{idx});
	xbt.probe_serial_num = ddeblank(parts{2});
end

% Recorder Type: Sippican MK-21
match = 'Recorder Type:';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.recorder_type = 'Unknown';
else
	parts = strsplit(match,theHdr{idx});
	xbt.recorder_type = ddeblank(parts{2});
end

% Recorder Code:  6
match = 'Recorder Code:';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.recorder_code = 'Unknown';
else
	parts = strsplit(match,theHdr{idx});
	xbt.recorder_code = ddeblank(parts{2});
end

% Bottom Depth: 9999 M
match = 'Bottom Depth:';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.ocean_bottom_depth = 'Unknown';
else
	parts = strsplit(match,theHdr{idx});
	xbt.ocean_bottom_depth = ddeblank(parts{2});
end

% GTS CRC: 27A0C3DF
match = 'GTS CRC:';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.gts_crc = 'Unknown';
else
	parts = strsplit(match,theHdr{idx});
	xbt.gts_crc = ddeblank(parts{2});
end

% SEAS ID: FEE2A519
match = 'SEAS ID:';
idx = strmatch(match,theHdr);
if isempty(idx)
	xbt.seas_id = 'Unknown';
else
	parts = strsplit(match,theHdr{idx});
	xbt.seas_id= ddeblank(parts{2});
end

%
% Finished parsing the meta data
% Now it remains to interpret the xbt.probe_code and apply the 
% fall rate equation
% This interpretation of the xbt.probe_code was given to me by Paul Chinn in Sep 2004
% It corresponds to the WMO Code Table
% Code Table 1770: http://www.nodc.noaa.gov/GTSPP/document/codetbls/wmocode.htm
% Ix Ix Ix - Instrument type for temperature, with fall rate equation coefficients for XBT 
% 
% The notation used to refer to the WMO code table is slightly different between 
% SEAS2000 (v6.0 or>) and SEAS version 4x or lower files.  
% In SEAS4x there are "Code" *and* "Eqn" fields immediately 
% following the "Probe Type" field.
% In SEAS200 a Sippican Deep Blue XBT with the updated Hanawa 1995 fall rate equation 
% applied to the full resolution data would appear with the following lines in the header
% section of the file.
% Probe Type: Sippican Deep Blue
% Probe Code: 52
%
% whereas in a SEAS v4x < file the following line would appear
% SEAS Version 4.00  Drop 035  Probe Type Deep Blue           Code 05 Equation 2
% The are to be interpreted identically.
xbt.infpts = infpts;
xbt.twometer = twometer;
xbt.te = fullres;
fallrate = wmocode1770(xbt.probe_code);
dpth = apply_fall_rate(fallrate.a,fallrate.b,length(fullres));
xbt.de = dpth(:);
xbt.castNum = NaN;
xbt.dataNames = {'te' 'de'};
xbt.dataUnits = {'degC' 'm'};

xbt.numDpthLevels = length(xbt.te);

if nargout > 0
        theResult = xbt;
else
        assignin('caller', 'ans', xbt)
end
return


return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fullres]  = get_fullres(fid, theFileName)
% get_fullres: Parse the full resolution data.
%
% Determine what line the XBT time series data begins on:
%     Search for string 'XBT' as the delimiter.
%     XBT temperature time series begins on the next line
%
frewind(fid);
cnt = 1;
file_line = fgetl (fid);
while ~feof(fid) & (isempty(file_line) | ~strcmp(file_line(1:3), 'XBT'))
     file_line = fgetl (fid);
     cnt = cnt + 1;
end

% Record the number of points expected in the full resolution trace.
nFullres = sscanf(file_line,'XBT: %f\n');
startFullres = cnt + 1;

% Should be at the start of the full res data stream
temp = fscanf(fid,'%4d');

if isempty(temp)
	fullres = [];
	lasterr('') % This may not be an error.  Reset lasterr and throw a warning.
	warning('XBT:NoFullResData',['No full resolution data available in ' theFileName])
	return
elseif length(temp) ~= nFullres
	fullres = [];
	warning('XBT:FullResDataCountMismatch',...
	['Expected ' num2str(nFullres) ' full resolution points and ' num2str(length(temp)) ... 
	' were read from ' theFileName]);
	return
else
	fullres = temp(:)./100; % Data stored in hundredths of degrees C;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xbt]  = get_twometer(fid, theFileName)
% get_twometer: Parse the two meter resolution data.
%
% Determine what line the XBT time series data begins on:
%     Search for string 'XBT' as the delimiter.
%     XBT temperature time series begins on the next line
%
frewind(fid);
cnt = 1;
file_line = fgetl (fid);
while ~feof(fid) & (isempty(file_line) | ~strcmp(file_line(1:8), 'TWOMETER'))
     file_line = fgetl (fid);
     cnt = cnt + 1;
end

% Record the number of points expected in the full resolution trace.
nExpected = sscanf(file_line,'TWOMETER: %f\n');

% Should be at the start of the full res data stream
temp = fscanf(fid,'%4d');
nRead = length(temp)./2;
dpth = temp(1:2:end-1)./10;
tmp = temp(2:2:end)./100;

% Output logic.
if isempty(temp)
	xbt = [];
	lasterr('') % This may not be an error.  Reset lasterr and throw a warning.
	warning('XBT:NoTwoMeterData',['No two meter resolution data available in ' theFileName])
	return
elseif nRead ~= nExpected
	xbt = [];
	warning('XBT:TwoMeterDataCountMismatch',...
	['Expected ' num2str(nExpected) ' two meter resolution data points and ' num2str(nRead) ' were read from ' theFileName]);
	return
else
	xbt.de = dpth;
	xbt.te = tmp;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xbt]  = get_infpts(fid, theFileName)
% get_twometer: Parse the two meter resolution data.
%
% Determine what line the XBT time series data begins on:
%     Search for string 'XBT' as the delimiter.
%     XBT temperature time series begins on the next line
%
frewind(fid);
cnt = 1;
file_line = fgetl (fid);
while ~feof(fid) & (isempty(file_line) | ~strcmp(file_line(1:6), 'INFPTS'))
     file_line = fgetl (fid);
     cnt = cnt + 1;
end

% Record the number of points expected in the full resolution trace.
nExpected = sscanf(file_line,'INFPTS:: %f\n');

% Should be at the start of the full res data stream
temp = fscanf(fid,'%4d');
nRead = length(temp)./2;
dpth = temp(1:2:end-1)./10;
tmp = temp(2:2:end)./100;

% Output logic.
if isempty(temp)
	xbt = [];
	lasterr('') % This may not be an error.  Reset lasterr and throw a warning.
	warning('XBT:NoInflectionPointData',['No inflection point data available in ' theFileName])
	return
elseif nRead ~= nExpected
	xbt = [];
	warning('XBT:InflectionPointDataCountMismatch',...
	['Expected ' num2str(nExpected) ' inflection data points and ' num2str(nRead) ' were read from ' theFileName]);
	return
else
	xbt.dpth = dpth;
	xbt.temperature = tmp;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function nHdrLines = count_header_lines(fid);
% Count the number of non-blank lines starting at the file beginning
%
file_line = 'begin';
frewind(fid);
cnt=0;
while ~feof(fid) & ~isempty(file_line)
     file_line = fgetl (fid);
     cnt = cnt + 1;
end
nHdrLines = cnt;
frewind(fid);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function UserData = parse_metadata(fid,nLines);
% Create a cell array and store meta-data for safekeeping.
%
% This routine only puts the first nLines header lines into
% UserData{1:nLines}
%
UserData = cell(nLines,1);
frewind(fid);
for i = 1:nLines
	UserData{i} = fliplr(deblank(fliplr(deblank(fgetl(fid)))));
end
return


