% ti_2006_read_ctd_data.m
% Tioga-171, Tioga-177, ... CTD data: Plot CTD data profiles
% Olga Kosnyreva, AOPE/WHOI, 5/5/2006

clear;
flops(0);
%*****************************************************************
%*************  C H A N G E  C O N S T A N T S   *****************
%*****************************************************************
% ---- Path - ctd data files.
% You MUST USE YOUR PATH -----------------------------------------
% Replace CTD data PATHs with your CTD data PATHs: 
path_name_ti171='~/tioga_2006/ti171/CTD/process/';
path_name_ti177='~/tioga_2006/ti177/ctd/process/';
% ---- Axes limits:
% Your OPTION: CHOOSE AXES LIMITS --------------------------------
% (1) One station Quantity and Depth limits
%showAxesLimits='station';
% or
% (2) All stations Quantity and Depth limits
showAxesLimits='all_stations';
% ---- Down or Up cast:
% Your OPTION: CHOOSE UP or DOWN CAST ---------------------------
% (1) DOWN cast:
f_ext='cdn';
% or
% (2) UP cast:
%f_ext='cup';
%******************************************************************
%****************************       MENU    ***********************
%******************************************************************
% ---- MENU: Tioga cruise
crNum=f_menu_of_numbers_title('Tioga cruise',[171 177]);
if isempty(crNum), return; end
if crNum==171
   path_name=path_name_ti171;
   ctdList=1:29;
   date_str='April (11 - 12), 2006'; 
   cruise_title='R/V Tioga TI171';
   cruise_id='ti171';
elseif crNum==177
   path_name=path_name_ti177;
   ctdList=1:29;
   date_str='April (25 - 26), 2006';
   cruise_title='R/V Tioga TI177';
   cruise_id='ti177';
end
%******************************************************************
%***********************     C O N S T A N T S   ******************
%******************************************************************
% ---- CTD datafile:
% Number of Header lines.
headerLastLineNum=95;
% Number of variables.
n_par = 14;
% Columns#.
colPres=2;
colTemp=3;
colCnd=4;
colFluor=5; % (from v1) Fluor
colPar=6;   % (from v2) Par
colTrans=7; % (from v3) Trans
colAltim=8; % (from v4) Altim (m)
colTurb=9;  % (from v5) Turb (FTU)
colOxy=10;  % (from v6) Oxy (ml/L)
colSpar=11; % (from v9) Spar
colSal=12;
colSig=13;
% Flag 'bad data'.
bad_flag=-9.990e-29;
% ---- CTD list file columns#:
col_st=1;    
col_year=2;
col_mon=3;
col_day=4;
col_h=5;
col_min=6;
col_sec=7;
col_lat=8;
col_lon=9;
% Months names list.
monNames=strvcat('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec');
% Font size.
fnt_axes=9;
fnt_label=10;
fnt_title=10;
%**************************************************************************
%***************************     MENU QUANTITY   **************************
%**************************************************************************
parNum = menu('Quantity',...   
                          'Temperature (deg C)',... 
                          'Conductivity (S/m)',... 
                          'Fluorescence (mkg/L)',...
                          'Irradiance PAR ()',...
                          'Transmissometer ()',...
                          'Altimeter (m)',...
                          'Seapoint Turbidity Meter (FTU)',...
                          'Oxygen (ml/L)',...
                          'Surface Irradiance SPAR ()',...
                          'Salinity (psu)',...
                          'SigmaT (kg/m^3)',...
                          'CALCULATE Oxygen (mkMol/kg)',...
                          'CALCULATE Oxygen saturation (%)',...
                          'CANCEL');
if parNum>13, return; end
% C O N S T A N T S 
% depending on Quantity name. 
if parNum==1
   colNum = colTemp; 
   j_label='Temperature (deg C)';
   jN='Temperature';
   x_lim_q=[4 8];
elseif parNum==2
   colNum = colCnd; 
   j_label='Conductivity (S/m)';
   jN='Conductivity';
   x_lim_q=[3 3.6];
elseif parNum==3
   colNum = colFluor; 
   j_label='Fluorescence (\mug/L)'; 
   jN='Fluorescence';
   x_lim_q=[0 40];
elseif parNum==4
   colNum = colPar; 
   j_label='Irradiance PAR ()';
   jN='Irradiance PAR';
   x_lim_q=[0 2050];
elseif parNum==5
   colNum = colTrans; 
   j_label='Transmissometer ()';
   jN='Transmissometer';
   x_lim_q=[50 110];
elseif parNum==6
   colNum = colAltim; 
   j_label='Altimeter (m)';
   jN='Altimeter';
   x_lim_q=[0 110];
