version on whoosher Apr 10. 2021
--- a/HISTORY Wed Jan 17 12:29:40 2018 -0500
+++ b/HISTORY Sat Apr 10 08:03:07 2021 -0400
@@ -1,9 +1,9 @@
======================================================================
H I S T O R Y
doc: Tue Aug 16 11:41:46 2005
- dlm: Wed Jan 17 12:19:26 2018
+ dlm: Wed Sep 4 18:13:56 2019
(c) 2005 A.M. Thurnherr
- uE-Info: 325 17 NIL 0 0 72 3 2 4 NIL ofnI
+ uE-Info: 343 0 NIL 0 0 72 3 2 4 NIL ofnI
======================================================================
NB: CHANGE VERSION IN [default.m] BEFORE UPLOADING
@@ -323,4 +323,29 @@
- updated version to IX_13 in [default.m] [.hg/hgrc]
- commit
- publish
-
+ ...
+ Sep 14, 2018:
+ - updated version to IX_14beta in [default.m] [.hg/hgrc]
+ - fixed (without much checking) bug introduced in 2008 code move
+ from [loadctd.m] to [loadnav.m]
+ ...
+ Feb 8, 2019:
+ - added pauses to figure-saving loop, because on TheThinMint
+ occasionally the wrong figures were saved
+ Feb 16, 2019:
+ - moved post-processing to step 17 to ensure it is done before
+ results are saved
+ ...
+ Sep 4, 2019:
+ - changed p.btrk_default from 3 to 2
+ - fixed [savearch.m] date bug
+ - swapped in Gerd's code [magdev.m] and made default in [loadnav.m]
+ - modifications suggested by Gerd:
+ - [getinv.m] rounded default ps.dz
+ - [prepinv.m] BUG related to weights
+ - [getbtrack.m] fixed minor typo
+ - [checkbtrk.m] fixed plotting error
+ - adapted [process_cast.m] to new [calc_shear3.m] from GK
+ - [calc_shear3.m] required some adaptation
+ - updated version in [default.m] to 15beta, skipping 14, which
+ was used for GO-SHIP work but never published
--- a/IGRF00.m Wed Jan 17 12:29:40 2018 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,91 +0,0 @@
-function [gh,G,H] = IGRF00
-% MATLAB routine to load Schmidt-normalized coefficients
-% retrieved from ftp://nssdc.gsfc.nasa.gov/pub/models/igrf/
-% igrf00.dat
-%++++++++++++++++++++++++++++++++++++++++++
-% The number 1 is added to ALL subscripts since MATLAB can't have an array
-% index of 0. Units of Tesla
-%
-% MARTIN VISBECK, LDEO Feb 2000
-% !!! CAUTION when updating the values NEVER us a 0. where use 0.1 instead!!!
-% I use the zeros to throw unused coeffs away....!
-
-G(2,1) = -29615e-9;
-G(2,2) = -1728e-9; H(2,2) = 5186e-9;
-G(3,1) = -2267e-9; H(3,1) = 0.0;
-G(3,2) = 3072e-9; H(3,2) = -2478e-9;
-G(3,3) = 1672e-9; H(3,3) = -458e-9;
-G(4,1) = 1341e-9; H(4,1) = 0.0;
-G(4,2) = -2290e-9; H(4,2) = -227e-9;
-G(4,3) = 1253e-9; H(4,3) = 296e-9;
-G(4,4) = 715e-9; H(4,4) = -492e-9;
-G(5,1) = 935e-9; H(5,1) = .0;
-G(5,2) = 787e-9; H(5,2) = 272e-9;
-G(5,3) = 251e-9; H(5,3) = -232e-9;
-G(5,4) = -405e-9; H(5,4) = 119e-9;
-G(5,5) = 110e-9; H(5,5) = -304e-9;
-G(6,1) = -217e-9; H(6,1) = .0;
-G(6,2) = 351e-9; H(6,2) = 44e-9;
-G(6,3) = 222e-9; H(6,3) = 172e-9;
-G(6,4) = -131e-9; H(6,4) = -134e-9;
-G(6,5) = -169e-9; H(6,5) = -40e-9;
-G(6,6) = -12e-9; H(6,6) = 107e-9;
-G(7,1) = 72e-9; H(7,1) = .0;
-G(7,2) = 68e-9; H(7,2) = -17e-9;
-G(7,3) = 74e-9; H(7,3) = 64e-9;
-G(7,4) = -161e-9; H(7,4) = 65e-9;
-G(7,5) = -5e-9; H(7,5) = -61e-9;
-G(7,6) = 17e-9; H(7,6) = 1e-9;
-G(7,7) = -91e-9; H(7,7) = 44e-9;
-G(8,1) = 79e-9; H(8,1) = .0;
-G(8,2) = -74e-9; H(8,2) = -65e-9;
-G(8,3) = 0.1e-9; H(8,3) = -24e-9;
-G(8,4) = 33e-9; H(8,4) = 6e-9;
-G(8,5) = 9e-9; H(8,5) = 24e-9;
-G(8,6) = 7e-9; H(8,6) = 15e-9;
-G(8,7) = 8e-9; H(8,7) = -25e-9;
-G(8,8) = -2e-9; H(8,8) = -6e-9;
-G(9,1) = 25e-9; H(9,1) = .0;
-G(9,2) = 6e-9; H(9,2) = 12e-9;
-G(9,3) = -9e-9; H(9,3) = -22e-9;
-G(9,4) = -8e-9; H(9,4) = 8e-9;
-G(9,5) = -17e-9; H(9,5) = -21e-9;
-G(9,6) = 9e-9; H(9,6) = 15e-9;
-G(9,7) = 7e-9; H(9,7) = 9e-9;
-G(9,8) = -8e-9; H(9,8) = -16e-9;
-G(9,9) = -7e-9; H(9,9) = -3e-9;
-G(10,1) = 5e-9; H(10,1) = .0;
-G(10,2) = 9e-9; H(10,2) = -20e-9;
-G(10,3) = 3e-9; H(10,3) = 13e-9;
-G(10,4) = -8e-9; H(10,4) = 12e-9;
-G(10,5) = 6e-9; H(10,5) = -6e-9;
-G(10,6) = -9e-9; H(10,6) = -8e-9;
-G(10,7) = -2e-9; H(10,7) = 9e-9;
-G(10,8) = 9e-9; H(10,8) = 4e-9;
-G(10,9) = -4e-9; H(10,9) = -8e-9;
-G(10,10) = -8e-9; H(10,10) = 5e-9;
-G(11,1) = -2e-9; H(11,1) = .0;
-G(11,2) = -6e-9; H(11,2) = 1e-9;
-G(11,3) = 2e-9; H(11,3) = 0.1e-9;
-G(11,4) = -3e-9; H(11,4) = 4e-9;
-G(11,5) = 0.1e-9; H(11,5) = 5e-9;
-G(11,6) = 4e-9; H(11,6) = -6e-9;
-G(11,7) = 1e-9; H(11,7) = -1e-9;
-G(11,8) = 2e-9; H(11,8) = -3e-9;
-G(11,9) = 4e-9; H(11,9) = 0.1e-9;
-G(11,10) = 0.1e-9; H(11,10) = -2e-9;
-G(11,11) = -1e-9; H(11,11) = -8e-9;
-
-% prepare compressed array
-
-g=reshape(G',11*11,1)*1e9;
-h=reshape(H',11*11,1)*1e9;
-
-gh=reshape([g,h]',2*11*11,1);
-% here is where zeros have a meaning....!!!!
-ii=find(gh==0);
-gh(ii)=[];
-
-return;
-
-
--- a/IGRF95.m Wed Jan 17 12:29:40 2018 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,95 +0,0 @@
-function [gh,G,H] = IGRF95
-% MATLAB routine to load Schmidt-normalized coefficients
-% retrieved from ftp://nssdc.gsfc.nasa.gov/pub/models/igrf/
-% igrf95.dat 1 Kb Mon Nov 13 00:00:00 1995
-% ? C.E. Barton, Revision of International Geomagnetic Reference
-% Field Released, EOS Transactions 77, #16, April 16, 1996.
-% The coefficients are from the 1995 International Geomagnetic Reference Field
-% Carlos Roithmayr, Jan. 22, 1997.
-%++++++++++++++++++++++++++++++++++++++++++
-% The number 1 is added to ALL subscripts since MATLAB can't have an array
-% index of 0. Units of Tesla
-%
-% MARTIN VISBECK, LDEO Feb 2000
-% !!! CAUTION when updating the values NEVER us a 0. where use 0.1 instead!!!
-% I use the zeros to throw unused coeffs away....!
-
-G(2,1) = -29682e-9;
-G(2,2) = -1789e-9; H(2,2) = 5318e-9;
-G(3,1) = -2197e-9; H(3,1) = 0.0;
-G(3,2) = 3074e-9; H(3,2) = -2356e-9;
-G(3,3) = 1685e-9; H(3,3) = -425e-9;
-G(4,1) = 1329e-9; H(4,1) = 0.0;
-G(4,2) = -2268e-9; H(4,2) = -263e-9;
-G(4,3) = 1249e-9; H(4,3) = 302e-9;
-G(4,4) = 769e-9; H(4,4) = -406e-9;
-G(5,1) = 941e-9; H(5,1) = .0;
-G(5,2) = 782e-9; H(5,2) = 262e-9;
-G(5,3) = 291e-9; H(5,3) = -232e-9;
-G(5,4) = -421e-9; H(5,4) = 98e-9;
-G(5,5) = 116e-9; H(5,5) = -301e-9;
-G(6,1) = -210e-9; H(6,1) = .0;
-G(6,2) = 352e-9; H(6,2) = 44e-9;
-G(6,3) = 237e-9; H(6,3) = 157e-9;
-G(6,4) = -122e-9; H(6,4) = -152e-9;
-G(6,5) = -167e-9; H(6,5) = -64e-9;
-G(6,6) = -26e-9; H(6,6) = 99e-9;
-G(7,1) = 66e-9; H(7,1) = .0;
-G(7,2) = 64e-9; H(7,2) = -16e-9;
-G(7,3) = 65e-9; H(7,3) = 77e-9;
-G(7,4) = -172e-9; H(7,4) = 67e-9;
-G(7,5) = 2e-9; H(7,5) = -57e-9;
-G(7,6) = 17e-9; H(7,6) = 4e-9;
-G(7,7) = -94e-9; H(7,7) = 28e-9;
-G(8,1) = 78e-9; H(8,1) = -.0;
-G(8,2) = -67e-9; H(8,2) = -77e-9;
-G(8,3) = 1e-9; H(8,3) = -25e-9;
-G(8,4) = 29e-9; H(8,4) = 3e-9;
-G(8,5) = 4e-9; H(8,5) = 22e-9;
-G(8,6) = 8e-9; H(8,6) = 16e-9;
-G(8,7) = 10e-9; H(8,7) = -23e-9;
-G(8,8) = -2e-9; H(8,8) = -3e-9;
-G(9,1) = 24e-9; H(9,1) = .0;
-G(9,2) = 4e-9; H(9,2) = 12e-9;
-G(9,3) = -1e-9; H(9,3) = -20e-9;
-G(9,4) = -9e-9; H(9,4) = 7e-9;
-G(9,5) = -14e-9; H(9,5) = -21e-9;
-G(9,6) = 4e-9; H(9,6) = 12e-9;
-G(9,7) = 5e-9; H(9,7) = 10e-9;
-G(9,8) = 0.1e-9; H(9,8) = -17e-9;
-G(9,9) = -7e-9; H(9,9) = -10e-9;
-G(10,1) = 4e-9; H(10,1) = .0;
-G(10,2) = 9e-9; H(10,2) = -19e-9;
-G(10,3) = 1e-9; H(10,3) = 15e-9;
-G(10,4) = -12e-9; H(10,4) = 11e-9;
-G(10,5) = 9e-9; H(10,5) = -7e-9;
-G(10,6) = -4e-9; H(10,6) = -7e-9;
-G(10,7) = -2e-9; H(10,7) = 9e-9;
-G(10,8) = 7e-9; H(10,8) = 7e-9;
-G(10,9) = 0.1e-9; H(10,9) = -8e-9;
-G(10,10) = -6e-9; H(10,10) = 1e-9;
-G(11,1) = -3e-9; H(11,1) = .0;
-G(11,2) = -4e-9; H(11,2) = 2e-9;
-G(11,3) = 2e-9; H(11,3) = 1e-9;
-G(11,4) = -5e-9; H(11,4) = 3e-9;
-G(11,5) = -2e-9; H(11,5) = 6e-9;
-G(11,6) = 4e-9; H(11,6) = -4e-9;
-G(11,7) = 3e-9; H(11,7) = 0.1e-9;
-G(11,8) = 1e-9; H(11,8) = -2e-9;
-G(11,9) = 3e-9; H(11,9) = 3e-9;
-G(11,10) = 3e-9; H(11,10) = -1e-9;
-G(11,11) = 0.1e-9; H(11,11) = -6e-9;
-
-% prepare compressed array
-
-g=reshape(G',11*11,1)*1e9;
-h=reshape(H',11*11,1)*1e9;
-
-gh=reshape([g,h]',2*11*11,1);
-% here is where zeros have a meaning....!!!!
-ii=find(gh==0);
-gh(ii)=[];
-
-return;
-
-
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/calc_shear3.m Sat Apr 10 08:03:07 2021 -0400
@@ -0,0 +1,462 @@
+function [ds,p,dr] = calc_shear3(d,p,ps,dr)
+% function [ds,p,dr] = calc_shear3(d,p,ps,dr)
+%
+% - compute shear profiles
+% - use only central difference
+% - use 2*std editing
+%
+% version 8 last change 06.08.2019
+
+% Martin Visbeck, LDEO, 3/7/97
+% some modifications and code cleanup GK, 16.05.2015 2-->3
+% reinstated the use of d.weight for the identification of shear pairs
+% GK, 03.12.2018 3-->4
+% reorganization of the dz handling GK, 12.12.2018 4-->5
+% added error output, catch different shear methods GK, 22.07.2019 5-->6
+% more text output
+% added new parameter ps.shear_throw_out_percent GK, 05.08.2019 6-->7
+% fixed bug in shear calculation GK, 06.08.2019 7-->8
+
+%======================================================================
+% C A L C _ S H E A R 3 . M
+% doc: Thu Sep 5 10:58:55 2019
+% dlm: Thu Sep 5 16:27:19 2019
+% (c) 2019 G. Krahmann
+% uE-Info: 31 1 NIL 0 0 72 0 2 4 NIL ofnI
+%======================================================================
+
+% CHANGES BY ANT:
+% Sep 5, 2019: - added local_weight to another case in if
+% - added defaults provided by Gerd
+% - removed 'messages' variable
+% - replaced all "replace" statements (required
+% removal of warning)
+% - added nmin() nmax() nmean() nstd()
+% - disabled code (replaced if 0 by if 1) calling
+% meanoutlier
+
+% Defaults from Gerd (email 9/5/2019):
+
+% average over data within how many standard deviations from median
+% this is the value used by the old calc_shear2.m
+% In calc_shear3.m the outlying fraction of data is discarded. Assuming a
+% normal distribution a value of stdf=2 converts to an outlier fraction of
+% 5%, stdf=3 converts to 1% and stdf=1 converts to 32%.
+% Internally in calc_shear3 this is converted so that shear_stdf does not
+% need to be changed.
+ps.shear_stdf = 2;
+
+
+% the minimum weight a bin must have to be accepted for shear
+ps.shear_weightmin = nan;
+
+
+% this one is similar to shear_weight_min, but gives the percentage of
+% values thrown out. The routine throws out the xx% values with the
+% lowest weight. The weights are the correlation based ones.
+ps.shear_throw_out_percent = 10;
+
+
+%
+% general function info
+%
+disp(' ')
+disp(['CALC_SHEAR3: calculate a baroclinic velocity profile based on shears only'])
+
+
+%
+% decide whether central differences or simple differences are to be used in
+% the shear calculation
+%
+% 1: simple differences
+% 2: central differences
+%
+diff_type = 2;
+
+
+%
+% resolution of final shear profile in meter
+%
+dz = ps.dz;
+shear_dz = dz * diff_type;
+disp([' Averaging shear profile over ',num2str(shear_dz),' m intervals'])
+
+
+%
+% Discard a certain amount of data as suspected outliers that
+% for relatively small numbers of values will skew the mean.
+% In the old calc_shear2.m the allowed range was the std of the shears
+% times the set factor stdf. Default was stdf=2 .
+% In calc_shear3.m the outlier determination is iterative and thus a bit
+% safer in calculation. The allowed range is set as a fraction of the whole
+% population. For gaussian distributions stdf converts to this fraction.
+% As usually stdf is not varied much, we have implemented a lookup list
+% here.
+%
+use_new_outlier = 1;
+stdf = ps.shear_stdf;
+if use_new_outlier==1
+ fracs = 1 - [31.8,13.4,4.5,1.3,0.3]/100;
+ stdfs = [1,1.5,2,2.5,3];
+ [dummy,ind] = min(abs(stdf-stdfs));
+ frac = fracs(ind);
+end
+disp([' Maximum allowed std within calculation intervals : ',num2str(stdf)])
+disp([' Data deviating more from the median will be discarded.'])
+
+
+%
+% check if only one istrument is to be used
+%
+if ps.up_dn_looker==2
+ % down looker only
+ d.weight(d.izu,:)=nan;
+elseif ps.up_dn_looker==3
+ % up looker only
+ d.weight(d.izd,:)=nan;
+end
+
+
+%
+% Apply a weight threshold to shear data.
+%
+% There are two variants.
+% First one is used when ps.shear_throw_out_percent is not NaN.
+% This one throws out the xx% shears with the lowest correlation-derived weights.
+% Second one is used when ps.shear_weightmin is not NaN.
+% This one sets a minimum weight to keep the shears.
+%
+% First one is now default with 10%.
+%
+disp([' Correlation-derived weights range from ',num2str(nmin(d.weight(:))),' to ',num2str(nmax(d.weight(:)))])
+if ~isnan(ps.shear_throw_out_percent)
+ local_weight = d.weight;
+ ind = find(isfinite(local_weight));
+ [dummy,ind2] = sort(local_weight(ind));
+ if length(ind2)>9
+ local_weight(ind(ind2(1:floor(length(ind2)/10)))) = nan;
+ end
+elseif ~isnan(ps.shear_weightmin)
+ disp([' Minimum weight ',num2str(ps.shear_weightmin),' for data to be used in shear calc.'])
+ local_weight = double(d.weight>ps.shear_weightmin);
+else
+ disp('> No weight criterion applied to raw shear data.')
+ disp('> You should set ps.shear_throw_out_percent or ps.shear_weightmin')
+ local_weight = d.weight;
+end
+disp([' Removed ',int2str(100-sum(isfinite(local_weight))/sum(isfinite(d.weight))*100),...
+ ' % of data with lowest weights from shear calculation.'])
+disp([' New weights range from ',num2str(nmin(local_weight(:))),' to ',num2str(nmax(local_weight(:)))])
+
+
+%
+% convert the weights to 1 and NaN
+%
+local_weight(find(local_weight <= 0)) = nan;
+local_weight(find(local_weight > 0)) = 1;
+
+%
+% compute shear
+%
+% two ways are offered here
+% first: central differences for the shears
+% second: single differences
+% the first is similar to the ways of the old calc_shear2.m
+%
+% central differences
+if diff_type==2
+ local_weight = [repmat(nan,1,size(local_weight,2));diff2(local_weight)+1;repmat(nan,1,size(local_weight,2))];
+ ushear = [NaN*d.ru(1,:);diff2(d.ru(:,:))./diff2(d.izm);NaN*d.ru(1,:)].*local_weight;
+ vshear = [NaN*d.rv(1,:);diff2(d.rv(:,:))./diff2(d.izm);NaN*d.rv(1,:)].*local_weight;
+ wshear = [NaN*d.rw(1,:);diff2(d.rw(:,:))./diff2(d.izm);NaN*d.rw(1,:)].*local_weight;
+ zshear = -d.izm;
+% single differences
+else
+ ushear = diff( d.ru.*local_weight )./diff(d.izm);
+ vshear = diff( d.rv.*local_weight )./diff(d.izm);
+ wshear = diff( d.rw.*local_weight )./diff(d.izm);
+ zshear = -(d.izm(1:end-1,:)+d.izm(2:end,:))/2;
+end
+ds.ushear = ushear;
+ds.vshear = vshear;
+ds.wshear = wshear;
+ds.zshear = zshear;
+
+
+%
+% set depth levels
+%
+z = dr.z;
+
+
+%
+% prepare shear solution result vectors
+%
+ds.usm = repmat(nan,length(z),1);
+ds.vsm = ds.usm;
+ds.wsm = ds.usm;
+ds.usmd = ds.usm;
+ds.vsmd = ds.usm;
+ds.use = ds.usm;
+ds.vse = ds.usm;
+ds.wse = ds.usm;
+ds.nn = ds.usm;
+ds.z = z;
+
+
+%
+% loop over depth levels and calculate the average shear at that level
+%
+% in the case of central differences this is oversampled here
+% but by sticking with the same resolution it makes the results easier
+% to work with
+%
+for n=[1:length(z)]
+
+ i1 = find( ( abs( zshear - z(n) ) <= shear_dz/2 ) & isfinite( ushear + vshear ) );
+ ds.nn(n) = length(i1);
+ if ds.nn(n) > 2
+
+ % two ways to select outliers
+ % first: select all that are beyond a fixed range around the median
+ % second: iteratively reject the worst (largest distance from mean)
+ % until a fixed fraction is rejected
+ % the second is usually the safer calculation but is a bit slower
+ if 1
+ usmm = median( ushear(i1) );
+ ussd1 = std( ushear(i1) );
+ vsmm = median( vshear(i1) );
+ vssd1 = std( vshear(i1) );
+ wsmm = median( wshear(i1) );
+ wssd1 = std( wshear(i1) );
+ ii1 = i1( find(abs(ushear(i1)-usmm)<stdf*ussd1) );
+ ii2 = i1( find(abs(vshear(i1)-vsmm)<stdf*vssd1) );
+ ii3 = i1( find(abs(wshear(i1)-wsmm)<stdf*wssd1) );
+ else
+ [dummy,ii1] = meanoutlier(ushear(i1),frac);
+ [dummy,ii2] = meanoutlier(vshear(i1),frac);
+ [dummy,ii3] = meanoutlier(wshear(i1),frac);
+ ii1 = i1(ii1);
+ ii2 = i1(ii2);
+ ii3 = i1(ii3);
+ end
+
+ % two ways of calculating the mean and std of the selected shears
+ % first: if there is a rejected one in any of u,v,w shears then use it
+ % for non of the calculations
+ % second: if there is a rejected one in any of u,v,w shears then use it
+ % only in u,v or w calculations
+ % the second one is the one used by the old calc_shear2.m
+ % but to me this does not make sense, GK May 2015
+ if 1
+ dummy = zeros(size(ushear));
+ dummy(ii1) = 1;
+ dummy(ii2) = dummy(ii2)+1;
+ dummy(ii3) = dummy(ii3)+1;
+ ii = find(dummy==3);
+ if length(ii)>1
+
+ ds.usm(n) = mean(ushear(ii));
+ ds.usmd(n) = median(ushear(ii));
+ ds.use(n) = std(ushear(ii));
+ ds.ii(n) = length(ii);
+
+ ds.vsm(n) = mean(vshear(ii));
+ ds.vsmd(n) = median(vshear(ii));
+ ds.vse(n) = std(vshear(ii));
+
+ ds.wsm(n) = mean(wshear(ii));
+ ds.wsmd(n) = median(wshear(ii));
+ ds.wse(n) = std(wshear(ii));
+
+ % debugging plot
+ if 0
+ figure(3)
+ clf
+ subplot(3,1,1)
+ hist(ushear(ii),30)
+ hold on
+ ax = axis;
+ plot([1,1]*mean(ushear(ii)),ax(3:4),'r')
+ plot([1,1]*median(ushear(ii)),ax(3:4),'m')
+ ind = find(dr.z==z(n));
+ if ind<length(dr.u)
+ plot([1,1]*(dr.u(ind-1)-dr.u(ind+1))/20,ax(3:4),'g')
+ end
+ title(int2str(z(n)))
+ subplot(3,1,2)
+ hist(vshear(ii),30)
+ hold on
+ ax = axis;
+ plot([1,1]*mean(vshear(ii)),ax(3:4),'r')
+ plot([1,1]*median(vshear(ii)),ax(3:4),'m')
+ ind = find(dr.z==z(n));
+ if ind<length(dr.u)
+ plot([1,1]*(dr.v(ind-1)-dr.v(ind+1))/20,ax(3:4),'g')
+ end
+ subplot(3,1,3)
+ hist( zshear(ii)-z(n), 30 )
+ pause
+ end
+
+ end
+ else
+ if length(ii1)>1
+ ds.usm(n) = mean(ushear(ii1));
+ ds.usmd(n) = median(ushear(ii1));
+ ds.use(n) = std(ushear(ii1));
+ end
+ if length(ii2)>1
+ ds.vsm(n) = mean(vshear(ii2));
+ ds.vsmd(n) = median(vshear(ii2));
+ ds.vse(n) = std(vshear(ii2));
+ end
+ if length(ii3)>1
+ ds.wsm(n) = mean(wshear(ii3));
+ ds.wsmd(n) = median(wshear(ii3));
+ ds.wse(n) = nstd(wshear(ii3));
+ end
+ end
+ end
+
+end
+
+
+%
+% a debugging figure
+%
+if 0
+sfigure(3);
+clf
+orient tall
+subplot(1,2,1)
+plot(ushear,zshear,'b.','markersize',3)
+hold on
+plot(ds.usm,ds.z,'r')
+plot(ds.usmd,ds.z,'k')
+inv_shear_u = -diff(dr.u-mean(dr.u))/dz;
+plot(inv_shear_u,(z(1:end-1)+z(2:end))/2,'g')
+set(gca,'ydir','reverse')
+
+subplot(1,2,2)
+plot(vshear,zshear,'b.','markersize',3)
+hold on
+plot(ds.vsm,ds.z,'r')
+plot(ds.vsmd,ds.z,'k')
+inv_shear_v = -diff(dr.v-mean(dr.v))/dz;
+plot(inv_shear_v,(z(1:end-1)+z(2:end))/2,'g')
+set(gca,'ydir','reverse')
+
+sfigure(2)
+end
+
+
+
+%
+% integrate shear profile (from bottom up)
+%
+ds.usm(find(isnan(ds.usm))) = 0;
+ds.vsm(find(isnan(ds.vsm))) = 0;
+ds.wsm(find(isnan(ds.wsm))) = 0;
+%if length(ind)/length(ds.usm)*100>5
+% disp(['> Found ',num2str(length(ind)/length(ds.usm)*100),...
+% '% Nan in shear data. Integration result might be problematic.'])
+%end
+if 1
+ ds.ur = flipud(cumsum(flipud(ds.usm)))*dz;
+ ds.vr = flipud(cumsum(flipud(ds.vsm)))*dz;
+ ds.wr = flipud(cumsum(flipud(ds.wsm)))*dz;
+else
+ ds.ur = flipud(cumsum(flipud(ds.usmd)))*dz;
+ ds.vr = flipud(cumsum(flipud(ds.vsmd)))*dz;
+ ds.wr = flipud(cumsum(flipud(ds.wsmd)))*dz;
+end
+ds.ur = ds.ur-mean(ds.ur);
+ds.vr = ds.vr-mean(ds.vr);
+ds.wr = ds.wr-mean(ds.wr);
+
+
+%
+% This is a peculiar place for the single ping error estimate. But
+% as it is based on the variability in the data itself, it makes sense.
+% The assumption is that there should be basically zero shear in the
+% vertical velocities. At least it is so small as to be not detectable
+% here. Thus any variability in the vertical shear is caused by the
+% errors/noise of the measurement. Together with an angular conversion
+% factor this gives an error/noise value for the horizontal velocities.
+%
+if ~isfield(d,'zd')
+ dz2 = abs(nmean(diff(d.z)));
+else
+ dz2 = diff_type*abs(mean(diff(d.zd)));
+end
+if isfield(d,'down')
+ fac = 1/tan(d.down.Beam_angle*pi/180)*sqrt(2)*dz2;
+else
+ fac = 1/tan(d.up.Beam_angle*pi/180)*sqrt(2)*dz2;
+end
+ds.ensemble_vel_err = ds.wse*fac;
+dr.ensemble_vel_err = ds.wse*fac;
+
+
+%
+% store result and give text output
+%
+dr.u_shear_method = ds.ur;
+dr.v_shear_method = ds.vr;
+dr.w_shear_method = ds.wr;
+uds = nstd(dr.u-mean(dr.u)-ds.ur);
+vds = nstd(dr.v-mean(dr.v)-ds.vr);
+uvds = sqrt(uds^2+vds^2);
+disp([' Inversion average error : ',num2str(nmean( dr.uerr ) ),' m/s'])
+if uvds>nmean(dr.uerr)*1.5
+ error_increase_factor = 1/nmean(dr.uerr)*uvds/1.5;
+ warn = ('> Increasing error estimate because of elevated shear - inverse difference');
+ disp(warn)
+ disp(['> by a factor of ',num2str(error_increase_factor)])
+ disp(['> std of difference between regular and shear profile : ',num2str(uvds),' m/s'])
+ p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;
+ dr.uerr = dr.uerr * error_increase_factor;
+end
+disp([' Final average error : ',num2str(nmean( dr.uerr ) ),' m/s'])
+
+%--------------------------------------------------
+
+function x = diff2(x,k,dn)
+%DIFF2 Difference function. If X is a vector [x(1) x(2) ... x(n)],
+% then DIFF(X) returns a vector of central differences between
+% elements [x(3)-x(1) x(4)-x(2) ... x(n)-x(n-2)]. If X is a
+% matrix, the differences are calculated down each column:
+% DIFF(X) = X(3:n,:) - X(1:n-2,:).
+% DIFF(X,n) is the n'th difference function.
+
+% J.N. Little 8-30-85
+% Copyright (c) 1985, 1986 by the MathWorks, Inc.
+
+if nargin < 2, k = 1; end
+if nargin < 3, dn = 2; end
+for i=1:k
+ [m,n] = size(x);
+ if m == 1
+ x = x(1+dn:n) - x(1:n-dn);
+ else
+ x = x(1+dn:m,:) - x(1:m-dn,:);
+ end
+end
+
+%----------------------------------------------------------------------
+
+function m = nmin(d)
+ m = min(d(find(isfinite(d))));
+
+function m = nmax(d)
+ m = max(d(find(isfinite(d))));
+
+function m = nmean(d)
+ m = mean(d(find(isfinite(d))));
+
+function m = nstd(d)
+ m = std(d(find(isfinite(d))));
+
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/calc_shear3.m.hacked Sat Apr 10 08:03:07 2021 -0400
@@ -0,0 +1,425 @@
+function [ds,p,dr] = calc_shear3(d,p,ps,dr)
+% function [ds,p,dr] = calc_shear3(d,p,ps,dr)
+%
+% - compute shear profiles
+% - use only central difference
+% - use 2*std editing
+%
+% version 8 last change 06.08.2019
+
+% Martin Visbeck, LDEO, 3/7/97
+% some modifications and code cleanup GK, 16.05.2015 2-->3
+% reinstated the use of d.weight for the identification of shear pairs
+% GK, 03.12.2018 3-->4
+% reorganization of the dz handling GK, 12.12.2018 4-->5
+% added error output, catch different shear methods GK, 22.07.2019 5-->6
+% more text output
+% added new parameter ps.shear_throw_out_percent GK, 05.08.2019 6-->7
+% fixed bug in shear calculation GK, 06.08.2019 7-->8
+
+%======================================================================
+% C A L C _ S H E A R 3 . M
+% doc: Wed Sep 4 17:55:22 2019
+% dlm: Wed Sep 4 18:55:36 2019
+% (c) 2019 G. Krahmann, M. Visbeck
+% uE-Info: 400 0 NIL 0 0 72 0 2 4 NIL ofnI
+%======================================================================
+
+% CHANGES BY ANT:
+% Sep 4, 2019: - replaced messages by LDEO warn mechanism
+
+%
+% general function info
+%
+disp(' ')
+disp(['CALC_SHEAR3: calculate a baroclinic velocity profile based on shears only'])
+
+
+%
+% decide whether central differences or simple differences are to be used in
+% the shear calculation
+%
+% 1: simple differences
+% 2: central differences
+%
+diff_type = 2;
+
+
+%
+% resolution of final shear profile in meter
+%
+dz = ps.dz;
+shear_dz = dz * diff_type;
+disp([' Averaging shear profile over ',num2str(shear_dz),' m intervals'])
+
+
+%
+% Discard a certain amount of data as suspected outliers that
+% for relatively small numbers of values will skew the mean.
+% In the old calc_shear2.m the allowed range was the std of the shears
+% times the set factor stdf. Default was stdf=2 .
+% In calc_shear3.m the outlier determination is iterative and thus a bit
+% safer in calculation. The allowed range is set as a fraction of the whole
+% population. For gaussian distributions stdf converts to this fraction.
+% As usually stdf is not varied much, we have implemented a lookup list
+% here.
+%
+use_new_outlier = 1;
+ps.shear_stdf = 4; % arbitrary guess
+stdf = ps.shear_stdf;
+if use_new_outlier==1
+ fracs = 1 - [31.8,13.4,4.5,1.3,0.3]/100;
+ stdfs = [1,1.5,2,2.5,3];
+ [dummy,ind] = min(abs(stdf-stdfs));
+ frac = fracs(ind);
+end
+disp([' Maximum allowed std within calculation intervals : ',num2str(stdf)])
+disp([' Data deviating more from the median will be discarded.'])
+
+
+%
+% check if only one istrument is to be used
+%
+if ps.up_dn_looker==2
+ % down looker only
+ d.weight(d.izu,:)=nan;
+elseif ps.up_dn_looker==3
+ % up looker only
+ d.weight(d.izd,:)=nan;
+end
+
+
+%
+% Apply a weight threshold to shear data.
+%
+% There are two variants.
+% First one is used when ps.shear_throw_out_percent is not NaN.
+% This one throws out the xx% shears with the lowest correlation-derived weights.
+% Second one is used when ps.shear_weightmin is not NaN.
+% This one sets a minimum weight to keep the shears.
+%
+% First one is now default with 10%.
+%
+disp([' Correlation-derived weights range from ',num2str(min(d.weight(:))),' to ',num2str(max(d.weight(:)))])
+if 0 % ~isnan(ps.shear_throw_out_percent)
+ local_weight = d.weight;
+ ind = find(isfinite(local_weight));
+ [dummy,ind2] = sort(local_weight(ind));
+ if length(ind2)>9
+ local_weight(ind(ind2(1:floor(length(ind2)/10)))) = nan;
+ end
+elseif 0 % ~isnan(ps.shear_weightmin)
+ disp([' Minimum weight ',num2str(ps.shear_weightmin),' for data to be used in shear calc.'])
+ local_weight = double(d.weight>ps.shear_weightmin);
+else
+ disp('> No weight criterion applied to raw shear data.')
+ disp('> You should set ps.shear_throw_out_percent or ps.shear_weightmin')
+end
+local_weight = d.weight;
+disp([' Removed ',int2str(100-sum(isfinite(local_weight))/sum(isfinite(d.weight))*100),...
+ ' % of data with lowest weights from shear calculation.'])
+disp([' New weights range from ',num2str(min(local_weight(:))),' to ',num2str(max(local_weight(:)))])
+
+
+%
+% convert the weights to 1 and NaN
+%
+%local_weight = replace( local_weight, local_weight<=0, nan );
+%local_weight = replace( local_weight, local_weight>0, 1 );
+
+local_weight(find(local_weight <= 0)) = nan;
+local_weight(find(local_weight > 0)) = 1;
+
+
+%
+% compute shear
+%
+% two ways are offered here
+% first: central differences for the shears
+% second: single differences
+% the first is similar to the ways of the old calc_shear2.m
+%
+% central differences
+if diff_type==2
+ local_weight = [repmat(nan,1,size(local_weight,2));diff2(local_weight)+1;repmat(nan,1,size(local_weight,2))];
+ ushear = [NaN*d.ru(1,:);diff2(d.ru(:,:))./diff2(d.izm);NaN*d.ru(1,:)].*local_weight;
+ vshear = [NaN*d.rv(1,:);diff2(d.rv(:,:))./diff2(d.izm);NaN*d.rv(1,:)].*local_weight;
+ wshear = [NaN*d.rw(1,:);diff2(d.rw(:,:))./diff2(d.izm);NaN*d.rw(1,:)].*local_weight;
+ zshear = -d.izm;
+% single differences
+else
+ ushear = diff( d.ru.*local_weight )./diff(d.izm);
+ vshear = diff( d.rv.*local_weight )./diff(d.izm);
+ wshear = diff( d.rw.*local_weight )./diff(d.izm);
+ zshear = -(d.izm(1:end-1,:)+d.izm(2:end,:))/2;
+end
+ds.ushear = ushear;
+ds.vshear = vshear;
+ds.wshear = wshear;
+ds.zshear = zshear;
+
+
+%
+% set depth levels
+%
+z = dr.z;
+
+
+%
+% prepare shear solution result vectors
+%
+ds.usm = repmat(nan,length(z),1);
+ds.vsm = ds.usm;
+ds.wsm = ds.usm;
+ds.usmd = ds.usm;
+ds.vsmd = ds.usm;
+ds.use = ds.usm;
+ds.vse = ds.usm;
+ds.wse = ds.usm;
+ds.nn = ds.usm;
+ds.z = z;
+
+
+%
+% loop over depth levels and calculate the average shear at that level
+%
+% in the case of central differences this is oversampled here
+% but by sticking with the same resolution it makes the results easier
+% to work with
+%
+for n=[1:length(z)]
+
+ i1 = find( ( abs( zshear - z(n) ) <= shear_dz/2 ) & isfinite( ushear + vshear ) );
+ ds.nn(n) = length(i1);
+ if ds.nn(n) > 2
+
+ % two ways to select outliers
+ % first: select all that are beyond a fixed range around the median
+ % second: iteratively reject the worst (largest distance from mean)
+ % until a fixed fraction is rejected
+ % the second is usually the safer calculation but is a bit slower
+ if 1
+ usmm = median( ushear(i1) );
+ ussd1 = std( ushear(i1) );
+ vsmm = median( vshear(i1) );
+ vssd1 = std( vshear(i1) );
+ wsmm = median( wshear(i1) );
+ wssd1 = std( wshear(i1) );
+ ii1 = i1( find(abs(ushear(i1)-usmm)<stdf*ussd1) );
+ ii2 = i1( find(abs(vshear(i1)-vsmm)<stdf*vssd1) );
+ ii3 = i1( find(abs(wshear(i1)-wsmm)<stdf*wssd1) );
+ else
+ [dummy,ii1] = meanoutlier(ushear(i1),frac);
+ [dummy,ii2] = meanoutlier(vshear(i1),frac);
+ [dummy,ii3] = meanoutlier(wshear(i1),frac);
+ ii1 = i1(ii1);
+ ii2 = i1(ii2);
+ ii3 = i1(ii3);
+ end
+
+ % two ways of calculating the mean and std of the selected shears
+ % first: if there is a rejected one in any of u,v,w shears then use it
+ % for non of the calculations
+ % second: if there is a rejected one in any of u,v,w shears then use it
+ % only in u,v or w calculations
+ % the second one is the one used by the old calc_shear2.m
+ % but to me this does not make sense, GK May 2015
+ if 1
+ dummy = zeros(size(ushear));
+ dummy(ii1) = 1;
+ dummy(ii2) = dummy(ii2)+1;
+ dummy(ii3) = dummy(ii3)+1;
+ ii = find(dummy==3);
+ if length(ii)>1
+
+ ds.usm(n) = mean(ushear(ii));
+ ds.usmd(n) = median(ushear(ii));
+ ds.use(n) = std(ushear(ii));
+ ds.ii(n) = length(ii);
+
+ ds.vsm(n) = mean(vshear(ii));
+ ds.vsmd(n) = median(vshear(ii));
+ ds.vse(n) = std(vshear(ii));
+
+ ds.wsm(n) = mean(wshear(ii));
+ ds.wsmd(n) = median(wshear(ii));
+ ds.wse(n) = std(wshear(ii));
+
+ % debugging plot
+ if 0
+ figure(3)
+ clf
+ subplot(3,1,1)
+ hist(ushear(ii),30)
+ hold on
+ ax = axis;
+ plot([1,1]*mean(ushear(ii)),ax(3:4),'r')
+ plot([1,1]*median(ushear(ii)),ax(3:4),'m')
+ ind = find(dr.z==z(n));
+ if ind<length(dr.u)
+ plot([1,1]*(dr.u(ind-1)-dr.u(ind+1))/20,ax(3:4),'g')
+ end
+ title(int2str(z(n)))
+ subplot(3,1,2)
+ hist(vshear(ii),30)
+ hold on
+ ax = axis;
+ plot([1,1]*mean(vshear(ii)),ax(3:4),'r')
+ plot([1,1]*median(vshear(ii)),ax(3:4),'m')
+ ind = find(dr.z==z(n));
+ if ind<length(dr.u)
+ plot([1,1]*(dr.v(ind-1)-dr.v(ind+1))/20,ax(3:4),'g')
+ end
+ subplot(3,1,3)
+ hist( zshear(ii)-z(n), 30 )
+ pause
+ end
+
+ end
+ else
+ if length(ii1)>1
+ ds.usm(n) = mean(ushear(ii1));
+ ds.usmd(n) = median(ushear(ii1));
+ ds.use(n) = std(ushear(ii1));
+ end
+ if length(ii2)>1
+ ds.vsm(n) = mean(vshear(ii2));
+ ds.vsmd(n) = median(vshear(ii2));
+ ds.vse(n) = std(vshear(ii2));
+ end
+ if length(ii3)>1
+ ds.wsm(n) = mean(wshear(ii3));
+ ds.wsmd(n) = median(wshear(ii3));
+ ds.wse(n) = nstd(wshear(ii3));
+ end
+ end
+ end
+
+end
+
+
+%
+% a debugging figure
+%
+if 0
+sfigure(3);
+clf
+orient tall
+subplot(1,2,1)
+plot(ushear,zshear,'b.','markersize',3)
+hold on
+plot(ds.usm,ds.z,'r')
+plot(ds.usmd,ds.z,'k')
+inv_shear_u = -diff(dr.u-mean(dr.u))/dz;
+plot(inv_shear_u,(z(1:end-1)+z(2:end))/2,'g')
+set(gca,'ydir','reverse')
+
+subplot(1,2,2)
+plot(vshear,zshear,'b.','markersize',3)
+hold on
+plot(ds.vsm,ds.z,'r')
+plot(ds.vsmd,ds.z,'k')
+inv_shear_v = -diff(dr.v-mean(dr.v))/dz;
+plot(inv_shear_v,(z(1:end-1)+z(2:end))/2,'g')
+set(gca,'ydir','reverse')
+
+sfigure(2)
+end
+
+
+
+%
+% integrate shear profile (from bottom up)
+%
+%%ds.usm = replace( ds.usm, isnan(ds.usm), 0 );
+%%ds.vsm = replace( ds.vsm, isnan(ds.vsm), 0 );
+%%ds.wsm = replace( ds.wsm, isnan(ds.wsm), 0 );
+ds.usm(find(isnan(ds.usm))) = 0;
+ds.vsm(find(isnan(ds.vsm))) = 0;
+ds.wsm(find(isnan(ds.wsm))) = 0;
+if length(ind)/length(ds.usm)*100>5
+ disp(['> Found ',num2str(length(ind)/length(ds.usm)*100),...
+ '% Nan in shear data. Integration result might be problematic.'])
+end
+if 1
+ ds.ur = flipud(cumsum(flipud(ds.usm)))*dz;
+ ds.vr = flipud(cumsum(flipud(ds.vsm)))*dz;
+ ds.wr = flipud(cumsum(flipud(ds.wsm)))*dz;
+else
+ ds.ur = flipud(cumsum(flipud(ds.usmd)))*dz;
+ ds.vr = flipud(cumsum(flipud(ds.vsmd)))*dz;
+ ds.wr = flipud(cumsum(flipud(ds.wsmd)))*dz;
+end
+ds.ur = ds.ur-mean(ds.ur);
+ds.vr = ds.vr-mean(ds.vr);
+ds.wr = ds.wr-mean(ds.wr);
+
+
+%
+% This is a peculiar place for the single ping error estimate. But
+% as it is based on the variability in the data itself, it makes sense.
+% The assumption is that there should be basically zero shear in the
+% vertical velocities. At least it is so small as to be not detectable
+% here. Thus any variability in the vertical shear is caused by the
+% errors/noise of the measurement. Together with an angular conversion
+% factor this gives an error/noise value for the horizontal velocities.
+%
+if ~isfield(d,'zd')
+ dz2 = abs(mean(diff(d.z)));
+else
+ dz2 = diff_type*abs(mean(diff(d.zd)));
+end
+if isfield(d,'down')
+ fac = 1/tan(d.down.Beam_angle*pi/180)*sqrt(2)*dz2;
+else
+ fac = 1/tan(d.up.Beam_angle*pi/180)*sqrt(2)*dz2;
+end
+ds.ensemble_vel_err = ds.wse*fac;
+dr.ensemble_vel_err = ds.wse*fac;
+
+
+%
+% store result and give text output
+%
+dr.u_shear_method = ds.ur;
+dr.v_shear_method = ds.vr;
+dr.w_shear_method = ds.wr;
+uds = std(dr.u-mean(dr.u)-ds.ur);
+vds = std(dr.v-mean(dr.v)-ds.vr);
+uvds = sqrt(uds^2+vds^2);
+disp([' Inversion average error : ',num2str(mean( dr.uerr ) ),' m/s'])
+if uvds>mean(dr.uerr)*1.5
+ error_increase_factor = 1/mean(dr.uerr)*uvds/1.5;
+ warn = ('> Increasing error estimate because of elevated shear - inverse difference');
+ disp(warn)
+ disp(['> by a factor of ',num2str(error_increase_factor)])
+ disp(['> std of difference between regular and shear profile : ',num2str(uvds),' m/s'])
+ p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;
+ dr.uerr = dr.uerr * error_increase_factor;
+end
+disp([' Final average error : ',num2str(mean( dr.uerr ) ),' m/s'])
+
+%--------------------------------------------------
+
+function x = diff2(x,k,dn)
+%DIFF2 Difference function. If X is a vector [x(1) x(2) ... x(n)],
+% then DIFF(X) returns a vector of central differences between
+% elements [x(3)-x(1) x(4)-x(2) ... x(n)-x(n-2)]. If X is a
+% matrix, the differences are calculated down each column:
+% DIFF(X) = X(3:n,:) - X(1:n-2,:).
+% DIFF(X,n) is the n'th difference function.
+
+% J.N. Little 8-30-85
+% Copyright (c) 1985, 1986 by the MathWorks, Inc.
+
+if nargin < 2, k = 1; end
+if nargin < 3, dn = 2; end
+for i=1:k
+ [m,n] = size(x);
+ if m == 1
+ x = x(1+dn:n) - x(1:n-dn);
+ else
+ x = x(1+dn:m,:) - x(1:m-dn,:);
+ end
+end
+
--- a/checkbtrk.m Wed Jan 17 12:29:40 2018 -0500
+++ b/checkbtrk.m Sat Apr 10 08:03:07 2021 -0400
@@ -5,9 +5,9 @@
%======================================================================
% C H E C K B T R K . M
% doc: Sat Jul 3 23:57:25 2004
-% dlm: Fri Oct 1 20:54:13 2010
+% dlm: Wed Sep 4 19:02:39 2019
% (c) 2004 ladcp@
-% uE-Info: 22 24 NIL 0 0 72 0 2 8 NIL ofnI
+% uE-Info: 24 0 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================
% CHANGES BY ANT:
@@ -20,6 +20,7 @@
% was set) but when all these bottom returns turned out
% to be outside the valid range (i.e. almost certainly
% false)
+% Sep 4, 2019: - fixed plotting BUGs reported by GK
if ~existf(p,'zbottom'), return , end
if ~isfinite(p.zbottom), return, end
@@ -254,7 +255,7 @@
subplot(222); hold on;
plot(rvabe,izm,'.k'); % all velocities
plot(rva,izm,'.g'); % edited velocities
- plot(rva(iz),izm(iz),'.r') % used-for-BT velocities
+ plot(rvabe(iz),izm(iz),'.r') % used-for-BT velocities
axis tight
ax=axis;
ax(1:2)=[-1 1]*0.8;
@@ -305,7 +306,7 @@
subplot(223); hold on;
plot(rwabe,izm,'.k'); % all velocities
plot(rwa,izm,'.g'); % edited velocities
- plot(rwa(iz),izm(iz),'.r'); % used-for-BT velocities
+ plot(rwabe(iz),izm(iz),'.r'); % used-for-BT velocities
axis tight
ax=axis;
ax(1:2)=[-1 1]*0.4;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/checkbtrk.m.orig Sat Apr 10 08:03:07 2021 -0400
@@ -0,0 +1,387 @@
+function p=checkbtrk(d,di,de,dr,p);
+% compare RDI to own bottom track
+%
+
+%======================================================================
+% C H E C K B T R K . M
+% doc: Sat Jul 3 23:57:25 2004
+% dlm: Wed Sep 4 18:19:17 2019
+% (c) 2004 ladcp@
+% uE-Info: 309 48 NIL 0 0 72 0 2 8 NIL ofnI
+%======================================================================
+
+% CHANGES BY ANT:
+% Jul 4, 2004: - adapted to [edit_data.m], i.e. plot unedited data
+% in black & edited in blue; also use unedited data
+% to calculate bias slope
+% Jul 16, 2004: - adapted to changes in [edit_data.m]
+% Mar 5, 2010: - debugging???
+% Oct 1, 2010: - BUG: code bombed when bottom was detected (i.e. p.zbottom
+% was set) but when all these bottom returns turned out
+% to be outside the valid range (i.e. almost certainly
+% false)
+% Sep 4, 2019: - fixed plotting BUGs reported by GK
+
+if ~existf(p,'zbottom'), return , end
+if ~isfinite(p.zbottom), return, end
+if ~existf(d,'bvel_own'), d.bvel_own=d.bvel*nan; d.hbot_own=d.hbot*nan; end
+if ~existf(d,'bvel_rdi'), d.bvel_rdi=d.bvel*nan; d.hbot_rdi=d.hbot*nan; end
+
+ ii=find(abs(d.z+p.zbottom)<max(p.btrk_range) & ...
+ abs(d.z+p.zbottom)>min(p.btrk_range));
+
+ if isempty(ii), p.zbottom = nan, return, end
+
+ if existf(d,'bvel_rdi')
+ bu_rdi=d.bvel_rdi(ii,1);
+ bv_rdi=d.bvel_rdi(ii,2);
+ bw_rdi=d.bvel_rdi(ii,3);
+ else
+ bu_rdi=d.bvel(ii,1)*nan;
+ bv_rdi=d.bvel(ii,2)*nan;
+ bw_rdi=d.bvel(ii,3)*nan;
+ end
+
+ if existf(d,'hbot_rdi')
+ bh_rdi=d.hbot_rdi(ii);
+ else
+ bh_rdi=d.hbot(ii)*nan;
+ end
+
+ u_o=interp1(dr.z,dr.u,-d.z(ii)');
+ u_c=-interp1(dr.tim',dr.uctd',d.time_jul(ii)');
+
+ bu_own=d.bvel_own(ii,1);
+ bv_own=d.bvel_own(ii,2);
+ bw_own=d.bvel_own(ii,3);
+ bh_own=d.hbot_own(ii);
+ v_o=interp1(dr.z,dr.v,-d.z(ii)');
+ v_c=-interp1(dr.tim',dr.vctd',d.time_jul(ii)');
+
+ rw=d.rw(d.izd,ii);
+ ij=find((d.izm(d.izd,ii)+p.zbottom)<d.zd(1));
+ rw(ij)=nan;
+ if existf(d,'wctd') ==-1
+ w_c=d.wctd(ii)';
+ disp(' use CTD pressure derived W')
+ else
+ w_c=medianan(rw,3)';
+ disp(' use reference layer W')
+ end
+ w_c2=w_c-[0;diff(w_c)/2 ];
+
+ bu_used=d.bvel(ii,1);
+ bv_used=d.bvel(ii,2);
+ bw_used=d.bvel(ii,3);
+ bh_used=d.hbot(ii);
+
+ blen=d.down.Cell_length/100;
+
+% extract raw velocity data near bottom
+rube = d.ru(d.izd,ii);
+rvbe = d.rv(d.izd,ii);
+rwbe = d.rw(d.izd,ii);
+
+% the same, but only velocities that passed editing
+ru = d.ru; ru(find(~isfinite(d.weight))) = NaN; ru = ru(d.izd,ii);
+rv = d.rv; rv(find(~isfinite(d.weight))) = NaN; rv = rv(d.izd,ii);
+rw = d.rw; rw(find(~isfinite(d.weight))) = NaN; rw = rw(d.izd,ii);
+
+% correct velocities with solved U_ctd
+
+rua=ru-repmat(u_c',[length(d.izd) 1]); ruabe=rube-repmat(u_c',[length(d.izd) 1]);
+rva=rv-repmat(v_c',[length(d.izd) 1]); rvabe=rvbe-repmat(v_c',[length(d.izd) 1]);
+rwa=rw-repmat(w_c',[length(d.izd) 1]); rwabe=rwbe-repmat(w_c',[length(d.izd) 1]);
+rwan=abs(rw)./abs(repmat(w_c',[length(d.izd) 1]));
+rwanbe=abs(rwbe)./abs(repmat(w_c',[length(d.izd) 1]));
+
+
+
+bua_rdi=bu_rdi-u_c;
+bva_rdi=bv_rdi-v_c;
+
+% try to use w inbetween pings for real bottom track
+bwa1_rdi=bw_rdi-w_c;
+bwa2_rdi=bw_rdi-w_c2;
+if stdnan(bwa1_rdi)>stdnan(bwa2_rdi)
+ bwa_rdi=bwa2_rdi;
+ disp(' use inbetween w-ref')
+else
+ bwa_rdi=bwa1_rdi;
+end
+
+
+bua_own=bu_own-u_c;
+bva_own=bv_own-v_c;
+bwa_own=bw_own-w_c;
+
+bua_used=bu_used-u_c;
+bva_used=bv_used-v_c;
+bwa_used=bw_used-w_c;
+
+% save bias and std of bottom track anomaly
+p.btrk_u_bias=medianan(bua_used,6);
+p.btrk_u_std=stdnan(bua_used);
+p.btrk_v_bias=medianan(bva_used,6);
+p.btrk_v_std=stdnan(bva_used);
+p.btrk_w_bias=medianan(bwa_used,6);
+p.btrk_w_std=stdnan(bwa_used);
+
+disp(['CHECKBTRK: check bottom track against U_ctd solution '])
+disp([' profiles within give acceptable range: ',int2str(length(ii))])
+disp([' U bias :',num3str(p.btrk_u_bias,6,3),' [m/s] std: ',...
+ num3str(p.btrk_u_std,5,3),' [m/s]'])
+disp([' V bias :',num3str(p.btrk_v_bias,6,3),' [m/s] std: ',...
+ num3str(p.btrk_v_std,5,3),' [m/s]'])
+disp([' W bias :',num3str(p.btrk_w_bias,6,3),' [m/s] std: ',...
+ num3str(p.btrk_w_std,5,3),' [m/s]'])
+
+if abs(p.btrk_u_bias)>0.1
+ warn=[' large U bottom track bias ',num2str(p.btrk_u_bias)];
+ disp(warn)
+ p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;
+end
+
+if abs(p.btrk_v_bias)>0.1
+ warn=[' large V bottom track bias ',num2str(p.btrk_v_bias)];
+ disp(warn)
+ p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;
+end
+
+% correct nominal depth with detected bottom depth
+izm=d.izm(d.izd,ii);
+zcor=(-d.hbot(ii)+d.z(ii));
+zcor_std=stdnan(zcor);
+p.btrk_rough=zcor_std;
+izm=izm-repmat(zcor,[length(d.izd) 1]);
+
+% compute bias W slope
+ dz=abs(diff(d.zd(1:2)));
+ wbias=[];
+ zbias=[];
+ nbias=[];
+ for i=1:6
+ dzi=dz*i/2;
+ ij=find(abs(izm+dzi)<(dz/4));
+ nok=sum(isfinite(rwanbe(ij)));
+ if nok>5
+ wbias=[wbias, 1-medianan(rwanbe(ij))];
+ zbias=[zbias, dzi];
+ nbias=[nbias, nok];
+ end
+ end
+ if exist('wbias')==1
+ p.wbslope=polyfit(zbias,wbias,1);
+ else
+ wbias=nan;
+ zbias=nan;
+ p.wbslope=[nan nan];
+ end
+
+disp([' W slope fact :',num3str(p.wbslope(1),6,4),' [1/m] lower W below bottom '])
+
+% compute 'wave' action
+p.btrk_wdiff=medianan(abs(diff(w_c)));
+disp([' W diff :',num3str(p.btrk_wdiff,6,4),' [m/s] ping to ping w difference '])
+
+
+
+disp([' H std :',num3str(zcor_std,6,1),' [m] large means bottom is rough/sloped'])
+
+p.btrk_tilt_mean=mean(d.tilt(ii));
+p.btrk_tilt_std=std(d.tilt(ii));
+disp([' Tilt mean :',num3str(p.btrk_tilt_mean,3,1),' +/- ',...
+ num3str(p.btrk_tilt_std,3,1),' [^o] '])
+
+
+% select bins used for own bottom track
+iz=find(abs(izm+blen*p.btrk_below)<blen);
+
+
+% plot results
+tx2=0.25;
+figure(13)
+clf
+
+ subplot(221); hold on;
+ plot(ruabe,izm,'.k'); % all velocities
+ plot(rua,izm,'.g'); % edited velocities
+ plot(ruabe(iz),izm(iz),'.r'); % used-for-BT velocities
+ axis tight;
+ ax=axis;
+ ax(1:2)=[-1 1]*0.8;
+ ax(3)=-160;
+ ax(4)=50;
+ vh=floor(ax(1)*20)/20+0.0125 : 0.025: ceil(ax(2)*20)/20;
+ plot(ax(1:2),[0 0],'-')
+ plot([0 0],ax(3:4),'k-','linewidth',1.5)
+ plot([0 0]+p.btrk_u_bias,ax(3:4),'r-','linewidth',1.5)
+ text(tx2,10,' Bottom')
+
+ if sum(isfinite(bua_rdi))>0
+ [y,x]=hist(bua_rdi,vh);
+ y([1,end])=0;
+ ys=sum(y);
+ [x,y]=stairs(x,y);
+ fill(x,y/max(y)*30-50,'-g')
+ plot([0 0]+medianan(bua_rdi,6),ax(3:4),'-g')
+ text(tx2,-40,['RDI n: ',int2str(ys)])
+ text(0.1,-10,[' bias ',num3str(medianan(bua_rdi,6),6,3)],'color','g')
+ end
+
+ [y,x]=hist(bua_own,vh);
+ y([1,end])=0;
+ ys=sum(y);
+ [x,y]=stairs(x,y);
+ fill(x,y/max(y)*30-100,'-b')
+ plot([0 0]+medianan(bua_own,6),ax(3:4),'-b')
+ text(tx2,-90,[' own n: ',int2str(ys)])
+ text(0.1,-60,[' bias ',num3str(medianan(bua_own,6),6,3)],'color','b')
+
+ [y,x]=hist(real(de.bvel-de.uctd(:,1).'),vh);
+ y([1,end])=0;
+ ys=sum(y);
+ [x,y]=stairs(x,y);
+ fill(x,y/max(y)*30-150,'-r')
+ text(tx2,-140,['Sup. ens. n: ',int2str(ys)])
+ text(0.1,-110,[' bias ',num3str(p.btrk_u_bias,6,3)],'color','r')
+ title([' U bot-tr mean U_{ctd} : ',num3str(-meannan(real(de.uctd(:,1))),4,2),...
+ ' [m/s] added'])
+ xlabel(' U [m/s]')
+ ylabel('depth [m]')
+ axis(ax)
+ text(-0.7,-154,['Bottom roughness: ',int2str(p.btrk_rough),' [m]'])
+
+subplot(222); hold on;
+ plot(rvabe,izm,'.k'); % all velocities
+ plot(rvabe,izm,'.g'); % edited velocities
+ plot(rvabe(iz),izm(iz),'.r') % used-for-BT velocities
+ axis tight
+ ax=axis;
+ ax(1:2)=[-1 1]*0.8;
+ ax(3)=-160;
+ ax(4)=50;
+ plot(ax(1:2),[0 0],'-')
+ plot([0 0],ax(3:4),'k-')
+ plot([0 0]+p.btrk_v_bias,ax(3:4),'r-','linewidth',1.5)
+
+ text(tx2,10,' Bottom')
+ if sum(isfinite(bva_rdi))>0
+ [y,x]=hist(bva_rdi,vh);
+ y([1,end])=0;
+ ys=sum(y);
+ [x,y]=stairs(x,y);
+ fill(x,y/max(y)*30-50,'-g')
+ plot([0 0]+medianan(bva_rdi,6),ax(3:4),'-g')
+ text(tx2,-40,[' RDI n: ',int2str(ys)])
+ text(0.1,-10,[' bias ',num3str(medianan(bva_rdi,6),6,3)],'color','g')
+ end
+
+ [y,x]=hist(bva_own,vh);
+ y([1,end])=0;
+ ys=sum(y);
+ [x,y]=stairs(x,y);
+ fill(x,y/max(y)*30-100,'-b')
+ plot([0 0]+medianan(bva_own,6),ax(3:4),'-b')
+ text(tx2,-90,[' own n: ',int2str(ys)])
+ text(0.1,-60,[' bias ',num3str(medianan(bva_own,6),6,3)],'color','b')
+
+ [y,x]=hist(imag(de.bvel-de.uctd(:,1).'),vh);
+ y([1,end])=0;
+ ys=sum(y);
+ [x,y]=stairs(x,y);
+ fill(x,y/max(y)*30-150,'-r')
+ text(tx2,-140,['Sup. ens. n: ',int2str(ys)])
+ text(0.1,-110,[' bias ',num3str(p.btrk_v_bias,6,3)],'color','r')
+ title([' V bot-tr mean V_{ctd} : ',num3str(-meannan(imag(de.uctd(:,1))),4,2),...
+ ' [m/s] added'])
+ xlabel(' V [m/s]')
+ ylabel('depth [m]')
+ text(-0.7,-154,['Tilt : ',num3str(p.btrk_tilt_mean,3,1),' +/- ',...
+ num3str(p.btrk_tilt_std,3,1),' [^o]'])
+ axis(ax)
+
+tx2=0.05;
+
+subplot(223); hold on;
+ plot(rwabe,izm,'.k'); % all velocities
+ plot(rwa,izm,'.g'); % edited velocities
+ plot(rwabe(iz),izm(iz),'.r'); % used-for-BT velocities
+ axis tight
+ ax=axis;
+ ax(1:2)=[-1 1]*0.4;
+ vh=floor(ax(1)*20)/20+0.005 : 0.01: ceil(ax(2)*20)/20;
+ ax(3)=-160;
+ ax(4)=50;
+ plot(ax(1:2),[0 0],'-')
+ plot([0 0],ax(3:4),'k-')
+ plot([0 0]+p.btrk_w_bias,ax(3:4),'r-','linewidth',1.5)
+
+ text(tx2,10,' Bottom')
+ if sum(isfinite(bwa_rdi))>0
+ [y,x]=hist(bwa_rdi,vh);
+ y([1,end])=0;
+ ys=sum(y);
+ [x,y]=stairs(x,y);
+ fill(x,y/max(y)*30-50,'-g')
+ plot([0 0]+medianan(bwa_rdi,6),ax(3:4),'-g')
+ text(tx2,-40,[' RDI n: ',int2str(ys)])
+ text(0.1,-10,[' bias ',num3str(medianan(bva_rdi,6),6,3)],'color','g')
+ end
+
+ [y,x]=hist(bwa_own,vh);
+ y([1,end])=0;
+ ys=sum(y);
+ [x,y]=stairs(x,y);
+ fill(x,y/max(y)*30-100,'-b')
+ plot([0 0]+medianan(bwa_own,6),ax(3:4),'-b')
+ text(tx2,-90,[' own n: ',int2str(ys)])
+ text(0.1,-60,[' bias ',num3str(medianan(bwa_own,6),6,3)],'color','b')
+
+ [y,x]=hist(di.bvel(3,:)-medianan(di.rw),vh);
+ y([1,end])=0;
+ ys=sum(y);
+ [x,y]=stairs(x,y);
+ fill(x,y/max(y)*30-150,'-r')
+ text(tx2,-140,['Sup. ens. n: ',int2str(ys)])
+ text(0.1,-110,[' bias ',num3str(p.btrk_w_bias,6,3)],'color','r')
+ text(-0.3,-154,['W diff: ',num3str(p.btrk_wdiff,6,3),' [m/s]'])
+
+ title([' W bottom track [m/s] '])
+ xlabel(' W [m/s]')
+ ylabel('depth [m]')
+ axis(ax)
+
+subplot(224); hold on;
+ plot(rwanbe,izm,'.k'); % all velocities
+ plot(rwan,izm,'.g'); % edited velocities
+ plot(rwan(iz),izm(iz),'.r'); % used-for-BT velocities
+ axis tight
+ ax=axis;
+ ax(1:2)=[0.5 1.2];
+ ax(3)=-60;
+ ax(4)=60;
+ title([' W bottom track factor '])
+ xlabel(' abs(W)/abs(W_{ref}) ')
+ ylabel('depth [m]')
+ axis(ax)
+ plot(ax(1:2),[0 0],'-')
+ plot([1 1],ax(3:4),'k-','linewidth',1.5)
+ text(tx2,10,' Bottom')
+ plot(1-polyval(p.wbslope,[40 0]),[-40 0],'-b')
+ plot(1-wbias,-zbias,'pb')
+ grid
+ text(0.54,-50,[' bias slope ',num3str(p.wbslope(1),6,3),' [1/m]'])
+ text(0.54,-58,[' offset ',num3str(p.wbslope(2),6,3) ])
+ text(0.54,50,[' beam angle ',int2str(d.down.Beam_angle)])
+ text(0.54,42,[' bin length ',int2str(d.down.Cell_length/100),' [m]'])
+ text(0.54,34,[' Pings/Ens ',int2str(d.down.Pings_per_Ensemble)])
+ text(0.54,26,[' Frequency ',int2str(d.down.Frequency),' [kHz]'])
+
+ text(0.54,10,[' Zbottom ',int2str(p.zbottom),' [m]'])
+
+
+ streamer([p.name,' Figure 13']);
+
+ orient tall
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/current Sat Apr 10 08:03:07 2021 -0400
@@ -0,0 +1,1 @@
+/Data/LADCP/Software/LDEO_IX/current
\ No newline at end of file
--- a/default.m Wed Jan 17 12:29:40 2018 -0500
+++ b/default.m Sat Apr 10 08:03:07 2021 -0400
@@ -1,9 +1,9 @@
%======================================================================
% D E F A U L T . M
% doc: Sat Jun 26 06:10:09 2004
-% dlm: Wed Jan 17 12:19:09 2018
+% dlm: Wed Sep 4 18:10:56 2019
% (c) 2004 ladcp@
-% uE-Info: 40 0 NIL 0 0 72 0 2 4 NIL ofnI
+% uE-Info: 26 48 NIL 0 0 72 0 2 4 NIL ofnI
%======================================================================
% CHANGES BY ANT:
@@ -22,6 +22,8 @@
% - changed version to IX_13beta
% Mar 29, 2017: - added saveplot_pdf
% Jan 17, 2018: - changed ersion to IX_13 and published
+% Sep 4, 2019: - changed p.btrk_mode from 3 to 2 (own only)
+% - updated to version IX_15beta
% LADCP processing software
% M. Visbeck. LDEO/2003
@@ -37,7 +39,7 @@
% the data
% structure ps.??? contains parameter for the solution
% structure att.??? contains attributes
-p.software='LDEO LADCP software: Version IX_13';
+p.software='LDEO LADCP software: Version IX_14beta2';
% file names
% f.ladcpdo is the ONLY required input
@@ -167,7 +169,7 @@
% 2 : use only own bottom track
% 3 : use RDI, if existent, own else (default)
% 0 : use not bottom track at all
-% p=setdefv(p,'btrk_mode',3);
+% p=setdefv(p,'btrk_mode',2);
% p.btrk_ts is in dB to detect bottom above bin1 level (for own btm track)
% p=setdefv(p,'btrk_ts',10);
--- a/getbtrack.m Wed Jan 17 12:29:40 2018 -0500
+++ b/getbtrack.m Sat Apr 10 08:03:07 2021 -0400
@@ -4,6 +4,17 @@
% create own bottom track in addition to the one used before
%
+%======================================================================
+% G E T B T R A C K . M
+% doc: Wed Sep 4 17:07:35 2019
+% dlm: Wed Sep 4 17:09:37 2019
+% (c) 2019 A.M. Thurnherr
+% uE-Info: 108 68 NIL 0 0 72 0 2 8 NIL ofnI
+%======================================================================
+
+% CHANGES BY ANT:
+% Sep 4, 2019: - fixed minor typo (GK suggestion)
+
% force own distance to RDI bottom track
p=setdefv(p,'bottomdist',0);
% "manual" bottom track for down looker
@@ -25,7 +36,7 @@
% at=0.039; % attenuation dB/m for 150 kHz
% at=0.062; % attenuation dB/m for 300 kHz
if d.down.Frequency==300,
- p=setdefv(p,'ts_att_dn',0.06);
+ p=setdefv(p,'ts_att_dn',0.062);
else
p=setdefv(p,'ts_att_dn',0.039);
end
@@ -50,7 +61,6 @@
if p.btrk_ts>0
disp(' using increased bottom echo amplitudes to create bottom track')
- % don't use first bin
fitb1=1;
zd=d.zd(fitb1:end);
ead=d.tg(d.izd(fitb1:end),:);
@@ -94,7 +104,7 @@
disp([' use ',num2str(p.btrk_below),...
' bins below maximum target strength for own bottom track velocity'])
- % locate acceptable bottom tracks (don't accept first/last bin)
+ % locate acceptable bottom tracks (don't accept first two and last bin)
ii=find(dts>p.btrk_ts & ...
zmead>min(p.btrk_range) & zmead<max(p.btrk_range) & ...
imeadbv<(nbin-1) & imeadbv>fitb1);
--- a/getdpthi.m Wed Jan 17 12:29:40 2018 -0500
+++ b/getdpthi.m Sat Apr 10 08:03:07 2021 -0400
@@ -1,13 +1,14 @@
%======================================================================
% G E T D P T H I . M
% doc: Wed Jan 7 16:25:26 2009
-% dlm: Fri Mar 5 15:51:15 2010
+% dlm: Fri Aug 30 12:40:38 2019
% (c) 2009 A.M. Thurnherr
-% uE-Info: 459 33 NIL 0 0 72 0 2 4 NIL ofnI
+% uE-Info: 12 0 NIL 0 0 72 0 2 4 NIL ofnI
%======================================================================
% CHANGES BY ANT:
% Jan 7, 2009: - tightened use of exist()
+% Aug 30, 2019: - BUG: missing pressure values cause problem in output
function [d,p]=getdpthi(d,p)
% function [d,p]=getdpthi(d,p)
@@ -221,7 +222,7 @@
d.z_ladcp=-zz;
dz=d.z_ladcp-d.z;
ii=find(isfinite(dz));
- p.ladcpr_CTD_depth_std=[mean(dz), std(dz)];
+ p.ladcpr_CTD_depth_std=[mean(dz(ii)), std(dz(ii))];
disp(' use CTD time series depth, will not do depth inversion ')
disp([' LADCP minus CTD depth mean: ',num2str(p.ladcpr_CTD_depth_std(1)),...
' std: ',num2str(p.ladcpr_CTD_depth_std(2))]);
--- a/geterr.m Wed Jan 17 12:29:40 2018 -0500
+++ b/geterr.m Sat Apr 10 08:03:07 2021 -0400
@@ -1,9 +1,9 @@
%======================================================================
% G E T E R R . M
% doc: Wed Jun 30 23:24:51 2004
-% dlm: Thu Dec 7 09:50:52 2017
+% dlm: Mon Jan 27 19:42:27 2020
% (c) 2004 ladcp@
-% uE-Info: 23 47 NIL 0 0 72 0 2 8 NIL ofnI
+% uE-Info: 25 73 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================
% MODIFICATIONS BY ANT:
@@ -21,6 +21,9 @@
% Jan 28, 2015: - BUG: figure legend typo
% Dec 7, 2017: - BUG: btmi was set to nan for P6 station 94; fixed with
% symptomatic work-around
+% Jan 27, 2020: - BUG: btmi was set to nan for JR195 stations 20 and 30;
+% fixed by skipping sub-plots if btmi is not finite
+
function l=geterr(ps,dr,d,iplot)
% function l=geterr(dr,d,iplot)
@@ -249,12 +252,14 @@
title(sprintf('U-err std: %.03f',meannan(stdnan(l.ru_err'))))
subplot(232)
- if ps.fig3_avgerr == 2
- plot(medianan(l.ru_err(:,1:btmi)')',-ib,'r',medianan(l.ru_err(:,btmi:end)')',-ib,'b');
- title('median(U-err) [r/b: down-/up-cast]')
- else
- plot(meanan(l.ru_err(:,1:btmi)')',-ib,'r',meanan(l.ru_err(:,btmi:end)')',-ib,'b');
- title('mean(U-err) [r/b: down-/up-cast]')
+ if isfinite(btmi)
+ if ps.fig3_avgerr == 2
+ plot(medianan(l.ru_err(:,1:btmi)')',-ib,'r',medianan(l.ru_err(:,btmi:end)')',-ib,'b');
+ title('median(U-err) [r/b: down-/up-cast]')
+ else
+ plot(meanan(l.ru_err(:,1:btmi)')',-ib,'r',meanan(l.ru_err(:,btmi:end)')',-ib,'b');
+ title('mean(U-err) [r/b: down-/up-cast]')
+ end
end
set(gca,'XLim',[-0.05 0.05]);
set(gca,'Ylim',[-ib(end) -ib(1)]);
@@ -315,12 +320,14 @@
title(sprintf('V-err std: %.03f',meannan(stdnan(l.rv_err'))))
subplot(235)
- if ps.fig3_avgerr == 2
- plot(medianan(l.rv_err(:,1:btmi)')',-ib,'r',medianan(l.rv_err(:,btmi:end)')',-ib,'b');
- title('median(V-err) [r/b: down-/up-cast]')
- else
- plot(meanan(l.rv_err(:,1:btmi)')',-ib,'r',meanan(l.rv_err(:,btmi:end)')',-ib,'b');
- title('mean(V-err) [r/b: down-/up-cast]')
+ if isfinite(btmi)
+ if ps.fig3_avgerr == 2
+ plot(medianan(l.rv_err(:,1:btmi)')',-ib,'r',medianan(l.rv_err(:,btmi:end)')',-ib,'b');
+ title('median(V-err) [r/b: down-/up-cast]')
+ else
+ plot(meanan(l.rv_err(:,1:btmi)')',-ib,'r',meanan(l.rv_err(:,btmi:end)')',-ib,'b');
+ title('mean(V-err) [r/b: down-/up-cast]')
+ end
end
set(gca,'XLim',[-0.05 0.05]);
set(gca,'Ylim',[-ib(end) -ib(1)]);
--- a/getinv.m Wed Jan 17 12:29:40 2018 -0500
+++ b/getinv.m Sat Apr 10 08:03:07 2021 -0400
@@ -8,9 +8,9 @@
%======================================================================
% G E T I N V . M
% doc: Thu Jun 17 15:36:21 2004
-% dlm: Thu Jul 23 14:06:09 2015
+% dlm: Wed Sep 4 16:55:50 2019
% (c) 2004 ladcp@
-% uE-Info: 292 1 NIL 0 0 72 0 2 8 NIL ofnI
+% uE-Info: 35 55 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================
% CHANGE HISTORY:
@@ -32,22 +32,17 @@
% Jul 28, 2014: - modified how to specify smallfac
% Aug 9, 2014: - modified how to specify dragfac (+ve for fixed; -ve for scaled Martin's default)
% Jul 23, 2015: - commented on bug below (#!#)
+% - rounded default ps.dz (GK suggestion)
if nargin<5, iplot=0; end
%### start by defining all processing parameter which are not set by the user
disp('GETINV: compute best velocity profile')
-% resolution of final profile in meter
-ps=setdefv(ps,'dz',medianan(abs(diff(di.izm(:,1)))));
-% how strong shall the bottom track influence the solution
+ps=setdefv(ps,'dz',round(medianan(abs(diff(di.izm(:,1))))));
ps=setdefv(ps,'botfac',1);
-% how strong shall the SADCP velocity influence the solution
ps=setdefv(ps,'sadcpfac',1);
-% how stong shall the GPS positon influence the barotropic solution
ps=setdefv(ps,'barofac',1);
-% how strong do you believe in the simple wire dynamics (not very good yet)
-%ps=setdefv(ps,'dragfac',tanh(p.maxdepth/3000)*0.2);
ps=setdefv(ps,'dragfac',0);
% negative values weigh Martin's original default
Binary file igrf12coeffs.xls has changed
--- a/loadctd.m Wed Jan 17 12:29:40 2018 -0500
+++ b/loadctd.m Sat Apr 10 08:03:07 2021 -0400
@@ -1,9 +1,9 @@
%======================================================================
% L O A D C T D . M
% doc: Sat Jun 26 15:56:43 2004
-% dlm: Thu May 28 07:51:39 2015
+% dlm: Tue Jan 28 13:22:51 2020
% (c) 2004 M. Visbeck & A. Thurnherr
-% uE-Info: 169 24 NIL 0 0 72 0 2 8 NIL ofnI
+% uE-Info: 94 45 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================
function [d,p]=loadctd(f,d,p)
@@ -85,6 +85,13 @@
% Mar 21, 2014: - BUG: f.ctd_time_base used p.ctd_time_base set as default
% May 27, 2015: - removed confusing diagnostic message regarding adjusting NAV time
% May 28, 2015: - added error message when there are no valid vertical velocities
+% Apr 18, 2018: - BUG: ADCP-time shift warning was meaningless with elapsed time_base
+% Sep 14, 2018: - BUG: code move in 2008 broke working with single GPS file for entire cruise
+% Jan 28, 2020: - I don't understand Sep 14, 2018 bug any more; fix for that bug
+% involved moving code to loadnav.m, which does not work because
+% during loadnav p.time_start and end are not known (LADCP turn on/off
+% times are used); present code works with SR1b repeat cruises, which
+% all have single gps files
% read SEABIRD ctd timeseries file
disp(['LOADCTD: load CTD time series ',f.ctd])
@@ -439,7 +446,7 @@
if abs(lag)<p.ctdmaxlag & co>p.ctdmincorr
disp(' adjust ADCP time to CTD time and shift depth record ')
d.time_jul=d.time_jul+lagdt;
- if lagdt*24*3600>10
+ if f.ctd_time_base~=0 && lagdt*24*3600>10
disp('WARNING WARNING WARNING')
warn=[' shifted ADCP timeseries by ',int2str(lagdt*24*3600),' seconds '];
disp(warn)
@@ -540,7 +547,7 @@
% at this stage
%----------------------------------------------------------------------
-if p.navdata
+if p.navdata %%%&& f.nav_time_base == 0
if min(d.navtime_jul)>max(d.time_jul) | max(d.navtime_jul)<min(d.time_jul)
disp('NAV timeseries does not overlap WRONG STATION????')
disp([' NAV-data : ',datestrj(d.navtime_jul(1)),' to ',datestrj(d.navtime_jul(end))])
@@ -649,7 +656,7 @@
%
% print julian date string
%
-a=datestr(b-julian([0 1 0 0 0 0]));
+a=datestr(b-julian([0 1 0 0 0 0]));
%
%==============================================================
--- a/loadnav.m Wed Jan 17 12:29:40 2018 -0500
+++ b/loadnav.m Sat Apr 10 08:03:07 2021 -0400
@@ -1,9 +1,9 @@
%======================================================================
% L O A D N A V . M
% doc: Thu Jun 17 18:01:50 2004
-% dlm: Thu Feb 18 13:38:50 2016
+% dlm: Tue Jan 28 13:24:01 2020
% (c) 2004 ladcp@
-% uE-Info: 182 7 NIL 0 0 72 0 2 8 NIL ofnI
+% uE-Info: 43 64 NIL 0 0 72 2 2 8 NIL ofnI
%======================================================================
% MODIFICATIONS BY ANT:
@@ -36,6 +36,11 @@
% Jan 3, 2011: - changed IGRF11 validity to end of 2015 (from 2010)
% Feb 18, 2016: - BUG: geomag year range check bombed in 2016
% - added p.interp_missing_GPS using code provided by Jay Hooper
+% Sep 14, 2018: - BUG: 2008 code move broke working with single GPS file for entire cruise
+% Sep 4, 2019: - adapted to GK new magdev.m
+% - added p.magdec_source
+% Jan 28, 2020: - 2008 code move is required because at the loadnav stage p.time_start and
+% end reflect the LADCP turn-on and -off times
function [d,p]=loadnav(f,d,p)
% function [d,p]=loadnav(f,d,p)
@@ -56,6 +61,12 @@
% INTERPOLATE MISSING GPS VALUES (Jay Hooper)
p=setdefv(p,'interp_missing_GPS',1);
+% MAGNETIC DECLINATION
+% There are three different sources for magnetic declination
+% - specify p.drot manually
+% - p.magdec_source = 1 % use external magdec program, if available, otherwise ...
+% - p.magdec_source = 2 % use [magdev.m] from Gerd Krahmann
+p=setdefv(p,'magdec_source',2);
% FILE LAYOUT
f = setdefv(f,'nav_header_lines',0);
@@ -182,6 +193,84 @@
end
end
+%----------------------------------------------------------------------
+% The following code was moved into [loadctd.m] in July 2008 to fix a
+% bug that most likely only occurs with elapsed times. In 2018, it was
+% discovered that the bug fix introduced another bug that only affects
+% data sets with a single GPS time series from which the per-station
+% information must be extracted. Therefore, the code was moved back
+% here but only called when not working with elapsed time.
+% In Jan 2020 I realized that the code cannot be here at all, because
+% p.time_start and end are still the LADCP turn-on and -off times.
+% Therefore, I disabled the code here and re-enabled it for all
+% time bases in [loadctd.m]. I do not remember what bug this now causes.
+%----------------------------------------------------------------------
+
+if 0
+ if min(d.navtime_jul)>max(d.time_jul) | max(d.navtime_jul)<min(d.time_jul)
+ disp('NAV timeseries does not overlap WRONG STATION????')
+ disp([' NAV-data : ',datestrj(d.navtime_jul(1)),' to ',datestrj(d.navtime_jul(end))])
+ disp([' ADCP-data : ',datestrj(d.time_jul(1)),' to ',datestrj(d.time_jul(end))])
+ disp(' will ignore nav data')
+ p.navdata = 0;
+ d.navtime_jul = d.time_jul; % make sure same length
+ d.slat = d.navtime_jul*NaN;
+ d.slon = d.slat;
+ else
+ if d.navtime_jul(1) > d.time_jul(1)
+ disp('NAV timeseries starts after ADCP timeseries: used first NAV value to patch ')
+ disp([' NAV-data : ',datestrj(d.navtime_jul(1)),' to ',datestrj(d.navtime_jul(end))])
+ disp([' ADCP-data : ',datestrj(d.time_jul(1)),' to ',datestrj(d.time_jul(end))])
+ d.navtime_jul(1) = d.time_jul(1);
+ end
+
+ if d.navtime_jul(end) < d.time_jul(end)
+ disp('NAV timeseries ends before ADCP timeseries: used last NAV value to patch ')
+ disp([' NAV-data : ',datestrj(d.navtime_jul(1)),' to ',datestrj(d.navtime_jul(end))])
+ disp([' ADCP-data : ',datestrj(d.time_jul(1)),' to ',datestrj(d.time_jul(end))])
+ d.navtime_jul(end) = d.time_jul(end);
+ end
+
+ % find valid
+ ii=find(diff(d.navtime_jul)>0);
+ ii=[ii;length(d.navtime_jul)];
+
+ % average over p.navtime_av days
+ dt2m=[-p.navtime_av:(1/3600/24):p.navtime_av]';
+ slon=medianan(interp1q(d.navtime_jul(ii),d.slon(ii),julian(p.time_start)+dt2m));
+ slat=medianan(interp1q(d.navtime_jul(ii),d.slat(ii),julian(p.time_start)+dt2m));
+keyboard
+ p.nav_start=[fix(slat), (slat-fix(slat))*60, fix(slon), (slon-fix(slon))*60];
+ elon=medianan(interp1q(d.navtime_jul(ii),d.slon(ii),julian(p.time_end)+dt2m));
+ elat=medianan(interp1q(d.navtime_jul(ii),d.slat(ii),julian(p.time_end)+dt2m));
+ p.nav_end=[fix(elat), (elat-fix(elat))*60, fix(elon), (elon-fix(elon))*60];
+
+ % interpolate on RDI data
+ % this also shortens vectors to length of d.time_jul, which may be the
+ % only thing that is really needed
+ d.slon=interp1q(d.navtime_jul(ii),d.slon(ii),d.time_jul')';
+ d.slat=interp1q(d.navtime_jul(ii),d.slat(ii),d.time_jul')';
+
+ p.poss = [NaN NaN NaN NaN];
+ p.pose = [NaN NaN NaN NaN];
+ end
+
+ p=setdefv(p,'poss',[NaN NaN NaN NaN]);
+ [slat,slon] = pos2str(p.poss(1)+p.poss(2)/60,p.poss(3)+p.poss(4)/60);
+ disp([' update start pos from:',slat,' ',slon])
+ p.poss=p.nav_start;
+ [slat,slon] = pos2str(p.poss(1)+p.poss(2)/60,p.poss(3)+p.poss(4)/60);
+ disp([' to:',slat,' ',slon])
+
+ p=setdefv(p,'pose',[NaN NaN NaN NaN]);
+ [slat,slon] = pos2str(p.pose(1)+p.pose(2)/60,p.pose(3)+p.pose(4)/60);
+ disp([' update end pos from:',slat,' ',slon])
+ p.pose=p.nav_end;
+ [slat,slon] = pos2str(p.pose(1)+p.pose(2)/60,p.pose(3)+p.pose(4)/60);
+ disp([' to:',slat,' ',slon])
+
+end % if p.navdata
+
% =================================================================
% - at this point nav data is in d.navtime_jul, d.slon, d.slat
% and p.navdata is set to 1
@@ -189,16 +278,22 @@
% in [loadctd.m]
% =================================================================
-if ~isfinite(p.drot) % set magdecl
- [s,o] = system('magdec');
- if s == 1
- p.drot = geomag(f,meannan(d.navtime_jul),medianan(d.slat),medianan(d.slon));
- else
- warn = sprintf('"magdec" not found; using old magdev code with IGRF00',f.IGRF);
- disp(['WARNING: ' warn]);
- p.warn(size(p.warn,1)+1,1:length(warn))=warn;
- p.drot = magdev(medianan(d.slat),medianan(d.slon));
+if ~isfinite(p.drot) % set magdecl
+ if p.magdec_source == 1 % external magdec program
+ [s,o] = system('magdec');
+ if s == 1
+ p.drot = geomag(f,meannan(d.navtime_jul),medianan(d.slat),medianan(d.slon));
+ else
+ warn = sprintf('"magdec" not found; using magdev Matlab code');
+ disp(['WARNING: ' warn]);
+ p.warn(size(p.warn,1)+1,1:length(warn))=warn;
+ end
end
+
+ if ~isfinite(p.drot)
+ p.drot = magdev(medianan(d.slat),medianan(d.slon),0,p.time_start(1));
+ end
+
[d.ru,d.rv]=uvrot(d.ru,d.rv,p.drot);
[d.bvel(:,1),d.bvel(:,2)]=uvrot(d.bvel(:,1),d.bvel(:,2),p.drot);
disp(sprintf(' corrected for magnetic declination of %.1f deg',p.drot));
--- a/loadrdi.m Wed Jan 17 12:29:40 2018 -0500
+++ b/loadrdi.m Sat Apr 10 08:03:07 2021 -0400
@@ -13,9 +13,9 @@
%======================================================================
% L O A D R D I . M
% doc: Fri Jun 18 18:21:56 2004
-% dlm: Tue Feb 23 16:44:19 2016
+% dlm: Mon Jan 27 18:41:35 2020
% (c) 2004 ladcp@
-% uE-Info: 52 52 NIL 0 0 72 2 2 8 NIL ofnI
+% uE-Info: 54 98 NIL 0 0 72 2 2 8 NIL ofnI
%======================================================================
% CHANGES BY ANT
@@ -50,6 +50,8 @@
% Apr 15, 2015: - modified ambiguity-velocity warning as suggested by Diana Cardoso
% May 27, 2015: - clarified time-related warnings
% Feb 23, 2016: - clarified header id error message
+% Jan 27, 2020: - moved magdev call further back, where start time is known
+% (with GK's new magdev code, this howto is miraculously correct again - I think)
% p=setdefv(p,'pg_save',[1 2 3 4]);
% Default =3 for loadctd_whoi.
@@ -61,22 +63,6 @@
if nargin<2, p.name='unknown'; end
-if existf(p,'poss')
- if isfinite(sum(p.poss))
- drot=magdev(p.poss(1)+p.poss(2)/60, p.poss(3)+p.poss(4)/60);
- if existf(p,'drot')
- if isfinite(p.drot)
- disp([' found drot:',num2str(p.drot),' should be ',num2str(drot)])
- else
- p.drot=drot;
- end
- else
- p.drot=drot;
- end
- end
-end
-
-
if existf(f,'ladcpdo')==0
error([' need file name f.ladcpdo !!!!! '])
end
@@ -342,12 +328,26 @@
% p.time_end=gregoria(l.tim(it(end)));
end
-
-d.soundc=0;
+if existf(p,'poss')
+ if isfinite(sum(p.poss))
+ drot=magdev(p.poss(1)+p.poss(2)/60, p.poss(3)+p.poss(4)/60,0,p.time_start(1));
+ if existf(p,'drot')
+ if isfinite(p.drot)
+ disp([' found drot:',num2str(p.drot),' should be ',num2str(drot)])
+ else
+ p.drot=drot;
+ end
+ else
+ p.drot=drot;
+ end
+ end
+end
+
%
% rotate for magnetic deviation
%
+d.soundc=0;
if isfinite(p.drot)
[d.ru,d.rv]=uvrot(l.u(:,it),l.v(:,it),p.drot);
[ub,vb]=uvrot(l.ub(it),l.vb(it),p.drot);
--- a/magdev.m Wed Jan 17 12:29:40 2018 -0500
+++ b/magdev.m Sat Apr 10 08:03:07 2021 -0400
@@ -1,245 +1,306 @@
-function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,nmax,igrf);
-% function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,nmax,igrf);
-%
-% compute magnetic deviation
-% input:
-% flat latitude (degree)
-% flon longitude (degree)
-% elevkm elevation above mean geoid (km)
-% nmax number of harmonics
-% igrf (need file, default 'IGRF95')
-%
-% output:
-% dev mag. deviation (degree)
-%
-%
-% based of FORTRAN ROUTINE GEOMAG.FOR
-% more info under http://fdd.gsfc.nasa.gov/IGRF.html
-%
-% M. Visbeck, LDEO FEB 2000
-%
-
-if nargin<5
- igrf='IGRF95';
- igrf='IGRF00';
-end
-eval(['gh=',igrf,';'])
-
-if nargin<4
- nmax=10;
-end
-
-if nargin<3
- elevkm=0;
-end
-
-
-
-[x,y,z] = shval3(flat,flon,elevkm,nmax,gh);
-[d,i,h,f]=dihf (x,y,z);
-
-dev=d*180/pi;
-
-if nargout<1
- disp([' model ',igrf,' harmonics: ',int2str(nmax)])
- disp([' lat: ',num2str(flat)])
- disp([' lon: ',num2str(flon)])
- disp([' mag dev [deg]: ',num2str(d*180/pi)])
-end
-
-%=====================================================================
+function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,year);
+% function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,year);
+%
+% compute magnetic deviation
+%
+% based on IGRF12 model (observed until end of 2014, predicted until end of 2019)
+%
+% input: flat - latitude (degree)
+% flon - longitude (degree)
+% elevkm [0] - elevation above mean geoid (km)
+% year - decimal year
+%
+% output: dev - mag. deviation (degree)
+%
+%
+% based of FORTRAN ROUTINE GEOMAG.FOR
+% more info under http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
+%
+%
+% version 0.5 last change 13.07.2015
+
+
+% M. Visbeck, LDEO FEB 2000
+% rewritten for other epochs, G. Krahmann, IFM-GEOMAR Mar 2007
+%
+% modified for >=R14 GK, Jun 2007 0.1-->0.2
+% switch to IGRF11 GK, 08.11.2010 0.2-->0.3
+% replaced sind/cosd with sind/cosd GK, 31.05.2011 0.3-->0.4
+% switch to IGRF12 GK, 12.07.2015 0.4-->0.5
+
+%======================================================================
+% M A G D E V . M
+% doc: Wed Sep 4 18:17:31 2019
+% dlm: Wed Sep 4 18:18:16 2019
+% (c) 2019 G. Krahmann
+% uE-Info: 40 49 NIL 0 0 72 0 2 8 NIL ofnI
+%======================================================================
-function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh)
-% function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh)
-% ================================================================
-%
-% version 1.01
-%
-% calculates field components from spherical harmonic (sh)
-% models.
-%
-% input:
-% flat - north latitude, in degrees
-% flon - east longitude, in degrees
-% elevkm - elevation above mean sea level
-% nmax - maximum degree and order of coefficients
-% gh - schmidt quasi-normal internal spherical
-% harmonic coefficients
-% iext - external coefficients flag (= 0 if none)
-% ext - the three 1st-degree external coefficients
-% (not used if iext = 0)
-%
-% output:
-% x - northward component
-% y - eastward component
-% z - vertically-downward component
-%
-% based on subroutine 'igrf' by d. r. barraclough and
-% s. r. c. malin, report no. 71/1, institute of geological
-% sciences, u.k.
-%
-% norman w. peddie, u.s. geological survey, mail stop 964,
-% federal center, box 25046, denver, colorado 80225
-%
-% ================================================================
-% the required sizes of the arrays used in this subroutine
-% depend on the value of nmax. the minimum dimensions
-% needed are indicated in the table below. (note that this
-% version is dimensioned for nmax of 14 or less).
-%
-% ================================================================
-dtr = pi/180;
-r=elevkm;
-erad=6371.2;
-a2=40680925.;
-b2=40408588.;
-
-
-
-slat = sin(flat*dtr);
-aa = min(89.999,max(-89.999,flat));
-clat = cos(aa*dtr);
-sl(1) = sin(flon*dtr);
-cl(1) = cos(flon*dtr);
-
-x=0.0;
-y=0.0;
-z=0.0;
-sd = 0.0;
-cd = 1.0;
-n=0;
-l=1;
-m=1;
+% CHANGES BY ANT:
+% SEp 4, 2019: - changed sin/cos_d to sin/cosd
-npq = (nmax*(nmax+3))/2;
-aa = a2*clat*clat;
-bb = b2*slat*slat;
-cc = aa+bb;
-dd = sqrt(cc);
-r=sqrt(elevkm*(elevkm+2.0*dd)+(a2*aa+b2*bb)/cc);
-cd = (elevkm+dd)/r;
-sd = (a2-b2)/dd*slat*clat/r;
-aa = slat;
-slat = slat*cd-clat*sd;
-clat = clat*cd+aa*sd;
-ratio = erad/r;
-aa = sqrt(3.0);
-p(1) = 2.0*slat;
-p(2) = 2.0*clat;
-p(3) = 4.5*slat*slat-1.5;
-p(4) = 3.0*aa*clat*slat;
-q(1) = -clat;
-q(2) = slat;
-q(3) = -3.0*clat*slat;
-q(4) = aa*(slat*slat-clat*clat);
-
-
-for k = 1: npq
-
-if (n<m)
- m=0;
- n=n+1;
- rr = ratio^(n+2);
- fn=n;
-end;
-
-fm=m;
-
-if (k>=5)
- if (m==n)
- aa = sqrt(1.0-.5/fm);
- j=k-n-1;
- p(k) = (1.0+1.0/fm)*aa*clat*p(j);
- q(k) = aa*(clat*q(j)+slat/fm*p(j));
- sl(m) = sl(m-1)*cl(1)+cl(m-1)*sl(1);
- cl(m) = cl(m-1)*cl(1)-sl(m-1)*sl(1);
- else
- aa = sqrt(fn*fn-fm*fm);
- bb = sqrt((fn-1.0)^2-fm*fm)/aa;
- cc = (2.0*fn-1.0)/aa;
- i=k-n;
- j=k-2*n+1;
- p(k) = (fn+1.0)*(cc*slat/fn*p(i)-bb/(fn-1.0)*p(j));
- q(k) = cc*(slat*q(i)-clat/fn*p(i))-bb*q(j);
- end;
-end;
-
-aa = rr*gh(l);
-
-if (m==0)
- x=x+aa*q(k);
- z=z-aa*p(k);
- l=l+1;
-else
- bb = rr*gh(l+1);
- cc = aa*cl(m)+bb*sl(m);
- x=x+cc*q(k);
- z=z-cc*p(k);
- if (clat>0.0)
- y=y+(aa*sl(m)-bb*cl(m))*fm*p(k)/((fn+1.0)*clat);
- else
- y=y+(aa*sl(m)-bb*cl(m))*q(k)*slat;
- end;
- l=l+2;
-end;
-
-m=m+1;
-
-end;
-
-aa=x;
-x=x*cd+z*sd;
-z=z*cd-aa*sd;
-
-return;
-
-%==================================================================
-function [d,i,h,f] = dihf(x,y,z)
-% function [d,i,h,f] = dihf(x,y,z)
-% ===============================================================
-%
-% version 1.01
-%
-% computes the geomagnetic elements d, i, h, and f from
-% x, y, and z.
-%
-% input:
-% x - northward component
-% y - eastward component
-% z - vertically-downward component
-%
-% output:
-% d - declination
-% i - inclination
-% h - horizontal intensity
-% f - total intensity
-%
-% a. zunde
-% usgs, ms 964, box 25046 federal center, denver, co 80225
-%
-% ===============================================================
-sn=0.0001;
-
-% ---------------------------------------------------------------
-% if d and i cannot be determined, set equal to NaN
-% ---------------------------------------------------------------
-h=sqrt(x*x+y*y);
-f=sqrt(x*x+y*y+z*z);
-if (f<sn)
- d=NaN;
- i=NaN;
-else
- i=atan2(z,h);
- if (h<sn)
- d=NaN;
- else
- hpx = h+x;
- if (hpx<200)
- d=pi;
- else
- d=2.0*atan2(y,hpx);
- end;
- end;
-end;
-return;
-
+
+%
+% read the coefficients
+%
+fname = 'igrf12coeffs.xls';
+warning off % avoid non-sensical warning in >=R14
+gh = xlsread(fname);
+warning on
+
+
+%
+% determine the maximum order of polynomials
+%
+% this might need modification in future versions of the file
+%
+if year<2000
+ nmax = 10;
+else
+ nmax = 13;
+end
+
+
+%
+% interpolate the coefficients in time
+%
+warning off % another workaround for annoying >=R14 warnings
+if year <= gh(1,end-1)
+ gh2 = interp1( gh(1,3:end-1),gh(2:end,3:end-1)',year)';
+else
+ gh2 = interp1( gh(1,3:end-1),gh(2:end,3:end-1)',gh(1,end-1))';
+ ghc = gh2*0;
+ dummy = gh(2:end,end);
+ good = find(~isnan(dummy));
+ ghc(good) = dummy(good);
+ gh2 = gh2+ghc*(year-gh(1,end-1));
+end
+warning on
+
+
+%
+% set default elevation
+%
+if nargin<3
+ elevkm=0;
+end
+
+
+%
+% calculate field
+%
+[x,y,z] = shval3(flat,flon,elevkm,nmax,gh2);
+[d,i,h,f]=dihf (x,y,z);
+
+
+%
+% convert units
+%
+dev=d*180/pi;
+
+
+%
+% more detailed screen output, if output is not stored
+%
+if nargout<1
+ disp([' model ',fname,' harmonics: ',int2str(nmax)])
+ disp([' lat: ',num2str(flat)])
+ disp([' lon: ',num2str(flon)])
+ disp([' tim: ',num2str(year)])
+ disp([' mag dev [deg]: ',num2str(d*180/pi)])
+end
+
+%=====================================================================
+
+function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh)
+% function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh)
+% ================================================================
+%
+% version 1.01
+%
+% calculates field components from spherical harmonic (sh)
+% models.
+%
+% input:
+% flat - north latitude, in degrees
+% flon - east longitude, in degrees
+% elevkm - elevation above mean sea level
+% nmax - maximum degree and order of coefficients
+% gh - schmidt quasi-normal internal spherical
+% harmonic coefficients
+% iext - external coefficients flag (= 0 if none)
+% ext - the three 1st-degree external coefficients
+% (not used if iext = 0)
+%
+% output:
+% x - northward component
+% y - eastward component
+% z - vertically-downward component
+%
+% based on subroutine 'igrf' by d. r. barraclough and
+% s. r. c. malin, report no. 71/1, institute of geological
+% sciences, u.k.
+%
+% norman w. peddie, u.s. geological survey, mail stop 964,
+% federal center, box 25046, denver, colorado 80225
+%
+% ================================================================
+% the required sizes of the arrays used in this subroutine
+% depend on the value of nmax. the minimum dimensions
+% needed are indicated in the table below. (note that this
+% version is dimensioned for nmax of 14 or less).
+%
+% ================================================================
+r=elevkm;
+erad=6371.2;
+a2=40680925.;
+b2=40408588.;
+
+
+
+slat = sind(flat);
+aa = min(89.999,max(-89.999,flat));
+clat = cosd(aa);
+sl(1) = sind(flon);
+cl(1) = cosd(flon);
+
+x=0.0;
+y=0.0;
+z=0.0;
+sd = 0.0;
+cd = 1.0;
+n=0;
+l=1;
+m=1;
+
+npq = (nmax*(nmax+3))/2;
+aa = a2*clat*clat;
+bb = b2*slat*slat;
+cc = aa+bb;
+dd = sqrt(cc);
+r=sqrt(elevkm*(elevkm+2.0*dd)+(a2*aa+b2*bb)/cc);
+cd = (elevkm+dd)/r;
+sd = (a2-b2)/dd*slat*clat/r;
+aa = slat;
+slat = slat*cd-clat*sd;
+clat = clat*cd+aa*sd;
+ratio = erad/r;
+aa = sqrt(3.0);
+p(1) = 2.0*slat;
+p(2) = 2.0*clat;
+p(3) = 4.5*slat*slat-1.5;
+p(4) = 3.0*aa*clat*slat;
+q(1) = -clat;
+q(2) = slat;
+q(3) = -3.0*clat*slat;
+q(4) = aa*(slat*slat-clat*clat);
+
+
+for k = 1: npq
+
+if (n<m)
+ m=0;
+ n=n+1;
+ rr = ratio^(n+2);
+ fn=n;
+end;
+
+fm=m;
+
+if (k>=5)
+ if (m==n)
+ aa = sqrt(1.0-.5/fm);
+ j=k-n-1;
+ p(k) = (1.0+1.0/fm)*aa*clat*p(j);
+ q(k) = aa*(clat*q(j)+slat/fm*p(j));
+ sl(m) = sl(m-1)*cl(1)+cl(m-1)*sl(1);
+ cl(m) = cl(m-1)*cl(1)-sl(m-1)*sl(1);
+ else
+ aa = sqrt(fn*fn-fm*fm);
+ bb = sqrt((fn-1.0)^2-fm*fm)/aa;
+ cc = (2.0*fn-1.0)/aa;
+ i=k-n;
+ j=k-2*n+1;
+ p(k) = (fn+1.0)*(cc*slat/fn*p(i)-bb/(fn-1.0)*p(j));
+ q(k) = cc*(slat*q(i)-clat/fn*p(i))-bb*q(j);
+ end;
+end;
+
+aa = rr*gh(l);
+
+if (m==0)
+ x=x+aa*q(k);
+ z=z-aa*p(k);
+ l=l+1;
+else
+ bb = rr*gh(l+1);
+ cc = aa*cl(m)+bb*sl(m);
+ x=x+cc*q(k);
+ z=z-cc*p(k);
+ if (clat>0.0)
+ y=y+(aa*sl(m)-bb*cl(m))*fm*p(k)/((fn+1.0)*clat);
+ else
+ y=y+(aa*sl(m)-bb*cl(m))*q(k)*slat;
+ end;
+ l=l+2;
+end;
+
+m=m+1;
+
+end;
+
+aa=x;
+x=x*cd+z*sd;
+z=z*cd-aa*sd;
+
+return;
+
+%==================================================================
+function [d,i,h,f] = dihf(x,y,z)
+% function [d,i,h,f] = dihf(x,y,z)
+% ===============================================================
+%
+% version 1.01
+%
+% computes the geomagnetic elements d, i, h, and f from
+% x, y, and z.
+%
+% input:
+% x - northward component
+% y - eastward component
+% z - vertically-downward component
+%
+% output:
+% d - declination
+% i - inclination
+% h - horizontal intensity
+% f - total intensity
+%
+% a. zunde
+% usgs, ms 964, box 25046 federal center, denver, co 80225
+%
+% ===============================================================
+sn=0.0001;
+
+% ---------------------------------------------------------------
+% if d and i cannot be determined, set equal to NaN
+% ---------------------------------------------------------------
+h=sqrt(x*x+y*y);
+f=sqrt(x*x+y*y+z*z);
+if (f<sn)
+ d=NaN;
+ i=NaN;
+else
+ i=atan2(z,h);
+ if (h<sn)
+ d=NaN;
+ else
+ hpx = h+x;
+ if (hpx<200)
+ d=pi;
+ else
+ d=2.0*atan2(y,hpx);
+ end;
+ end;
+end;
+return;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/magdev.m.pre2019 Sat Apr 10 08:03:07 2021 -0400
@@ -0,0 +1,245 @@
+function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,nmax,igrf);
+% function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,nmax,igrf);
+%
+% compute magnetic deviation
+% input:
+% flat latitude (degree)
+% flon longitude (degree)
+% elevkm elevation above mean geoid (km)
+% nmax number of harmonics
+% igrf (need file, default 'IGRF95')
+%
+% output:
+% dev mag. deviation (degree)
+%
+%
+% based of FORTRAN ROUTINE GEOMAG.FOR
+% more info under http://fdd.gsfc.nasa.gov/IGRF.html
+%
+% M. Visbeck, LDEO FEB 2000
+%
+
+if nargin<5
+ igrf='IGRF95';
+ igrf='IGRF00';
+end
+eval(['gh=',igrf,';'])
+
+if nargin<4
+ nmax=10;
+end
+
+if nargin<3
+ elevkm=0;
+end
+
+
+
+[x,y,z] = shval3(flat,flon,elevkm,nmax,gh);
+[d,i,h,f]=dihf (x,y,z);
+
+dev=d*180/pi;
+
+if nargout<1
+ disp([' model ',igrf,' harmonics: ',int2str(nmax)])
+ disp([' lat: ',num2str(flat)])
+ disp([' lon: ',num2str(flon)])
+ disp([' mag dev [deg]: ',num2str(d*180/pi)])
+end
+
+%=====================================================================
+
+function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh)
+% function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh)
+% ================================================================
+%
+% version 1.01
+%
+% calculates field components from spherical harmonic (sh)
+% models.
+%
+% input:
+% flat - north latitude, in degrees
+% flon - east longitude, in degrees
+% elevkm - elevation above mean sea level
+% nmax - maximum degree and order of coefficients
+% gh - schmidt quasi-normal internal spherical
+% harmonic coefficients
+% iext - external coefficients flag (= 0 if none)
+% ext - the three 1st-degree external coefficients
+% (not used if iext = 0)
+%
+% output:
+% x - northward component
+% y - eastward component
+% z - vertically-downward component
+%
+% based on subroutine 'igrf' by d. r. barraclough and
+% s. r. c. malin, report no. 71/1, institute of geological
+% sciences, u.k.
+%
+% norman w. peddie, u.s. geological survey, mail stop 964,
+% federal center, box 25046, denver, colorado 80225
+%
+% ================================================================
+% the required sizes of the arrays used in this subroutine
+% depend on the value of nmax. the minimum dimensions
+% needed are indicated in the table below. (note that this
+% version is dimensioned for nmax of 14 or less).
+%
+% ================================================================
+dtr = pi/180;
+r=elevkm;
+erad=6371.2;
+a2=40680925.;
+b2=40408588.;
+
+
+
+slat = sin(flat*dtr);
+aa = min(89.999,max(-89.999,flat));
+clat = cos(aa*dtr);
+sl(1) = sin(flon*dtr);
+cl(1) = cos(flon*dtr);
+
+x=0.0;
+y=0.0;
+z=0.0;
+sd = 0.0;
+cd = 1.0;
+n=0;
+l=1;
+m=1;
+
+npq = (nmax*(nmax+3))/2;
+aa = a2*clat*clat;
+bb = b2*slat*slat;
+cc = aa+bb;
+dd = sqrt(cc);
+r=sqrt(elevkm*(elevkm+2.0*dd)+(a2*aa+b2*bb)/cc);
+cd = (elevkm+dd)/r;
+sd = (a2-b2)/dd*slat*clat/r;
+aa = slat;
+slat = slat*cd-clat*sd;
+clat = clat*cd+aa*sd;
+ratio = erad/r;
+aa = sqrt(3.0);
+p(1) = 2.0*slat;
+p(2) = 2.0*clat;
+p(3) = 4.5*slat*slat-1.5;
+p(4) = 3.0*aa*clat*slat;
+q(1) = -clat;
+q(2) = slat;
+q(3) = -3.0*clat*slat;
+q(4) = aa*(slat*slat-clat*clat);
+
+
+for k = 1: npq
+
+if (n<m)
+ m=0;
+ n=n+1;
+ rr = ratio^(n+2);
+ fn=n;
+end;
+
+fm=m;
+
+if (k>=5)
+ if (m==n)
+ aa = sqrt(1.0-.5/fm);
+ j=k-n-1;
+ p(k) = (1.0+1.0/fm)*aa*clat*p(j);
+ q(k) = aa*(clat*q(j)+slat/fm*p(j));
+ sl(m) = sl(m-1)*cl(1)+cl(m-1)*sl(1);
+ cl(m) = cl(m-1)*cl(1)-sl(m-1)*sl(1);
+ else
+ aa = sqrt(fn*fn-fm*fm);
+ bb = sqrt((fn-1.0)^2-fm*fm)/aa;
+ cc = (2.0*fn-1.0)/aa;
+ i=k-n;
+ j=k-2*n+1;
+ p(k) = (fn+1.0)*(cc*slat/fn*p(i)-bb/(fn-1.0)*p(j));
+ q(k) = cc*(slat*q(i)-clat/fn*p(i))-bb*q(j);
+ end;
+end;
+
+aa = rr*gh(l);
+
+if (m==0)
+ x=x+aa*q(k);
+ z=z-aa*p(k);
+ l=l+1;
+else
+ bb = rr*gh(l+1);
+ cc = aa*cl(m)+bb*sl(m);
+ x=x+cc*q(k);
+ z=z-cc*p(k);
+ if (clat>0.0)
+ y=y+(aa*sl(m)-bb*cl(m))*fm*p(k)/((fn+1.0)*clat);
+ else
+ y=y+(aa*sl(m)-bb*cl(m))*q(k)*slat;
+ end;
+ l=l+2;
+end;
+
+m=m+1;
+
+end;
+
+aa=x;
+x=x*cd+z*sd;
+z=z*cd-aa*sd;
+
+return;
+
+%==================================================================
+function [d,i,h,f] = dihf(x,y,z)
+% function [d,i,h,f] = dihf(x,y,z)
+% ===============================================================
+%
+% version 1.01
+%
+% computes the geomagnetic elements d, i, h, and f from
+% x, y, and z.
+%
+% input:
+% x - northward component
+% y - eastward component
+% z - vertically-downward component
+%
+% output:
+% d - declination
+% i - inclination
+% h - horizontal intensity
+% f - total intensity
+%
+% a. zunde
+% usgs, ms 964, box 25046 federal center, denver, co 80225
+%
+% ===============================================================
+sn=0.0001;
+
+% ---------------------------------------------------------------
+% if d and i cannot be determined, set equal to NaN
+% ---------------------------------------------------------------
+h=sqrt(x*x+y*y);
+f=sqrt(x*x+y*y+z*z);
+if (f<sn)
+ d=NaN;
+ i=NaN;
+else
+ i=atan2(z,h);
+ if (h<sn)
+ d=NaN;
+ else
+ hpx = h+x;
+ if (hpx<200)
+ d=pi;
+ else
+ d=2.0*atan2(y,hpx);
+ end;
+ end;
+end;
+return;
+
--- a/prepinv.m Wed Jan 17 12:29:40 2018 -0500
+++ b/prepinv.m Sat Apr 10 08:03:07 2021 -0400
@@ -1,9 +1,9 @@
%======================================================================
% P R E P I N V . M
% doc: Wed Jan 7 16:46:29 2009
-% dlm: Fri Apr 1 08:52:25 2016
+% dlm: Wed Sep 4 17:03:23 2019
% (c) 2009 A.M. Thurnherr
-% uE-Info: 17 29 NIL 0 0 72 0 2 8 NIL ofnI
+% uE-Info: 18 93 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================
% CHANGES BY ANT:
@@ -15,6 +15,7 @@
% available
% Jun 9, 2014: - improved messages
% Apr 1, 2016: - cosmetics
+% Sep 4, 2019: - BUG: superens std was calculated without removing means! (reported by GK)
function [di,p,d]=prepinv(d,p,dr)
% function [di,p,d]=prepinv(d,p,dr)
@@ -530,7 +531,7 @@
iav=round(length(ur)/200*p.avpercent);
ur=meshgrid(ur,ibin);
di.ru(:,im)=medianan([d.ru(:,i1).*w2-ur]',iav)'+ruav;
- rus=stdnan([d.ru(:,i1).*w2]')';
+ rus=stdnan([d.ru(:,i1).*w2-ur]')';
% V
vr=medianan(d.rv(izr,i1).*w);
rvav=meannan(vr);
@@ -540,7 +541,7 @@
vr=meshgrid(vr,ibin);
di.rv(:,im)=medianan([d.rv(:,i1).*w2-vr]',iav)'+rvav;
% estimate mean STD of U and V
- di.ruvs(:,im)=sqrt(rus.^2+stdnan([d.rv(:,i1).*w2]')'.^2);
+ di.ruvs(:,im)=sqrt(rus.^2+stdnan([d.rv(:,i1).*w2-vr]')'.^2);
% W
wr=medianan(d.rw(izr,i1).*w);
rwav=meannan(wr);
--- a/process_cast.m Wed Jan 17 12:29:40 2018 -0500
+++ b/process_cast.m Sat Apr 10 08:03:07 2021 -0400
@@ -40,9 +40,9 @@
%======================================================================
% P R O C E S S _ C A S T . M
% doc: Thu Jun 24 16:54:23 2004
-% dlm: Wed Mar 29 12:57:24 2017
+% dlm: Wed Sep 4 17:59:11 2019
% (c) 2004 A.M. Thurnherr
-% uE-Info: 90 39 NIL 0 0 72 0 2 8 NIL ofnI
+% uE-Info: 467 12 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================
% NOTES:
@@ -88,6 +88,12 @@
% Mar 29, 2017: - added att and da to saveres
% - added saveplot_pdf
% - made fignums 2-digit
+% Sep 15, 2018: - disabled serial-number code
+% Feb 8, 2019: - added pause before saving figures (TheThinMint requires this)
+% Feb 16, 2019: - move cast post-processing to step 17 so that post-processing
+% is done before results are saved
+% Aug 30, 2019: - changed error message about p.getdepth
+% Sep 4, 2019: - replaced [getshear2.m] by GK's new [calc_shear3.m]
%----------------------------------------------------------------------
% STEP 0: EXECUTE ALWAYS
@@ -173,7 +179,7 @@
[d,p]=loadrdi(f,p);
% get instrument serial number
- p=getserial(f,p);
+ % p=getserial(f,p);
end_processing_step;
end % OF STEP 1: LOAD DATA
@@ -291,7 +297,7 @@
if p.getdepth==2
[d,p]=getdpthi(d,p);
if length(find(~isfinite(d.izm(1,:))))
- error('Non-finite values in d.izm --- try processing with p.getdepth == 1');
+ error('Missing values in d.izm --- likely missing values in CTD file (processing with p.getdepth = 1; might work');
end
else
[d,p]=getdpth(d,p);
@@ -456,9 +462,9 @@
if ps.shear>0
if ps.shear==1
- [ds,dr,ps,p]=getshear2(d,p,ps,dr);
+ [ds,p,dr]=calc_shear3(d,p,ps,dr); % use all data (d)
else
- [ds,dr,ps,p]=getshear2(di,p,ps,dr);
+ [ds,p,dr]=calc_shear3(di,p,ps,dr); % use superensemble data (di)
end
end
@@ -531,6 +537,10 @@
dr.ctd_ss=interp1q(d.ctdprof_z,d.ctdprof_ss,dr.z);
end
+ if exist('post_process_cast','file') % cruise-specific post-processing
+ post_process_cast;
+ end
+
if length(f.res)>1
%
@@ -550,7 +560,7 @@
disp(sprintf(' figure %d...',j));
ok = 1; eval(sprintf('h = get(%d);',j),'ok = 0;');
if ok
- figure(j);
+ figure(j); pause(1);
eval(sprintf('print -dpsc %s_%02d.ps',f.res,j))
end
end
@@ -560,7 +570,7 @@
if any(ismember(j,pcs.update_figures))
ok = 1; eval(sprintf('h = get(%d);',j),'ok = 0;');
if ok
- figure(j)
+ figure(j); pause(1);
eval(sprintf('print -dpng %s_%02d.png',f.res,j))
end
end
@@ -570,7 +580,7 @@
if any(ismember(j,pcs.update_figures))
ok = 1; eval(sprintf('h = get(%d);',j),'ok = 0;');
if ok
- figure(j)
+ figure(j); pause(1);
eval(sprintf('print -dpdf %s_%02d.pdf',f.res,j))
end
end
@@ -592,10 +602,6 @@
% FINAL STEP: CLEAN UP
%----------------------------------------------------------------------
-if exist('post_process_cast','file') % cruise-specific post-processing
- post_process_cast;
-end
-
fclose('all'); % close all files just to make sure
disp(' ') % final message
--- a/savearch.m Wed Jan 17 12:29:40 2018 -0500
+++ b/savearch.m Sat Apr 10 08:03:07 2021 -0400
@@ -1,14 +1,15 @@
%======================================================================
% S A V E A R C H . M
% doc: Wed Jan 7 16:51:58 2009
-% dlm: Thu Aug 15 12:17:03 2013
+% dlm: Wed Sep 4 16:11:00 2019
% (c) 2009 A.M. Thurnherr
-% uE-Info: 217 16 NIL 0 0 72 0 2 8 NIL ofnI
+% uE-Info: 12 60 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================
% CHANGES BY ANT:
% Jan 7, 2009: - tightened use of exist()
% Aug 15, 2013: - adapted to new [ladcp2cdf.m]
+% Sep 4, 2019: - fixed julian bug reported by Eric Firing
function [da]=savearch(dr,d,p,ps,f,att)
% function [da]=savearch(dr,d,p,ps,f,att)
@@ -19,7 +20,7 @@
p=setdefv(p,'ladcp_cast',1);
g=gregoria(d.time_jul(1));
p=setdefv(p,'ref_year',g(1));
-year0=julian([p.ref_year,0,0,0,0,0]);
+year0=julian([p.ref_year,1,1,0,0,0]);
da.GEN_Velocity_Units = 'm/s';