%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Read .NDC files by (1) read_ndc.m, or (2) read_seas4x.m %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% (1)
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

%
% Read in temperature time series until EOF reached.
%
file_line = fgetl (fid);
temp = [];
while ~feof(fid) 
     temp = [temp ; str2num(reshape(file_line, 4, length(file_line)/4)')/100. ];
     file_line = fgetl (fid);
end

fullres = temp;

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% (2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xbt  = read_seas4x (theFileName,varargin)
% read_aoml: Read an ascii file of XBT data in the so-called SEAS/NODC format.
%
% USAGE:  dataStruct = read_seas4x (thePathAndFilename, [fillValue])
% 
% INPUT:
%          thePathAndFilename: string containing the pathand filename to be read
%          fillValue: optional value representing fill value.  If passed, fillValue will be 
%                         replaced by NaNs in the output.
% 
% OUTPUT:  
%          dataStruct:  Structure containing the data in the input file.  See examples below.
% 
% DESCRIPTION: 
%     IF YOU FIND ANY BUGS, please do not fix them on your own without at least telling me.
%     I'm trying to create a general PhOD toolbox that is available to everyone in a central 
%     location (/phodnet/share/matlab or on the intranet) Therefore I want these programs to be
%     as general and well tested as possible. If you fix a bug please contact me.
% 
% The SEAS/NODC ascii format is used for files containing 
% XBT profile data.  Typically, these files are output from the XBT autolaunching software.
% Originally, this was the ascii format produced by version 4 of the SEAS software used 
% for hand launchers.  Several "versions" of the SEAS software changed slightly the 
% format of the date string in particular.
% The formatting rules implemented in this program are:
% 1. Each file contains one and only one profile.
% 2. The date/time/postion is found on the third line.  They are read using
%     c style formatted read statements which means, I don't count columns, but the presence 
%     of the slashes in the date are expected.  White space, no matter the size, is treated as 
%     one space.  
% 3. The other header lines are typically not used and for this program will be discarded.
%    Hopefully, Derrick or someone very good at this sort of thing will change that.  
%    For now, the header lines will just be read in as character strings.
% 4. There is a section that begins with the string "INFPTS" which lists 
%    depth/tempearture pairs produced as inflection points to the full resolution profile.  
%    Typically, these values are the ones sent over GOES/GTS in real-time.  For now this 
%    program will ignore these values as well.
% 5. The XBT temperature time series follows the section preceeded by the string "XBT".  
%    Temperature data is stored in a 4-digit format with 2 decimal places (no ".")  
%    e.g. 2037 = 20.37 degrees C.
%    There are 20 temperature values on each line (80 columns wide)
% 6. Time is returned in a string of the form 
%     obsTime = 'dd-mmm-yyyy HH:MM:SS'   e.g. 01-Mar-2000 15:45:17
%     which is dateform #0 in the Matlab datestr  documentation.  
%     See below for help in dealing with this form.
%   
%     e.g. dtNum = datenum(obsTime,0) or [yr,mnth,day,hh,mm,ss] = datevec(dtNum);
% 
% 7. obsTimeZone contains the time zone information in column 4 of line 3.
% 8. lat and lon are assumed to be in degrees and decimal-minute format with lon in 
%    the western hemisphere marked as "W", latitudes in the southern hemisphere 
%    marked as "S" on line 3.
% 
%
% Depth is determined from the fall rate equation assuming 
%    iprobe_type = 3 = Deep Blue = T7s
% In the future you could certainly read in the probe type number.
% However, I've found that the probe number in the file is often incorrect.
% I would recommend reading in the probe number if and only if (IFF)
% the probe number is not an optional input variable.
%
% Output field UserData contains the ascii lines found in the header of the SEAS4x file.
% UserData{1:2} empty
% UserData{3} filename of original SEAS format file
% UserData{4:10} the first 7 lines in the original Seas4X file
% UserData{11:12} empty
%
% ROUTINES CALLED:
%   sw_pres - slightly improved version of sw_prs from CSIRO oceans toolbox
%   determine_depth_xbt - included below.
%
% AUTHOR: Molly Baringer
%         NOAA/AOML/PhOD
% 
% SEE ALSO:
%    read_aoml.m written by Derrick Snowden - which this program is directly modelled.
%    fortran routines read_seas4_xbt and determine_depth_xbt that were taken from 
%         Molly's hdenxbt page   http://www.aoml.noaa.gov/phod/hdenxbt/README
%    
% DEPENDENCIES:
%
% EXAMPLES:
% aoml = read_seas4x('/mobster/m1/data/xbt/seas/ax7/apr04/ax70404.181')

% 
% REFERENCES:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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).
%
missingData = -9;
if nargin == 2
	missingData = varargin{1};
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])
	outProf = [];
	return