elseif parNum==7
   colNum = colTurb; 
   j_label='Seapoint Turbidity Meter ()';
   jN='Seapoint Turbidity Meter';
   x_lim_q=[0 3];
elseif parNum==8
   colNum = colOxy; 
   j_label='Oxygen (ml/L)';
   jN='Oxygen (ml/L)';
   %x_lim_q=[-1 10.5];
   x_lim_q=[4 11];
elseif parNum==9
   colNum = colSpar; 
   j_label='Surface Irradiance SPAR ()';
   jN='Surface Irradiance SPAR';
   x_lim_q=[200 2000];
elseif parNum==10
   colNum = colSal;
   j_label='Salinity (psu)';
   jN='Salinity';
   x_lim_q=[30 34.5];
elseif parNum==11
   colNum = colSig; 
   j_label='SigmaT (kg/m^3)';
   jN='SigmaT';
   x_lim_q=[23.5 27];
elseif parNum==12 % CALCULATE Oxygen in mkMol/kg
   colNum = colOxy; 
   j_label='Oxygen (\muMol/kg)';
   jN='Oxygen (\muMol/kg)';
   %x_lim_q=[-40 460];
   x_lim_q=[150 460];
elseif parNum==13 % CALCULATE Oxygen Saturation in %
   colNum = colOxy; 
   j_label='Oxygen saturation (%)';
   jN='Oxygen saturation (%)';
   %x_lim_q=[-20 150];
   x_lim_q=[50 150];
end
%**************************************************************************
%***************************     MENU STATION#   **************************
%**************************************************************************
stNum=f_menu_of_numbers_title('Station #',ctdList);
if isempty(stNum), return; end
%**************************************************************************
%***************************  LOAD CTD LIST DATAFILE   ********************
%**************************************************************************
p=['load ti' num2str(crNum) '_ctd_list;'];
eval(p);
p=['c_list=ti' num2str(crNum) '_ctd_list;'];
eval(p);
%**************************************************************************
%************************   STATION DATE/TIME/CAST  ***********************
%**************************************************************************
% Date.
r_st=find(c_list(:,col_st)==stNum);
yyyy_st=c_list(r_st,col_year);
mm_st=c_list(r_st,col_mon);
dd_st=c_list(r_st,col_day);
% Time.
hh_st=c_list(r_st,col_h);
minmin_st=c_list(r_st,col_min);
% Date and time - string represetation.
mm_st = deblank(monNames(mm_st,:));
str_date=[num2str(dd_st) '-' mm_st '-' num2str(yyyy_st)];
str_time=[num2str(hh_st) ':' num2str(minmin_st) ' UTC'];
% DOWN or UP cast ID.
if strcmp(f_ext,'cdn')
  cast_id='DOWN';
elseif strcmp(f_ext,'cup')
  cast_id='UP';
else
   return;
end
% String - plot title.
str_title(1)={[cruise_title ', Station ' num2str(stNum) ' ' cast_id ' cast, ' jN]};
str_title(2)={[str_date ' ' str_time]};
%**************************************************************************
%***************************        FIGURE      ***************************
%**************************************************************************
figure;
set(gcf,'color',[1 1 1],'position',[413   553   303   376],...
     'paperPositionMode','auto');
