0001 function [cal_coeff,cal_stats]=calibration_oxygen(cruiseid)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029 maxiter=10;
0030 tcor = 0.0022;
0031 pcor = 1.350e-4;
0032 Boc=0.0;
0033 Soc=0.4088;
0034 Voffset=-0.4884;
0035
0036
0037
0038
0039 group = group_name(cruiseid);
0040
0041
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
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
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
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;
0078
0079
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
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
0129
0130
0131
0132 staa=input('Input first station: ');
0133 stab=input('Input last station: ');
0134 stddev=input('Enter # of standard deviations for data filtering : ');
0135
0136
0137
0138
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
0171 stn=sta(s_index);
0172 nis=bot(s_index);
0173 o2_bot_num = num_o2(s_index);
0174 prs=cprs(s_index);
0175
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
0180
0181 oxv0=coxv(s_index);
0182 coxy0=coxy(s_index);
0183 boxy0=boxy(s_index);
0184 boxysav=boxy(s_index);
0185 [tot,tot2]=size(boxy0);
0186
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
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
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
0208
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
0227
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
0235 ylabel('Bottle - CTD Oxygen (ml/l)');
0236 xlabel('Sample Bottle Number');
0237 hold on
0238
0239
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
0250
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
0258 ylabel('Bottle - CTD Oxygen (ml/l)');
0259 xlabel('Pressure');
0260 hold on
0261
0262
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
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
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
0305
0306 d=boxy0-new_oxy;
0307 m=mean(d);
0308 m_abs=mean(abs(d)./boxy0)*100;
0309 s=std(d);
0310
0311 sub=length(boxy0);
0312
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
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
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
0420 sub=length(boxy0(ii));
0421
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
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
0480
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
0488 ylabel('Bottle - CTD Oxygen (ml/l)');
0489 xlabel('Pressure');
0490 hold on
0491 end
0492 return