carter

 Function to apply Carter table correction to echosounder
 data that was collected assuming sound speed = 1500 m/s.

 Inputs are (latitude, longitude, uncorrected depth)
 Inputs should have South and East = -ve
 Inputs can be row or column arrays, but not matrices (function is
 designed for single-beam echosounder data)

 Uses lookup tables BOUNDARY2.DAT and CORRECTN2.DAT

 Outputs are corrected depth, and Carter correction area
 (Note that function returns correction_area = NaN if no correction
 is applied because input depth is <200m)

 Author: Mike Meredith, British Antarctic Survey, Dec 2005

 Stuck in ice onboard RRS James Clark Ross, off Adelaide Island, 
   Antarctic Peninsula

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [cordepth,correction_area] = carter(lat,lon,uncdepth);
0002 % Function to apply Carter table correction to echosounder
0003 % data that was collected assuming sound speed = 1500 m/s.
0004 %
0005 % Inputs are (latitude, longitude, uncorrected depth)
0006 % Inputs should have South and East = -ve
0007 % Inputs can be row or column arrays, but not matrices (function is
0008 % designed for single-beam echosounder data)
0009 %
0010 % Uses lookup tables BOUNDARY2.DAT and CORRECTN2.DAT
0011 %
0012 % Outputs are corrected depth, and Carter correction area
0013 % (Note that function returns correction_area = NaN if no correction
0014 % is applied because input depth is <200m)
0015 %
0016 % Author: Mike Meredith, British Antarctic Survey, Dec 2005
0017 %
0018 % Stuck in ice onboard RRS James Clark Ross, off Adelaide Island,
0019 %   Antarctic Peninsula
0020 %
0021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0022 
0023 % CHECK INPUTS
0024 
0025 if nargin ~= 3
0026    error('Error: Must pass 3 parameters ')
0027 end %if
0028 
0029 % CHECK input dimensions and verify consistent
0030 [mlat,nlat] = size(lat);
0031 [mlon,nlon] = size(lon);
0032 [mdepth,ndepth] = size(uncdepth);
0033   
0034 if (mlat~=1 & nlat~=1) | (mlon~=1 & nlon~=1) | (mdepth~=1 & ndepth~=1)
0035     error('Error: must pass arrays not matrices');
0036 end
0037 
0038 % CHECK THAT INPUTS HAVE SAME LENGTH
0039 slat = length(lat);
0040 slon = length(lon);
0041 sdepth = length(uncdepth);
0042 if (slat~=slon) | (slat~=sdepth) | (slon~=sdepth)
0043    error('Error: inputs must have same lengths')
0044 end %if
0045 
0046 % open files containing area definitions and corrections
0047 fid = fopen('BOUNDARY2.DAT');
0048 fid2 = fopen('CORRECTN2.DAT');
0049 
0050 %% main loop
0051 
0052 for ii = 1:slat;
0053 % Check that lats and lons are reasonable
0054 if lat(ii)<-90 | lat(ii)>90 | lon(ii) < -180 | lon(ii) > 180;
0055     error('Error: Navigation out of range')
0056 end; %if
0057 
0058 % Check that inputs are not NaNs, and depth is in range of tables (will
0059 % return NaN in first case, and will leave depth uncorrected in second)
0060 
0061 if isnan(lat(ii)) ~= 1 & isnan(lon(ii)) ~= 1 & isnan(uncdepth(ii)) ~= 1 & uncdepth(ii) > 200;
0062 
0063 % round down lat and lon
0064 lat_use = floor(lat(ii));
0065 lon_use = floor(lon(ii));
0066 
0067 lat_file = 999;
0068 
0069 % scan through file until find line with appropriate latitude
0070 while lat_file ~= lat_use;
0071 
0072 a = fgetl(fid);
0073 b = str2num(a);
0074 
0075 num_records = b(4);
0076 
0077 if num_records <=9;
0078     a1 = fgetl(fid); b1 = str2num(a1);
0079     lon_file = b1(4:2:end);
0080     area = b1(5:2:end);
0081 elseif num_records > 9 & num_records <= 18;
0082     a1 = fgetl(fid); b1 = str2num(a1);
0083     a2 = fgetl(fid); b2 = str2num(a2);
0084     lon_file = [b1(4:2:end) b2(4:2:end)];
0085     area = [b1(5:2:end) b2(5:2:end)]; 
0086 elseif num_records > 18 & num_records <= 27;
0087     a1 = fgetl(fid); b1 = str2num(a1);
0088     a2 = fgetl(fid); b2 = str2num(a2);
0089     a3 = fgetl(fid); b3 = str2num(a3);
0090     lon_file = [b1(4:2:end) b2(4:2:end) b3(4:2:end)];
0091     area = [b1(5:2:end) b2(5:2:end) b3(5:2:end)]; 
0092 elseif num_records > 27 & num_records <= 36;
0093     a1 = fgetl(fid); b1 = str2num(a1);
0094     a2 = fgetl(fid); b2 = str2num(a2);
0095     a3 = fgetl(fid); b3 = str2num(a3);
0096     a4 = fgetl(fid); b4 = str2num(a4);
0097     lon_file = [b1(4:2:end) b2(4:2:end) b3(4:2:end) b4(4:2:end)];
0098     area = [b1(5:2:end) b2(5:2:end) b3(5:2:end) b4(5:2:end)]; 
0099 else;
0100     a1 = fgetl(fid); b1 = str2num(a1);
0101     a2 = fgetl(fid); b2 = str2num(a2);
0102     a3 = fgetl(fid); b3 = str2num(a3);
0103     a4 = fgetl(fid); b4 = str2num(a4);
0104     a5 = fgetl(fid); b5 = str2num(a5);
0105     lon_file = [b1(4:2:end) b2(4:2:end) b3(4:2:end) b4(4:2:end) b5(4:2:end)];
0106     area = [b1(5:2:end) b2(5:2:end) b3(5:2:end) b4(5:2:end) b5(5:2:end)];
0107 end; % if num_records
0108 
0109 lat_file = b(2);
0110 end; % while
0111 
0112 
0113 % now scan through longitudes to find correction area
0114 
0115 for i = 1:length(lon_file)-1;
0116     if lon_use >= lon_file(i) & lon_use < lon_file(i+1);
0117         correction_area(ii) = area(i);
0118     elseif lon_use > lon_file(end);
0119         correction_area(ii) = area(end); % if is on last record in group
0120     end;
0121 end;
0122 
0123 frewind(fid);
0124 
0125 % now read corrections table to find corrected depth
0126 
0127 area_file = 999;
0128 
0129 % scan through file until find line with appropriate latitude
0130 while area_file ~= correction_area(ii);
0131 
0132 a = fgetl(fid2);
0133 b = str2num(a);
0134 
0135 area_file = b(2);
0136 num_records = b(4);
0137 
0138 if num_records <=12;
0139     a1 = fgetl(fid2); b1 = str2num(a1);
0140     data = b1(4:end);
0141 elseif num_records > 12 & num_records <= 24;
0142     a1 = fgetl(fid2); b1 = str2num(a1);
0143     a2 = fgetl(fid2); b2 = str2num(a2);
0144     data = [b1(4:end) b2(4:end)]; 
0145 elseif num_records > 24 & num_records <= 36;
0146     a1 = fgetl(fid2); b1 = str2num(a1);
0147     a2 = fgetl(fid2); b2 = str2num(a2);
0148     a3 = fgetl(fid2); b3 = str2num(a3);
0149     data = [b1(4:end) b2(4:end) b3(4:end)]; 
0150 elseif num_records > 36 & num_records <= 48;
0151     a1 = fgetl(fid2); b1 = str2num(a1);
0152     a2 = fgetl(fid2); b2 = str2num(a2);
0153     a3 = fgetl(fid2); b3 = str2num(a3);
0154     a4 = fgetl(fid2); b4 = str2num(a4);
0155     data = [b1(4:end) b2(4:end) b3(4:end) b4(4:end)]; 
0156 elseif num_records > 48 & num_records <= 60;
0157     a1 = fgetl(fid2); b1 = str2num(a1);
0158     a2 = fgetl(fid2); b2 = str2num(a2);
0159     a3 = fgetl(fid2); b3 = str2num(a3);
0160     a4 = fgetl(fid2); b4 = str2num(a4);
0161     a5 = fgetl(fid2); b5 = str2num(a5);
0162     data = [b1(4:end) b2(4:end) b3(4:end) b4(4:end) b5(4:end)]; 
0163 elseif num_records > 60 & num_records <= 72;
0164     a1 = fgetl(fid2); b1 = str2num(a1);
0165     a2 = fgetl(fid2); b2 = str2num(a2);
0166     a3 = fgetl(fid2); b3 = str2num(a3);
0167     a4 = fgetl(fid2); b4 = str2num(a4);
0168     a5 = fgetl(fid2); b5 = str2num(a5);
0169     a6 = fgetl(fid2); b6 = str2num(a6);
0170     data = [b1(4:end) b2(4:end) b3(4:end) b4(4:end) b5(4:end) b6(4:end)]; 
0171 elseif num_records > 72 & num_records <= 84;
0172     a1 = fgetl(fid2); b1 = str2num(a1);
0173     a2 = fgetl(fid2); b2 = str2num(a2);
0174     a3 = fgetl(fid2); b3 = str2num(a3);
0175     a4 = fgetl(fid2); b4 = str2num(a4);
0176     a5 = fgetl(fid2); b5 = str2num(a5);
0177     a6 = fgetl(fid2); b6 = str2num(a6);
0178     a7 = fgetl(fid2); b7 = str2num(a7);
0179     data = [b1(4:end) b2(4:end) b3(4:end) b4(4:end) b5(4:end) b6(4:end) b7(4:end)];
0180 elseif num_records > 84 & num_records <= 96;
0181     a1 = fgetl(fid2); b1 = str2num(a1);
0182     a2 = fgetl(fid2); b2 = str2num(a2);
0183     a3 = fgetl(fid2); b3 = str2num(a3);
0184     a4 = fgetl(fid2); b4 = str2num(a4);
0185     a5 = fgetl(fid2); b5 = str2num(a5);
0186     a6 = fgetl(fid2); b6 = str2num(a6);
0187     a7 = fgetl(fid2); b7 = str2num(a7);
0188     a8 = fgetl(fid2); b8 = str2num(a8);
0189     data = [b1(4:end) b2(4:end) b3(4:end) b4(4:end) b5(4:end) b6(4:end) b7(4:end) b8(4:end)];
0190 elseif num_records > 96 & num_records <= 108;
0191     a1 = fgetl(fid2); b1 = str2num(a1);
0192     a2 = fgetl(fid2); b2 = str2num(a2);
0193     a3 = fgetl(fid2); b3 = str2num(a3);
0194     a4 = fgetl(fid2); b4 = str2num(a4);
0195     a5 = fgetl(fid2); b5 = str2num(a5);
0196     a6 = fgetl(fid2); b6 = str2num(a6);
0197     a7 = fgetl(fid2); b7 = str2num(a7);
0198     a8 = fgetl(fid2); b8 = str2num(a8);
0199     a9 = fgetl(fid2); b9 = str2num(a9);
0200     data = [b1(4:end) b2(4:end) b3(4:end) b4(4:end) b5(4:end) b6(4:end) b7(4:end) b8(4:end) b9(4:end)];
0201 else;  % for num_records > 108
0202     a1 = fgetl(fid2); b1 = str2num(a1);
0203     a2 = fgetl(fid2); b2 = str2num(a2);
0204     a3 = fgetl(fid2); b3 = str2num(a3);
0205     a4 = fgetl(fid2); b4 = str2num(a4);
0206     a5 = fgetl(fid2); b5 = str2num(a5);
0207     a6 = fgetl(fid2); b6 = str2num(a6);
0208     a7 = fgetl(fid2); b7 = str2num(a7);
0209     a8 = fgetl(fid2); b8 = str2num(a8);
0210     a9 = fgetl(fid2); b9 = str2num(a9);
0211     a10 = fgetl(fid2); b10 = str2num(a10);
0212     data = [b1(4:end) b2(4:end) b3(4:end) b4(4:end) b5(4:end) b6(4:end) b7(4:end) b8(4:end) b9(4:end) b10(4:end)];
0213 end; % if num_records
0214 
0215 end; % while
0216 
0217 uncdepth_file = 200:100:10000;
0218 uncdepth_file = uncdepth_file(1:length(data));
0219 
0220 for i = 1:length(data)-1;
0221     if uncdepth(ii) >= uncdepth_file(i) & uncdepth(ii) < uncdepth_file(i+1);
0222         cordepth_high = data(i+1);
0223         cordepth_low = data(i);
0224         uncdepth_low = uncdepth_file(i);
0225     end;
0226 end;
0227 
0228 % now linearly interpolate to find corrected depth
0229 
0230 offset = uncdepth(ii) - uncdepth_low;
0231 cordepth(ii) = cordepth_low + ((cordepth_high-cordepth_low)*offset)./100;
0232 
0233 frewind(fid2);
0234 
0235 else;  % (if one of inputs is NaN or depth < 200)
0236     
0237 cordepth(ii) = uncdepth(ii); 
0238 correction_area(ii) = NaN;
0239 
0240 end; %(if one of inputs is NaN or depth < 200)
0241 
0242 end; % ii, main loop
0243 
0244 fclose(fid);
0245 fclose(fid2);
0246 
0247 return;