%**************************************************************************
%***************************   READ  CTD DATAFILE    **********************
%**************************************************************************
str_disp='"BAD" DATA INFO:';
         % Write input file names to string.
         if stNum < 10
            j = ['ti' num2str(crNum) '00' num2str(stNum) '.' f_ext];
         elseif (stNum<100) & (stNum>9) 
            j = ['ti' num2str(crNum) '0' num2str(stNum) '.' f_ext];
         else 
            j = ['ti' num2str(crNum) num2str(stNum) '.' f_ext];
         end
         % Open file for reading.
         pfid = fopen([path_name j]);
         % Read first lines (headers).
         for lineNum=1:1:headerLastLineNum
            headerLine=fgets(pfid);
         end
         % Read formated data from file.
         da=fscanf(pfid,'%f',[1 inf])';
         % Reshape data matrix.
         data = (reshape(da,n_par,size(da,1)/n_par))';
         % Close file.
         fclose(pfid);
         % Current station Pressure.
         presSt=data(:,colPres);
         % Current station quantity.
         if parNum==12
            % ----------Oxygen Concentration in mkMol/kg.
            OX = data(:,colOxy);
            T=data(:,colTemp);
            S=data(:,colSal);
            P=zeros(size(T));
            % (2) Find and remove "bad" data.
            r_bad=find(P(:,1)==bad_flag | T(:,1)==bad_flag | S(:,1)==bad_flag | OX(:,1)==bad_flag);
            if ~isempty(r_bad)
               str_disp=strvcat(str_disp,['Station No ' num2str(stNum) ', ' jN ': ' num2str(length(r_bad)) ' line(s) with "bad" data (flag=-9.990e-29)'],...
                 'Pressure  Quantity',num2str([presSt(r_bad,:) parSt(r_bad,:)]));
               T(r_bad,:)=nan;
               S(r_bad,:)=nan;
               P(r_bad,:)=nan;
               OX(r_bad,:)=nan;
               presSt(r_bad,:)=nan;
            end
            [swan,sigT]=swstate(S,T,P);
            parSt = OX * 44.66 ./ (1000 + sigT) * 1000; % Oxygen (mkMol/kg)
       elseif parNum==13
            % ----------Oxygen Saturation in %.
            OX = data(:,colOxy);
            T=data(:,colTemp);
            S=data(:,colSal);
            % Find and remove "bad" data.
            r_bad=find(T(:,1)==bad_flag | S(:,1)==bad_flag | OX(:,1)==bad_flag);
            if ~isempty(r_bad)
               str_disp=strvcat(str_disp,['Station No ' num2str(stNum) ', ' jN ': ' num2str(length(r_bad)) ' line(s) with "bad" data (flag=-9.990e-29)'],...
                 'Pressure  Quantity',num2str([presSt(r_bad,:) parSt(r_bad,:)]));
               T(r_bad,:)=nan;
               S(r_bad,:)=nan;
               OX(r_bad,:)=nan;
               presSt(r_bad,:)=nan;
            end
            % Oxygen solubility in ml/l (Ocean Toolbox, Phil Morgan): 
            soly=sw_satO2(S,T);
            % Oxygen Saturation in %.
            parSt = OX ./ soly *100;
         else
            parSt=data(:,colNum);
            % Find and remove "bad" data.
            r_bad=find(parSt(:,1)==bad_flag);
            if ~isempty(r_bad)
               str_disp=strvcat(str_disp,['Station No ' num2str(stNum) ', ' jN ': ' num2str(length(r_bad)) ' line(s) with "bad" data (flag=-9.990e-29)'],...
                 'Pressure  Quantity',num2str([presSt(r_bad,:) parSt(r_bad,:)]));
               parSt(r_bad,:)=nan;
               presSt(r_bad,:)=nan;
            end
         end
%**************************************************************************
%***************************    Plot profiles   ***************************
%**************************************************************************
        plot(parSt,presSt,'-','color',[1 0 0],'linewidth',1.5);
%**************************************************************************
%*******************************    AXES LIMITS   *************************
%**************************************************************************
if strcmp(showAxesLimits,'all_stations')
   x_lim=x_lim_q;
   y_lim=[0 275];
   str_t(1)={['Min=' num2str(min(parSt))]};
   str_t(2)={['Max=' num2str(max(parSt))]};
   text(x_lim(1)+range(x_lim)*0.5,y_lim(1)+range(y_lim)*0.9,str_t);
else
   % Quantity and pressure range - for AXES limits.
   min_max_info=[min(parSt) max(parSt) min(presSt) max(presSt)];
   min_max_info=[min(parSt) max(parSt) min(presSt) max(presSt)];
   par_range=[min(min_max_info(:,1)) max(min_max_info(:,2))];
   pres_range=[min(min_max_info(:,3)) max(min_max_info(:,4))];
   % X-axes limits.
   if min(parSt)==0
      x_lim_1=0;
   else
      x_lim_1=par_range(1)-0.05*range(par_range);
   end
   x_lim_2 = par_range(2)+0.05*range(par_range);
   x_lim=[x_lim_1 x_lim_2];
   % Y-axes limits.
   y_lim_1=0;
   y_lim_2 = pres_range(2)+0.05*range(pres_range);
   y_lim=[y_lim_1 y_lim_2];
end
% Set Axes limits.
axis([x_lim y_lim]);
%**************************************************************************
%****************************    AXES PROPERTIES   ************************
%**************************************************************************
grid on;
set(gca,'ydir','reverse','fontsize',fnt_axes);
xlabel(j_label,'fontsize',fnt_label);
ylabel('Pressure (db)','fontsize',fnt_label);
title(str_title,'fontsize',fnt_title,'color',[0 0 1]);
% Change axes position.
ah=get(gca,'position');
set(gca,'position',[ah(1)+0.02*ah(3) ah(2)-0*ah(4) ah(3) ah(4)*0.95]);
if size(str_disp,1)>1
   disp(str_disp);    
end
