Home > ctdcal > calibration_oxygen.m

calibration_oxygen

PURPOSE ^

calibration_cond - Main oxygen calibration calling routine.

SYNOPSIS ^

function [cal_coeff,cal_stats]=calibration_oxygen(cruiseid)

DESCRIPTION ^

 calibration_cond - Main oxygen calibration calling routine.

 CTD Calibration toolbox.

 USAGE:
 [cal_coeff,cal_stats]=calibration_oxy(cruiseid);

 INPUT: cruiseid: Cruise identification name.
   
 OUTPUT: .... 
          
 DESCRIPTION: ...  
 
 AUTHORS:
            Carlos Fonseca 
         UM/CIMAS
         Derrick Snowden
         NOAA/AOML/PhOD
         Tue May 11 11:23:37 EDT 2004

 CHANGELOG: 
   23-Jul-2004, Version 1.0
        * Initial version. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%Parameters**************************
H_COLOR=1000; %depth for color coding the output

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [cal_coeff,cal_stats]=calibration_oxygen(cruiseid)
0002 % calibration_cond - Main oxygen calibration calling routine.
0003 %
0004 % CTD Calibration toolbox.
0005 %
0006 % USAGE:
0007 % [cal_coeff,cal_stats]=calibration_oxy(cruiseid);
0008 %
0009 % INPUT: cruiseid: Cruise identification name.
0010 %
0011 % OUTPUT: ....
0012 %
0013 % DESCRIPTION: ...
0014 %
0015 % AUTHORS:
0016 %            Carlos Fonseca
0017 %         UM/CIMAS
0018 %         Derrick Snowden
0019 %         NOAA/AOML/PhOD
0020 %         Tue May 11 11:23:37 EDT 2004
0021 %
0022 % CHANGELOG:
0023 %   23-Jul-2004, Version 1.0
0024 %        * Initial version.
0025 %
0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0027 %%%%%%%%%%%%Parameters**************************
0028 %H_COLOR=1000; %depth for color coding the output
0029 maxiter=10; %  maximum number of iterations
0030 tcor = 0.0022; %nominal value - 12/17/01 sensor
0031 pcor = 1.350e-4; %nominal value - 12/17/01 sensor
0032 Boc=0.0;%nominal value - 12/17/01 sensor
0033 Soc=0.4088;%nominal value - 12/17/01 sensor
0034 Voffset=-0.4884;%nominal value - 12/17/01 sensor
0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0037 
0038 % Retrieve the necessary preferences for this cruise
0039 group = group_name(cruiseid);
0040 
0041 % cruise_dir is the main calibration data directory
0042 
0043 if ispref(group,'cal_data_dir')
0044     cruise_dir=getpref(group,'cal_data_dir');
0045 else
0046     error(['The cal_data_dir preference was not set for cruise: ' cruiseid '.  Run register_cruise.m']);
0047 end
0048 
0049 % Load data from the cruise database
0050 db_file = fullfile(cruise_dir,[cruiseid '_db.mat']);
0051 if exist(db_file)==2
0052     load(db_file);
0053 else
0054     error(['The cruise data base file was not found on the path.  ' db_file]);
0055 end
0056 n=who('oxy_cal');
0057 
0058 if isempty(n)==1
0059     error(['oxy_cal not found. Run make_oxy_sumfile.m ']);
0060     return;
0061 end
0062 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0063 % coefficients for the initialization of the station fit
0064 cal=0;
0065 gama=0;
0066 p1=0;
0067 p2=0;
0068 p3=0;
0069 p4=0;
0070 p5=0;
0071 p6=0;
0072 %%%%%Seting the Variables%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0073 sta=oxy_cal.stnnbr;
0074 bot=oxy_cal.niskinbtlpos;
0075 num_o2=oxy_cal.sampnbr;
0076 cprs=oxy_cal.ctdprs;
0077 boxy=oxy_cal.btloxy;  % if O2 values already in ml/l
0078 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0079 % Choose the sensor to calibrate pri/sec
0080 disp(' 0 = use the primary sensor (DEFAULT)' )
0081 disp(' 1 = use the secondary sensor' )
0082 pri_sec=input('Do you want to use the primary or the secondary sensor: ');
0083 while (isempty(pri_sec)==1 | pri_sec<0 | pri_sec>1)
0084     disp('Please answer 0 (zero) for no or  1(one) for yes ')
0085     pri_sec=input('Do you want to use the primary or the secondary sensor: ');
0086 end
0087 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0088 %%%%%%%%%%Seting other variables%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0089  if (pri_sec==0)
0090  ctem=oxy_cal.ctdtmp1;
0091  csal=oxy_cal.ctdsal1;
0092  coxv=oxy_cal.ctdoxv1;
0093  coxy=oxy_cal.ctdoxy1;
0094  ccon=oxy_cal.ctdcon1;
0095  else
0096  ctem=oxy_cal.ctdtmp2;
0097  csal=oxy_cal.ctdsal2;
0098  coxv=oxy_cal.ctdoxv2;
0099  coxy=oxy_cal.ctdoxy2;
0100  ccon=oxy_cal.ctdcon2;
0101  end
0102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0103 min_sta=min(sta);
0104 max_sta=max(sta);
0105 disp(' ')
0106 disp('PROCESSING CRUISE')
0107 disp('****************************************** ')
0108 disp(' This program needs the Seawater toolbox')
0109 disp('The data found in oxy_cal goes from') 
0110 disp(['the station ',num2str(min_sta),' to the station ',num2str(max_sta),'.'] )
0111 disp(['The maximum pressure is: ',num2str(max(oxy_cal.ctdprs))])
0112 disp('**************************************************************')
0113 %
0114 disp('You have an option to color code the output selecting a depth')
0115 col_code=input('Do you want select a depth? (0=no 1=yes)');
0116 
0117 while (isempty(col_code)==1 | col_code<0 | col_code>1)
0118     disp('Please answer 0 (zero) for no or  1(one) for yes')
0119     col_code=input('Do you want select a depth? (0=no 1=yes)');
0120 end
0121 if col_code==0
0122     H_COLOR=0;
0123 end
0124 
0125 if col_code==1
0126     H_COLOR=input('Enter the depth for color coding the output:');
0127 end
0128 %Careful! If O2 bottle values in the master bottle file are in uM/kg
0129 %they must first be converted to ml/l, so be sure to choose the
0130 %correct line to comment out below (look where boxy0 and coxy0 are assigned)
0131 
0132 staa=input('Input first station: ');
0133 stab=input('Input last station: ');
0134 stddev=input('Enter # of standard deviations for data filtering : ');
0135 %good choice for this is 2.5 for final fit.  Use a higher number
0136 %to see most of the cal points and do a quick initial fit (say, 5).
0137 %lower than 2.5 takes some time to do the fit and rejects a lot
0138 %of points.  Shouldn't use less than about 2 for final fit.
0139 counter=1;
0140 flag_bottle=0;
0141 bad_btl(counter)=input('Enter # of bottle sample to remove or 0 to skip: ');
0142  while bad_btl(counter)~=0
0143  counter=counter+1;
0144  bad_btl(counter)=input('Enter oxygen bottle sample number to remove or 0 to skip: ');
0145  end
0146  if (counter > 1)
0147  flag_bottle=1;
0148  bad_bottle=bad_btl(1:length(bad_btl)-1);  
0149  disp('Sample Bottles marked as bad')
0150  bad_bottle
0151  end
0152 disp('Fit as a function of station number?: ')
0153 disp('  0 = none no fit for the station number (DEFAULT)')
0154 disp('  1 = exp fits an extra term [cal*exp(gama*sta)] ')
0155 disp('  2 = First Order:  linear fits [p1*sta] ')
0156 disp('  3 = Second Order: poly2 fits  [p1*sta+p2*sta^2] ')
0157 disp('  4 = Third Order:  poly3 fits  [p1*sta+p2*sta^2+p3*sta^3] ')
0158 disp('  5 = Four Order:   poly4 fits  [p1*sta+p2*sta^2+p3*sta^3+p4*sta^4] ')
0159 disp('  6 = Fifth Order:  poly5 fits  [p1*sta+p2*sta^2+p3*sta^3+p4*sta^4+p5*sta^5] ')
0160 disp('  7 = Sixth Order:  poly6 fits  [p1*sta+p2*sta^2+p3*sta^3+p4*sta^4+p5*sta^5+p6*sta^6] ')
0161 method=input('What kind of fit for the station number: ');
0162 
0163 if (isempty(method)==1);
0164     method=0;
0165 end
0166     stnrange=[staa stab];
0167 
0168 s_index=find(sta >= stnrange(1) & sta  <= stnrange(2));
0169 
0170 % subsample for selected station numbers
0171 stn=sta(s_index);
0172 nis=bot(s_index);
0173 o2_bot_num = num_o2(s_index);
0174 prs=cprs(s_index);
0175 %tem0 = ctem(s_index);
0176 tem0=sw_ptmp(csal(s_index),ctem(s_index),cprs(s_index),0);
0177 ccon0=ccon(s_index);
0178 csal0=csal(s_index);
0179 %oxc0=coxc2(s_index);
0180 %oxt0=coxt2(s_index);
0181 oxv0=coxv(s_index);
0182 coxy0=coxy(s_index);  % if O2 values already in ml/l
0183 boxy0=boxy(s_index);  % if O2 values already in ml/l
0184 boxysav=boxy(s_index);
0185 [tot,tot2]=size(boxy0);
0186 % removing the sample bottles marked as bad bottles
0187 if (flag_bottle==1);
0188     for b=1:length(bad_bottle)
0189     bad_index=find(o2_bot_num==bad_bottle(b))
0190     boxy0(bad_index)=NaN;
0191     end
0192 end
0193 % saving the rejected bottles
0194 tt=find(isnan(boxy0)==1);
0195 trash=[stn(tt) nis(tt) prs(tt) o2_bot_num(tt)] 
0196 save trash trash -ascii
0197 
0198 % get rid of values where bottle O2 = NaN
0199 jj_nan=find(isnan(boxy0)==0);
0200 stn=stn(jj_nan);
0201 nis=nis(jj_nan);
0202 o2_bot_num=o2_bot_num(jj_nan);
0203 prs=prs(jj_nan);
0204 tem0 = tem0(jj_nan);
0205 ccon0=ccon0(jj_nan);
0206 csal0=csal0(jj_nan);
0207 %oxc0=oxc0(jj_nan);
0208 %oxt0=oxt0(jj_nan);
0209 oxv0=oxv0(jj_nan);
0210 coxy0=coxy0(jj_nan);  
0211 boxy0=boxy0(jj_nan);  
0212 boxysav=boxysav(jj_nan);
0213 ii=tot;
0214 
0215 figure;
0216 subplot(2,1,1);
0217 deep_idx=find(prs>=H_COLOR);
0218 shallow_idx=find(prs<H_COLOR);
0219 plot(nis(deep_idx),boxy0(deep_idx)-coxy0(deep_idx),'b+');
0220 hold on
0221 plot(nis(shallow_idx),boxy0(shallow_idx)-coxy0(shallow_idx),'r+');
0222 grid on
0223 title('Dissolved Oxygen Calibration - raw data');
0224 ylabel('Bottle - Primary CTD Oxygen (ml/l)');
0225 xlabel('Niskin Bottle');
0226 %set(gca,'YLim', [-0.4 0.4]);
0227 %set(gca,'Ytick', [-0.4 -0.2 0 0.2 0.4]);
0228 hold on
0229 subplot(2,1,2);
0230 plot(o2_bot_num(deep_idx),boxy0(deep_idx)-coxy0(deep_idx),'b+');
0231 hold on
0232 plot(o2_bot_num(shallow_idx),boxy0(shallow_idx)-coxy0(shallow_idx),'r+');
0233 grid on
0234 %title('Bottle-Secondary CTD Oxygen vs. pressure - all stations');
0235 ylabel('Bottle - CTD Oxygen (ml/l)');
0236 xlabel('Sample Bottle Number');
0237 hold on
0238 %set(gca,'YLim', [-0.4 0.4]);
0239 %set(gca,'Ytick', [-0.4 -0.2 0 0.2 0.4]);
0240 figure
0241 subplot(2,1,1);
0242 plot(stn(deep_idx),boxy0(deep_idx)-coxy0(deep_idx),'b+');
0243 hold on
0244 plot(stn(shallow_idx),boxy0(shallow_idx)-coxy0(shallow_idx),'r+');
0245 grid on
0246 title('Dissolved Oxygen Calibration - raw data');
0247 ylabel('Bottle - Primary CTD Oxygen (ml/l)');
0248 xlabel('Station number');
0249 %set(gca,'YLim', [-0.4 0.4]);
0250 %set(gca,'Ytick', [-0.4 -0.2 0 0.2 0.4]);
0251 hold on
0252 subplot(2,1,2);
0253 plot(prs(deep_idx),boxy0(deep_idx)-coxy0(deep_idx),'b+');
0254 hold on
0255 plot(prs(shallow_idx),boxy0(shallow_idx)-coxy0(shallow_idx),'r+');
0256 grid on
0257 %title('Bottle-Secondary CTD Oxygen vs. pressure - all stations');
0258 ylabel('Bottle - CTD Oxygen (ml/l)');
0259 xlabel('Pressure');
0260 hold on
0261 
0262 %constants for oxsat calculation
0263 A1 = -173.4292;
0264 A2 = 249.6339;
0265 A3 = 143.3438;
0266 A4 = -21.8492;
0267 B1 = -0.033096;
0268 B2 = 0.014259;
0269 B3 = -0.00170;
0270 T  = tem0 + 273.15;
0271 
0272 disp(' ');
0273 disp('First guess SOC & Voffset (linear regression):');
0274 
0275 [msize,nsize] = size(boxy0);
0276 oxy=zeros(msize,nsize);
0277 oxsat=zeros(msize,nsize);
0278 oxsat = exp(A1 + A2*(100./T) + A3*log(T./100) + A4*(T./100) + csal0 .* (B1+B2*(T./100)+ B3*(T./100).*(T./100)));
0279 phi=exp(tcor.*tem0).*oxsat.*exp(pcor.*prs);
0280 oxones=1+0*oxy;
0281 c=[oxv0 oxones]\(boxy0./phi);
0282 Soc=c(1);
0283 Voffset=c(2)/c(1);
0284 new_oxy=Soc.*(oxv0+Voffset).*phi;
0285 % Evaluating the fit
0286 figure;
0287 
0288 plot(oxv0(deep_idx),boxy0(deep_idx)./phi(deep_idx),'b*');
0289 hold on
0290 plot(oxv0(shallow_idx),boxy0(shallow_idx)./phi(shallow_idx),'r*');
0291 hold on
0292 
0293 %text(oxv0,boxy0./phi,num2str(stn));
0294 grid on
0295 hold on
0296 plot(oxv0(deep_idx),new_oxy(deep_idx)./phi(deep_idx),'bo')
0297 hold on
0298 plot(oxv0(shallow_idx),new_oxy(shallow_idx)./phi(shallow_idx),'ro')
0299 hold on
0300 title('First Guess: Soc and Voffset only ');
0301 ylabel('Bottle Oxygen (ml/l)/phi');
0302 xlabel('SBE voltage');
0303 legend('raw deep','raw shallow','fitted depp','fitted shallow') 
0304 %set(gca,'YLim', [-0.4 0.4]);
0305 %set(gca,'Ytick', [-0.4 -0.2 0 0.2 0.4]);
0306 d=boxy0-new_oxy;
0307 m=mean(d);
0308 m_abs=mean(abs(d)./boxy0)*100;
0309 s=std(d);
0310 % display some statistics
0311 sub=length(boxy0);
0312 %disp(['Forming initial estimates for Soc and Boc'])
0313 disp([' Primary Sensor']);
0314 disp(['Number of points used   ',int2str(sub)]);
0315 disp(['Total number of points  ',int2str(tot)]);
0316 disp(['Percent of points used  ',num2str(100*sub/tot,4)]);
0317 disp(['fit standard deviation  ',num2str(s,4)]);
0318 disp(['fit mean residual  ',num2str(m,4)]);
0319 disp(['fit mean absolute residual  ',num2str(mean(abs(d)),4)]);
0320 disp(['fit mean absolute residual (%) ',num2str(mean(abs(d)./boxy0)*100,4)]);
0321 disp(['First guess Soc = ', num2str(Soc)]);
0322 disp(['First guess Voffset = ', num2str(Voffset)]);
0323 m_old=m;
0324 m_abs_old=m_abs;
0325 s_old=s;
0326 d_old=d;
0327 %Soc= boxy0 ./(oxsat .* exp(tcor*(tem0+wt*(oxt0-tem0))+pcor*prs).*(oxc0 - zoc));
0328 %Boc= -Soc*zoc;
0329 %for i=1:tot
0330 %  oxsat(i) = exp(A1 + A2*(100/T(i)) + A3*log(T(i)/100) + A4*(T(i)/100) + csal0(i)*(B1+B2*(T(i)/100)+ B3*(T(i)/100)*(T(i)/100)));
0331 %  Soc(i)= boxy0(i)/(oxsat(i)*exp(tcor*(tem0(i)+wt*(oxt0(i)-tem0(i)))+pcor*prs(i))*(oxc0(i) - zoc));
0332 %  Boc(i)= -Soc(i)*zoc;
0333 %end
0334 
0335 %check for outliers and get a best estimate for Soc
0336 %stddev should be set to about 2 or 3.
0337 %jj=find(abs(d)>stddev*s);
0338 %boxy0(jj)=NaN*jj;
0339 jj=1;
0340 
0341 iter=0;
0342 iold=1e29;
0343 %%%%
0344 new_oxy_old=new_oxy;
0345 ii_old=find(finite(boxy0)==1);
0346 Soc_old=Soc ;
0347 Voffset_old=Voffset;
0348 tcor_old = tcor; 
0349 pcor_old = pcor; 
0350 Boc_old = Boc;
0351 cal_old = cal;
0352 gama_old = gama;
0353 p1_old=p1;
0354 p2_old=p2;
0355 p3_old=p3;
0356 p4_old=p4;
0357 p5_old=p5;
0358 p6_old=p6;
0359 
0360 %%%
0361 while length(jj)>0;
0362    iter=iter+1;
0363    if iter>maxiter
0364        final_results
0365        break
0366    end
0367    
0368    
0369 disp('Looping through fit and rejecting outliers...'); 
0370 disp(' ');
0371 
0372 ii=find(finite(boxy0)==1);
0373 if (method==0);
0374 [new_oxy,Soc,Voffset,Boc,tcor,pcor]=...
0375     oxy_fit_none(oxv0(ii),oxsat(ii),tem0(ii),prs(ii),stn(ii),boxy0(ii),Soc,Voffset,Boc,tcor,pcor);
0376 end
0377 if (method==1);
0378    [new_oxy,Soc,Voffset,Boc,tcor,pcor,cal,gama]=...
0379     oxy_fit_exp(oxv0(ii),oxsat(ii),tem0(ii),prs(ii),stn(ii),boxy0(ii),Soc,Voffset,Boc,tcor,pcor,cal,gama);
0380 end
0381 if (method==2);
0382      [new_oxy,Soc,Voffset,Boc,tcor,pcor,p1]=...
0383     oxy_fit_linear(oxv0(ii),oxsat(ii),tem0(ii),prs(ii),stn(ii),boxy0(ii),Soc,Voffset,Boc,tcor,pcor,p1);
0384 end
0385 if (method==3);
0386         [new_oxy,Soc,Voffset,Boc,tcor,pcor,p1,p2]=...
0387     oxy_fit_poly2(oxv0(ii),oxsat(ii),tem0(ii),prs(ii),stn(ii),boxy0(ii),Soc,Voffset,Boc,tcor,pcor,p1,p2);
0388 end        
0389 if (method==4);
0390      [new_oxy,Soc,Voffset,Boc,tcor,pcor,p1,p2,p3]=...
0391     oxy_fit_poly3(oxv0(ii),oxsat(ii),tem0(ii),prs(ii),stn(ii),boxy0(ii),Soc,Voffset,Boc,tcor,pcor,p1,p2,p3);
0392 end
0393 if (method==5);
0394      [new_oxy,Soc,Voffset,Boc,tcor,pcor,p1,p2,p3,p4]=...
0395     oxy_fit_poly4(oxv0(ii),oxsat(ii),tem0(ii),prs(ii),stn(ii),boxy0(ii),Soc,Voffset,Boc,tcor,pcor,p1,p2,p3,p4);
0396 end
0397 if (method==6);
0398      [new_oxy,Soc,Voffset,Boc,tcor,pcor,p1,p2,p3,p4,p5]=...
0399     oxy_fit_poly5(oxv0(ii),oxsat(ii),tem0(ii),prs(ii),stn(ii),boxy0(ii),Soc,Voffset,Boc,tcor,pcor,p1,p2,p3,p4,p5);
0400 end
0401 if (method==7);
0402      [new_oxy,Soc,Voffset,Boc,tcor,pcor,p1,p2,p3,p4,p5,p6]=...
0403     oxy_fit_poly6(oxv0(ii),oxsat(ii),tem0(ii),prs(ii),stn(ii),boxy0(ii),Soc,Voffset,Boc,tcor,pcor,p1,p2,p3,p4,p5,p6);
0404 end
0405 d=boxy0(ii)-new_oxy; 
0406 m=mean(d);
0407 m_abs=mean(abs(d)./boxy0(ii))*100;
0408 s=std(d);
0409 if m_abs>m_abs_old
0410        disp('The new mean is bigger...'); 
0411        final_results
0412        break
0413 elseif s>s_old    
0414     disp('The new std is bigger...');
0415     final_results
0416     break    
0417 end
0418 
0419 % display some statistics
0420 sub=length(boxy0(ii));
0421 %disp(['Forming initial estimates for Soc and Boc'])
0422 disp([' Primary Sensor']);
0423 disp(['Number of points used   ',int2str(sub)]);
0424 disp(['Total number of points  ',int2str(tot)]);
0425 disp(['Percent of points used  ',num2str(100*sub/tot,4)]);
0426 disp(['fit standard deviation  ',num2str(s,4)]);
0427 disp(['fit mean residual  ',num2str(m,4)]);
0428 disp(['fit mean absolute residual  ',num2str(mean(abs(d)),4)]);
0429 disp(['fit mean absolute residual (%)  ',num2str(mean(abs(d)./boxy0(ii))*100,4)]);
0430 disp(['First guess Soc = ', num2str(Soc)]);
0431 disp(['First guess Voffset = ', num2str(Voffset)]);
0432 disp(['First guess Boc = ', num2str(Boc)]);
0433 disp(['First guess tcor = ', num2str(tcor)]);
0434 disp(['First guess pcor = ', num2str(pcor)]);
0435 disp(['First guess cal = ', num2str(cal)]);
0436 disp(['First guess gama = ', num2str(gama)]);
0437 disp(['First guess p1 = ', num2str(p1)]);
0438 disp(['First guess p2 = ', num2str(p2)]);
0439 disp(['First guess p3 = ', num2str(p3)]);
0440 disp(['First guess p4 = ', num2str(p4)]);
0441 disp(['First guess p5 = ', num2str(p5)]);
0442 disp(['First guess p6 = ', num2str(p6)]);
0443 rejects=[stn(jj) nis(jj) prs(jj) o2_bot_num(jj)] 
0444 eval(['save rejs' num2str(iter) ' rejects -ascii'])
0445 jj=find(abs(d)>stddev*s);
0446 boxy0(jj)=NaN*jj; 
0447 % refreshing the variables for the next iteration
0448 m_old=m;
0449 m_abs_old=m_abs;
0450 s_old=s;
0451 d_old=d;
0452 new_oxy_old=new_oxy;
0453 ii_old=ii;
0454 Soc_old=Soc ;
0455 Voffset_old=Voffset;
0456 tcor_old = tcor; 
0457 pcor_old = pcor; 
0458 Boc_old = Boc;
0459 cal_old = cal;
0460 gama_old = gama;
0461 p1_old=p1;
0462 p2_old=p2;
0463 p3_old=p3;
0464 p4_old=p4;
0465 p5_old=p5;
0466 p6_old=p6;
0467 %
0468 figure;
0469 subplot(2,1,1);
0470 deep_idx=find(prs(ii)>=H_COLOR);
0471 shallow_idx=find(prs(ii)<H_COLOR);
0472 plot(stn(ii(deep_idx)),d(deep_idx),'bo');
0473 hold on
0474 plot(stn(ii(shallow_idx)),d(shallow_idx),'ro');
0475 grid on
0476 title(['Dissolved Oxygen Calibration - loop: ',num2str(iter)]);
0477 ylabel('Bottle - Primary CTD Oxygen (ml/l)');
0478 xlabel('Station number');
0479 %set(gca,'YLim', [-0.4 0.4]);
0480 %set(gca,'Ytick', [-0.4 -0.2 0 0.2 0.4]);
0481 hold on
0482 subplot(2,1,2);
0483 plot(prs(ii(deep_idx)),d(deep_idx),'bo');
0484 hold on
0485 plot(prs(ii(shallow_idx)),d(shallow_idx),'ro');
0486 grid on
0487 %title('Bottle-Secondary CTD Oxygen vs. pressure - all stations');
0488 ylabel('Bottle - CTD Oxygen (ml/l)');
0489 xlabel('Pressure');
0490 hold on
0491 end
0492 return

Generated on Fri 08-Oct-2004 11:57:17 by m2html © 2003