mcoxyhyst_reverse that handles missing data

source/mextras/mcoxyhyst_reverse.m

-function [oxygen_rev]=mcoxyhyst_reverse(oxygen_sbe,time,press,H1,H2,H3)
-% function [oxygen_rev]=mcoxyhyst_reverse(oxygen_sbe,time,press,H1,H2,H3)
-%
-% gdm on di346
-% function to reverse the adjustment for hysteresis in the oxygen sensor
-% applied by SeaBird 
-% 
-% the default values and ranges or the constants are:
-% H1     -0.033     [-0.02 to -0.05]
-% H2      5000      
-% H3      1450      [1200  to  2000]
-
-oxygen_rev=oxygen_sbe;
-
-for k=2:length(time)
-    D=1+H1*(exp(press(k)/H2)-1);
-    C=exp(-1*(time(k)-time(k-1))/H3);
-
-    oxygen_rev(k)=D*(oxygen_sbe(k)-C*oxygen_sbe(k-1))+C*oxygen_rev(k-1);
-    
+function [oxygen_rev]=mcoxyhyst_reverse(oxygen_sbe,time,press,H1,H2,H3)
+% function [oxygen_rev]=mcoxyhyst_reverse(oxygen_sbe,time,press,H1,H2,H3)
+%
+% gdm on di346
+% function to reverse the adjustment for hysteresis in the oxygen sensor
+% applied by SeaBird 
+% 
+% the default values and ranges or the constants are:
+% H1     -0.033     [-0.02 to -0.05]
+% H2      5000      
+% H3      1450      [1200  to  2000]
+
+% DAS chnaged to deal with absent data but note this could give slightly different results from original data as 
+% the sbe software may deal with absent data differently
+
+oxygen_rev=oxygen_sbe;
+
+kfirst = min(find(isfinite(oxygen_sbe)));
+klastgood = kfirst; % keep track of most recent good cycle
+
+for k=2:length(time)
+    D=1+H1*(exp(press(k)/H2)-1);
+    C=exp(-1*(time(k)-time(k-1))/H3);
+
+    oxygen_rev(k)=D*(oxygen_sbe(k)-C*oxygen_sbe(k-1))+C*oxygen_rev(k-1);
+    if isnan(oxygen_sbe(k)+press(k))
+        oxygen_rev(k) = nan; %already the case because of initialisation of oxygen_out
+    else
+        if press(k) < 0; press(k) = 0; end
+        D=1+H1*(exp(press(k)/H2)-1);
+        C=exp(-1*(time(k)-time(klastgood))/H3);
+		oxygen_rev(k)=D*(oxygen_sbe(k)-C*oxygen_sbe(klastgood))+C*oxygen_rev(klastgood);
+        klastgood = k;
+    end    
 end;
\ No newline at end of file