Description of calibration

0001 function [cal_coeff,cal_stats]=calibration_cond(cruiseid)
0002 % calibration_cond - Main conductivity calibration calling routine.
0003 %
0004 % CTD Calibration toolbox.
0005 %
0006 % USAGE:
0007 % [cal_coeff,cal_stats]=calibration_cond(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 
0028 % Retrieve the necessary preferences for this cruise
0029 group = group_name(cruiseid);
0030 
0031 % cruise_dir is the main calibration data directory
0032 
0033 if ispref(group,'cal_data_dir')
0034     cruise_dir=getpref(group,'cal_data_dir');
0035 else
0036     error(['The cal_data_dir preference was not set for cruise: ' cruiseid '.  Run register_cruise.m']);
0037 end
0038 
0039 % Load data from the cruise database
0040 db_file = fullfile(cruise_dir,[cruiseid '_db.mat']);
0041 if exist(db_file)==2
0042     load(db_file);
0043 else
0044     error(['The cruise data base file was not found on the path.  ' db_file]);
0045 end
0046 n=who('cond_cal');
0047 
0048 if isempty(n)==1
0049     error(['cond_cal not found. Run make_cond_sumfile.m ']);
0050     return;
0051 end
0052 
0053 stations=unique(cond_cal.stnnbr);
0054 station_min=min(stations);
0055 station_max=max(stations);
0056 disp('**************************************************************')
0057 disp('The data found in cond_cal goes from') 
0058 disp(['the station ',num2str(station_min),' to the station ',num2str(station_max),'.'] )
0059 disp(['The maximum pressure is: ',num2str(max(cond_cal.ctdprs))])
0060 disp('**************************************************************')
0061 %
0062 type_fit={'calcos0' 'calcos1' 'calcos2' 'calcos3' 'calcos4' 'calcos5' 'calcos6'...
0063           'calcop0' 'calcop1' 'calcop2' 'calcop3' 'calcop4' 'calcop5' 'calcop6'...
0064           'calcodcdp0' 'calcodcdp1' 'calcodcdp2' 'calcodcdp3' 'calcodcdp4' 'calcodcdp5' 'calcodcdp6'...
0065           'calcopdcdp0' 'calcopdcdp1' 'calcopdcdp2' 'calcopdcdp3' 'calcopdcdp4' 'calcopdcdp5' 'calcopdcdp6'};
0066 
0067 % Choose the sensor to calibrate pri/sec
0068 disp(' 0 = use the primary sensor (DEFAULT)' )
0069 disp(' 1 = use the secondary sensor' )
0070 pri_sec=input('Do you want to use the primary or the secondary sensor: ');
0071 while (isempty(pri_sec)==1 | pri_sec<0 | pri_sec>1)
0072     disp('Please answer 0 (zero) for no or  1(one) for yes ')
0073     pri_sec=input('Do you want to use the primary or the secondary sensor: ');
0074 end
0075 
0076 % Station interval
0077 st1=input('Enter first station number: ');
0078 st2=input('Enter last station number: ');
0079 
0080 stddev=input('Enter number of standard deviations for fit: ');
0081 
0082 %
0083 disp('You have an option to color code the output selecting a depth')
0084 disp('This depth will be also used in the criteria to throw out data' )
0085 col_code=input('Do you want select a depth? (0=no 1=yes)');
0086 
0087 while (isempty(col_code)==1 | col_code<0 | col_code>1)
0088     disp('Please answer 0 (zero) for no or  1(one) for yes')
0089     col_code=input('Do you want select a depth? (0=no 1=yes)');
0090 end
0091 if col_code==0
0092     hcolor=0;
0093 end
0094 
0095 if col_code==1
0096     hcolor=input('Enter the depth for color coding the output:');
0097 end
0098 
0099 %
0100 deep_all=input('Do you want to make the fit only using deeper values? 0=no 1=yes ');
0101 while (isempty(deep_all)==1 | deep_all<0 | deep_all>1)
0102     disp('Please answer 0 (zero) for no or  1(one) for yes ')
0103     deep_all=input('Do you want to make the fit only using deeper values? 0=no 1=yes ');
0104 end
0105 
0106 %
0107 fig_save=input('Do you want to keep all figures on screen? 0=no 1=yes ');
0108 while (isempty(deep_all)==1 | deep_all<0 | deep_all>1)
0109     disp('Please answer 0 (zero) for no or  1(one) for yes ')
0110     fig_save=input('Do you want to keep all figures on screen? 0=no 1=yes ');
0111 end
0112 
0113 fig=1;
0114 
0115 % Loading the possible models for the conductivity calibration
0116 load_model_types
0117 
0118 % seting up some variables
0119 idx_stats=zeros(2,2);
0120 cal_coeff=zeros(28,11);
0121 cal_stats=zeros(28,4);
0122 cal_stats_deep=zeros(28,4);
0123 fits=1:28;
0124 cal_coeff(:,1)=fits';
0125 cal_stats(:,1)=fits';
0126 cal_stats_deep(:,1)=fits';
0127 counter_fit=1;
0128 max_idx=0;
0129 min_idx=1:10000;
0130 
0131 if pri_sec==0
0132     ii=find(cond_cal.stnnbr>=st1 & cond_cal.stnnbr<=st2);
0133     plottitle='Primary';
0134     ctdcon=cond_cal.ctdcon1;
0135     ctddcdp=cond_cal.dcdp1;
0136     btlcon=cond_cal.btlcon1;
0137     for i=1:28
0138      if (i<8)
0139         ldx=['X=s',num2str(i-1),'(cond_cal.ctdcon1,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.dcdp1);'];
0140         fit=['s',num2str(i-1)];
0141      elseif(i>7 &i<15)
0142         ldx=['X=p',num2str(i-8),'(cond_cal.ctdcon1,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.dcdp1);'];  
0143         fit=['p',num2str(i-8)];
0144      elseif(i>14 & i<22)
0145         ldx=['X=dcdp',num2str(i-15),'(cond_cal.ctdcon1,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.dcdp1);'];
0146         fit=['dcdp',num2str(i-15)];
0147      else
0148         ldx=['X=pdcdp',num2str(i-22),'(cond_cal.ctdcon1,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.dcdp1);'];  
0149         fit=['pdcdp',num2str(i-22)];
0150      end
0151      eval(ldx);
0152      %
0153      [coeff,stats,stats_deep,newbotco,newctdco]=...
0154      make_fit_cond(X,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.btlcon1,stddev,fit,plottitle,hcolor,deep_all,fig);
0155     
0156      [cal_coeff,cal_stats,cal_stats_deep,max_idx,min_idx,idx_stats]=...
0157      save_fit(coeff,stats,stats_deep,newbotco,cal_coeff,cal_stats,cal_stats_deep,max_idx,min_idx,idx_stats,counter_fit);
0158      counter_fit=counter_fit+1;
0159       disp('Hit return to continue...')
0160       pause
0161     
0162      if fig_save==1
0163         fig=fig+1;
0164      end%
0165      %
0166     end
0167 else
0168  ii=find(cond_cal.stnnbr>=st1 & cond_cal.stnnbr<=st2);
0169  plottitle='Secondary';
0170  ctdcon=cond_cal.ctdcon2;
0171  ctddcdp=cond_cal.dcdp2;
0172  btlcon=cond_cal.btlcon2;
0173  for i=1:28
0174      if (i<8)
0175      ldx=['X=s',num2str(i-1),'(cond_cal.ctdcon2,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.dcdp2);'];
0176      fit=['s',num2str(i-1)];
0177      elseif(i>7 &i<15)
0178      ldx=['X=p',num2str(i-8),'(cond_cal.ctdcon2,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.dcdp2);'];  
0179      fit=['p',num2str(i-8)];
0180      elseif(i>14 & i<22)
0181      ldx=['X=dcdp',num2str(i-15),'(cond_cal.ctdcon2,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.dcdp2);'];
0182      fit=['dcdp',num2str(i-15)];
0183      else
0184      ldx=['X=pdcdp',num2str(i-22),'(cond_cal.ctdcon2,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.dcdp2);'];  
0185      fit=['pdcdp',num2str(i-22)];
0186      end
0187      eval(ldx)
0188      %
0189      [coeff,stats,stats_deep,newbotco,newctdco]=...
0190      make_fit_cond(X,cond_cal.stnnbr,cond_cal.ctdprs,cond_cal.btlcon2,stddev,fit,plottitle,hcolor,deep_all,fig);
0191      %
0192      [cal_coeff,cal_stats,cal_stats_deep,max_idx,min_idx,idx_stats]=...
0193      save_fit(coeff,stats,stats_deep,newbotco,cal_coeff,cal_stats,cal_stats_deep,max_idx,min_idx,idx_stats,counter_fit);
0194      counter_fit=counter_fit+1;
0195      disp('Hit return to continue...')
0196      pause
0197      
0198      if fig_save==1
0199      fig=fig+1;
0200      end%
0201    end
0202 end
0203 aux_cal_stats=sortrows(cal_stats,2);
0204 aux_cal_stats=flipud(aux_cal_stats);
0205 fid = fopen('calibration_statistics.dat','w');
0206 fprintf(fid,' fit #  pct. data used   std. dev   residual\n ');
0207 fprintf(fid,'%2d %6.3f %9.6f %6.3e\n',aux_cal_stats');
0208 fclose(fid);
0209 aux_cal_stats1=sortrows(cal_stats,3);
0210 fid = fopen('calibration_statistics_std_dev.dat','w');
0211 fprintf(fid,' fit #  pct. data used   std. dev   residual\n ');
0212 fprintf(fid,'%2d %6.3f %9.6f %6.3e\n',aux_cal_stats1');
0213 fclose(fid);
0214 for i=3:length(aux_cal_stats)+2
0215 out1{i}=sprintf('%2d %11s  %2.4f %1.3e %1.4e ',aux_cal_stats1(i-2,1),type_fit{aux_cal_stats(i-2,1)},aux_cal_stats(i-2,2),aux_cal_stats(i-2,3),aux_cal_stats(i-2,4));
0216 end
0217 if fig_save==1
0218 out1{2}=[ '  figure  ' 'fit  ' '  %data' '   std_dev' '   mean_residual'];
0219 else
0220  out1{2}=['fit  ' '  fit  ' '  %data' '   std_dev' '   mean_residual'];
0221 end
0222 out1{1}=['fits sorted by % of data used'];
0223 s=textwrap(out1,100);
0224 msgbox(s,'fits sorted by % of data used')
0225 
0226 for i=3:length(aux_cal_stats1)+2
0227 out2{i}=sprintf('%2d %11s %2.4f %1.3e %1.4e ',aux_cal_stats1(i-2,1),type_fit{aux_cal_stats1(i-2,1)},aux_cal_stats1(i-2,2),aux_cal_stats1(i-2,3),aux_cal_stats1(i-2,4));
0228 end
0229 if fig_save==1
0230 out2{2}=['  figure  ' 'fit  '  '  %data' '   std_dev' '   mean_residual'];
0231 else
0232  out2{2}=['fit  ' '  fit  ' '  %data' '   std_dev' '   mean_residual'];
0233 end
0234 out2{1}=['fits sorted by standard deviation'];
0235 s=textwrap(out2,100);
0236 msgbox(s,'fits sorted by standard deviation')
0237 %
0238 figure
0239 plot(aux_cal_stats(:,3),aux_cal_stats(:,2),'k.')
0240 hold on
0241 ylabel(['% Data Used'])
0242 xlabel(['Standard Deviation'])
0243 text(aux_cal_stats(:,3),aux_cal_stats(:,2),num2str(aux_cal_stats(:,1)))
0244 hold off
0245 %
0246 [max_stats]=eval_fit(cond_cal.stnnbr(ii),ctdcon(ii),cond_cal.ctdprs(ii),ctddcdp(ii),btlcon(ii),cal_coeff(:,2:11),max_idx);
0247 aux_max_stats=sortrows(max_stats,2);
0248 fid = fopen('calibration_statistics_maximum_data.dat','w');
0249 fprintf(fid,' fit #  std. dev   residual\n ');
0250 fprintf(fid,'%2d %9.6f %6.3e\n',aux_max_stats');
0251 fclose(fid);
0252 for i=3:length(aux_max_stats)+2
0253 out3{i}=sprintf('%2d  %11s %1.3e %1.4e ',aux_max_stats(i-2,1),type_fit{aux_max_stats(i-2,1)}, aux_max_stats(i-2,2),aux_max_stats(i-2,3));
0254 end
0255 if fig_save==1
0256 out3{2}=['  figure  ' 'fit  '   '   std_dev' '   mean_residual'];
0257 else
0258  out3{2}=['fit  ' '  fit   ' '     std_dev' '   mean_residual'];
0259 end
0260 out3{1}=['fits sorted by standard deviation (maximum data)'];
0261 s=textwrap(out3,100);
0262 msgbox(s,'fits sorted by standard deviation (maximum data)')
0263 
0264 [min_stats]=eval_fit(cond_cal.stnnbr(ii),ctdcon(ii),cond_cal.ctdprs(ii),ctddcdp(ii),btlcon(ii),cal_coeff(:,2:11),min_idx);
0265 aux_min_stats=sortrows(min_stats,2);
0266 fid = fopen('calibration_statistics_minimum_data.dat','w');
0267 fprintf(fid,' fit #   std. dev   residual\n ');
0268 fprintf(fid,'%2d  %9.6f %6.3e\n',aux_min_stats');
0269 fclose(fid);
0270 for i=3:length(aux_min_stats)+2
0271 out4{i}=sprintf(' % 2d %11s  %1.3e  %1.4e  ',aux_min_stats(i-2,1),type_fit{aux_min_stats(i-2,1)},aux_min_stats(i-2,2),aux_min_stats(i-2,3));
0272 end
0273 if fig_save==1
0274 out4{2}=['  figure  ' 'fit  '   '   std_dev' '   mean_residual'];
0275 else
0276  out4{2}=['fit  ' '  fit   '   '   std_dev' '   mean_residual'];
0277 end
0278 out4{1}=['fits sorted by standard deviation (minimum data)'];
0279 s=textwrap(out4,100);
0280 msgbox(s,'fits sorted by standard deviation (minimum data)')
0281 
0282 aux_sum_stats=aux_max_stats(:,1:2);
0283 for i=1:length(aux_sum_stats)
0284 fit_idx=find(aux_min_stats(:,1)==aux_sum_stats(i));
0285 aux_sum_stats(i,3)=aux_min_stats(fit_idx,2);
0286 end
0287 figure
0288 plot(aux_sum_stats(:,2),aux_sum_stats(:,3),'k.')
0289 hold on
0290 ylabel(['Standard Deviation (minimum data)'])
0291 xlabel(['Standard Deviation (maximum data)'])
0292 text(aux_sum_stats(:,2),aux_sum_stats(:,3),num2str(aux_sum_stats(:,1)))
0293 hold off
0294 return
calibration_cond

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
