From 8cbdaaff14f7bd7c5f009d9475ebb100ca4fae1f Mon Sep 17 00:00:00 2001 From: Gael Alory <gael.alory@legos.obs-mip.fr> Date: Thu, 8 Oct 2009 09:53:57 +0000 Subject: [PATCH] bug dans calcul d'ecart aux donnees externes quand seulement 2-3 donnees externes + calcul d'erreur instrument quand pas de donnees de conductivite --- tsg_util/plot_Correction.m | 4 +- tsg_util/plot_Validation.m | 2 +- tsg_util/tsg_accuracy.m | 113 ++++++++++++++++++++----------------- 3 files changed, 63 insertions(+), 56 deletions(-) diff --git a/tsg_util/plot_Correction.m b/tsg_util/plot_Correction.m index 9e766c2..de189c5 100644 --- a/tsg_util/plot_Correction.m +++ b/tsg_util/plot_Correction.m @@ -30,7 +30,7 @@ if ~isempty( tsg.([SAMPLE '_EXT']) ) % Plot bucket samples % ------------------- - ind = 1: length( tsg.([SAMPLE '_EXT_TYPE'])); + ind = 1: length( tsg.([SAMPLE '_EXT'])); indWS = strmatch( 'WS', tsg.([SAMPLE '_EXT_TYPE']), 'exact'); if ~isempty(indWS) plot_Tsg( hMainFig, hPlotAxes, 1, tsg.DAYD_EXT(indWS),... @@ -68,7 +68,7 @@ end % Plot SAMPLE and SSPS (SSJT, SSTP) with code color, on axe 2 % ----------------------------------------------------------- if ~isempty( tsg.([SAMPLE '_EXT']) ) - ind = 1: length( tsg.([SAMPLE '_EXT_TYPE'])); + ind = 1: length( tsg.([SAMPLE '_EXT'])); indWS = strmatch( 'WS', tsg.([SAMPLE '_EXT_TYPE']), 'exact'); if ~isempty(indWS) plot_Tsg( hMainFig, hPlotAxes, 2, tsg.DAYD_EXT(indWS),... diff --git a/tsg_util/plot_Validation.m b/tsg_util/plot_Validation.m index 54f7ff0..24692d7 100644 --- a/tsg_util/plot_Validation.m +++ b/tsg_util/plot_Validation.m @@ -44,7 +44,7 @@ switch nPlot % Plot squares for WS data % ----------------------- - ind = 1: size( tsg.([SAMPLE '_EXT_TYPE']),1); + ind = 1: length( tsg.([SAMPLE '_EXT'])); indWS = strmatch( 'WS', tsg.([SAMPLE '_EXT_TYPE']), 'exact'); if ~isempty(ind) plot_Tsg( hMainFig, hPlotAxes, 1,... diff --git a/tsg_util/tsg_accuracy.m b/tsg_util/tsg_accuracy.m index 8f3af21..181cdf6 100644 --- a/tsg_util/tsg_accuracy.m +++ b/tsg_util/tsg_accuracy.m @@ -52,59 +52,66 @@ if strcmp( PARA, 'SSPS') if ~isempty( tsg.CNDC_CAL) C = tsg.CNDC_CAL(dtTsg); else - C = tsg.CNDC(dtTsg); - end + if ~isempty( tsg.CNDC) + + C = tsg.CNDC(dtTsg); + + % salinity error is computed from temperature/conductivity errors + % by error propagation in the salinity equation + % see Emery and Thomson, Data analysis methods in phys. oceano., p.273 + + R=C/sw_c3515(); + rt = sw_salrt(T); + Rt=R./rt; + + %Rt is a function of C,T + %error on Rt : err(Rt)^2=(d(Rt)/dC*errC)^2+(d(Rt)/dT*errT)^2 + + dRtdC=1./(rt*sw_c3515()); + + c68 = 1.00024; + c0 = 0.6766097; + c1 = 2.00564e-2; + c2 = 1.104259e-4; + c3 = -6.9698e-7; + c4 = 1.0031e-9; + dRtdT=c68*R.*(c1+(2*c2+(3*c3+4*c4*T).*T*c68).*T*c68)./(rt.^2); + + errRt=sqrt((dRtdC*errC).^2+(dRtdT*errT).^2); + + %S is a function of Rt,T + %error on S: err(S)^2=(dS/dRt*errRt)^2+(dS/dT*errT)^2 + + del_T68 = T * 1.00024 - 15; + Rtx = sqrt(Rt); + a0 = 0.0080; + a1 = -0.1692; + a2 = 25.3851; + a3 = 14.0941; + a4 = -7.0261; + a5 = 2.7081; + b0 = 0.0005; + b1 = -0.0056; + b2 = -0.0066; + b3 = -0.0375; + b4 = 0.0636; + b5 = -0.0144; + k = 0.0162; + dSdRt=(a1+(2*a2+(3*a3+(4*a4+5*a5*Rtx).*Rtx).*Rtx).*Rtx)./(2*Rtx)... + +(del_T68./(1+k*del_T68)).*(b1+(2*b2+(3*b3+(4*b4+5*b5*Rtx).*Rtx).*Rtx).*Rtx)./(2*Rtx); + + dSdT=(b0+(b1+(b2+(b3+(b4+b5.*Rtx).*Rtx).*Rtx).*Rtx).*Rtx)./((1+k*del_T68).^2); + + errS=sqrt((dSdRt.*errRt).^2+(dSdT.*errT).^2); + + minerror=errS; + + else + + minerror=0.02; + + end - % salinity error is computed from temperature/conductivity errors - % by error propagation in the salinity equation - % see Emery and Thomson, Data analysis methods in phys. oceano., p.273 - - R=C/sw_c3515(); - rt = sw_salrt(T); - Rt=R./rt; - - %Rt is a function of C,T - %error on Rt : err(Rt)^2=(d(Rt)/dC*errC)^2+(d(Rt)/dT*errT)^2 - - dRtdC=1./(rt*sw_c3515()); - - c68 = 1.00024; - c0 = 0.6766097; - c1 = 2.00564e-2; - c2 = 1.104259e-4; - c3 = -6.9698e-7; - c4 = 1.0031e-9; - dRtdT=c68*R.*(c1+(2*c2+(3*c3+4*c4*T).*T*c68).*T*c68)./(rt.^2); - - errRt=sqrt((dRtdC*errC).^2+(dRtdT*errT).^2); - - %S is a function of Rt,T - %error on S: err(S)^2=(dS/dRt*errRt)^2+(dS/dT*errT)^2 - - del_T68 = T * 1.00024 - 15; - Rtx = sqrt(Rt); - a0 = 0.0080; - a1 = -0.1692; - a2 = 25.3851; - a3 = 14.0941; - a4 = -7.0261; - a5 = 2.7081; - b0 = 0.0005; - b1 = -0.0056; - b2 = -0.0066; - b3 = -0.0375; - b4 = 0.0636; - b5 = -0.0144; - k = 0.0162; - dSdRt=(a1+(2*a2+(3*a3+(4*a4+5*a5*Rtx).*Rtx).*Rtx).*Rtx)./(2*Rtx)... - +(del_T68./(1+k*del_T68)).*(b1+(2*b2+(3*b3+(4*b4+5*b5*Rtx).*Rtx).*Rtx).*Rtx)./(2*Rtx); - - dSdT=(b0+(b1+(b2+(b3+(b4+b5.*Rtx).*Rtx).*Rtx).*Rtx).*Rtx)./((1+k*del_T68).^2); - - errS=sqrt((dSdRt.*errRt).^2+(dSdT.*errT).^2); - - minerror=errS; - -end + end end -- GitLab