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