end

%
% 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
%
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

%
% Read in temperature time series until EOF reached.
%
file_line = fgetl (fid);
temp = [];
while ~feof(fid) 
     temp = [temp ; str2num(reshape(file_line, 4, length(file_line)/4)')/100. ];
     file_line = fgetl (fid);
end

%
% Rewind and
% Save up to 12 "UserData" lines of the original file
%
% Save all meta-data lines in UserData
UserData = parse_metadata(fid,min(cnt,12));

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

%
% Note that Date/Time for Seas4X occurs on line 3 and looks like:
%  Date/Time 01/05/2004 22:16 GMT  Latitude  26 21.5 N     Longitude 073 06.5 W
% While the Date/Time for Seas 2000 *>NDC files is on line 5 and looks like:
%  Date/Time(dd/mm/yyyy): 02/27/2004 17:03 GMT
% Hence a generic processing step would be to search the file for "Date" at the
% start of each line and then read 
%
% The steps below work on both file types.
%
frewind(fid);
file_line = fgetl (fid);
while ~feof(fid) & (isempty(file_line) | ~strcmp(file_line(2:5), 'Date'))
     file_line = fgetl (fid);
end
%
% Get rid of special characters characters
% Some examples to test:
% file_line = ' Date/Time(dd/mm/yyyy): 27/02/2004 17:03 GMT'
% file_line = '  Date/Time 01/05/2004 22:16 GMT  Latitude  26 21.5 N     Longitude 073 06.5 W'
file_line=strrep (file_line, ':', ' ');
file_line=strrep (file_line, '(dd/mm/yyyy)', '            ');
file_line=strrep (file_line, '/', ' ');
%
% Break line up into discrete units
%
i = 1;
while file_line
	[colNames{i},file_line] = strtok(file_line,' ');
	i = i + 1;
end
[year month day junk junk junk] = datevec([colNames{4} '/' colNames{3} '/' colNames{5}]);
[junk junk junk  hr  minute  sec]  = datevec([colNames{6} ':' colNames{7}]);
xbt.obsDateTime = datestr ([year month day hr minute sec],0);
xbt.obsTimeZone = colNames{8};

%
% Note: now to get lat/lon/station number from header information
% will depend on the exact format of the files.  Up to this point
% this routine can be used as is to read SEAS2000 *.NDC format 
% files (e.g. the ascii files that Paul produces that contain inflection 
% points and two-meter data. 
%

%
% Parse line 1 for XBT number
%   Note SEAS 2000 *NDC files have no station number in them!
%   Will need to make one up.
%
file_line = UserData{4};
%
% Break line up into discrete units
%
i = 1;
while file_line
 	[colNames{i},file_line] = strtok(file_line,' ');
	i = i + 1;
end
%
% Find value just after specified column name
%
for I = 1: length(colNames)
     if strcmp (colNames{I}, 'Drop')
         xbt.castNum = str2num(colNames{I+1});
     end
end

% 
% Parse line 3 for lat/lon/date/time etc.
%
file_line = UserData{6};
%
% Get rid of ":" characters
%
file_line=strrep (file_line, ':', ' ');
%
% Break line up into discrete units
%
i = 1;
while file_line
	[colNames{i},file_line] = strtok(file_line,' /');
	i = i + 1;
end

xbt.lat  = str2num(colNames{10}) + str2num(colNames{11})/60.;
xbt.lon = str2num(colNames{14}) + str2num(colNames{15})/60.;

%
% Correct the sign of the latitude/longitude
%
if colNames{12} == 'S' 
     xbt.lat = -1.*xbt.lat;
end
if colNames{16} == 'W' 
     xbt.lon = -1.*xbt.lon;
end

xbt.numDpthLevels = length(temp);

xbt.te = temp(:)';

%
% Determine Depth from fall rate equation
%
% where iprobe_type = 3 = Deep Blue = T7s
% In the future you could certainly read in the probe type number.
% However, I've found that the probe number in the file is often incorrect.
% I would recommend reading in the probe number if and only if (IFF)
% the probe number is not an optional input variable.
%
xbt.de = determine_depth_xbt(length(temp), 3);

xbt.pr = sw_pres (xbt.de', xbt.lat)';

xbt.dataUnits = {'degC' 'm' 'dbar'};
xbt.dataNames = {'te' 'de' 'pr'};

xbt.UserData{3} = ['Wrote Seas4x file ' theFileName ' ' date];  

return

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dpth]  = determine_depth_xbt(npts, iprobe_type)
%
% determine_depth_xbt: Apply fall rate equation to time series of samples.
%
% Usage:  [dpth]  = determine_depth_xbt(npts, iprobe_type)
%
%  Based on the following subroutine
% c**********************************************************************
%     subroutine determine_depth_xbt (p, npts, iprobe_type)
%     real p(1)
%     real aconst(8), bconst(8)
%
%     data aconst / 6.828,   6.346,   6.691,   6.3301,   1.779,  
%     +              6.301,   5.861,   6.472   /
%     data bconst / 0.00182, 0.00182, 0.00225, 0.00216, 0.000255, 
%     +              0.00216,.0000904, 0.00216 /
% c
% c This routine determines the depth of an XBT probe based on its fall 
% c rate.  
% c   aconst, bconst are the coefficents for the fall rate equations for
% c   various probe types.  
% c
% c iprobe_type = probe type number
% c            = 1 = Sippican T-5, (XBT-5, XBT-5DB)
% c            = 2 = Sippican Fast Deep
% c            = 3 = Sippican (T-4, T-7, T-6, Deep Blue), 
% c                  TSK (T-4, T-6, T-7)
% c            = 4 = T-10
% c            = 5 = T-11
% c            = 6 = Spartan (XBT-1, XBT-10)
% c            = 7 = Spartan (XBT-3)
% c            = 8 = Spartan (XBT-4, XBT-6, XBT-7, XBT-7DB, XBT-20)
% c                           (...XBT-20DB)
% c Note:
% c   These values obtained from http://dbcp.nos.noaa.gov/seas/jjyy.txt
% c on Jan 29, 1998.  There are some errors with the numbers quoted 
% c on this page, namely bconst is off by a factor of 100.
% c-
%     iprobe_type = 3
%     do j = 1, npts
%       p(j) = aconst(iprobe_type) * 0.1*j - 
%     +                  bconst(iprobe_type) * 0.01*j*j
%     end do
%     return
%     end
% c end determine_depth_xbt (p, npts, iprobe_type)

aconst = [6.828,   6.346,   6.691,   6.3301,   1.779,    6.301,   5.861,   6.472  ];
bconst = [0.00182, 0.00182, 0.00225, 0.00216, 0.000255, 0.00216,.0000904, 0.00216 ];
j = 0.1.*[1:npts];
dpth = aconst(iprobe_type).*j - bconst(iprobe_type).*j.*j;

return