% Display theta-S and theta-O2 for CTD oxygen and bottles
% A. Ganachaud, 2014 + Glazeo 2015
%---------------------------------------------------------
% load oxygen coefficients (and btl0 variable)
strf=sprintf('load Califiles/ox_cals_%s_is%d_ws%d_%s.mat',cname,iswitch,wswitch,csuf);
disp(strf)
eval(strf) %loads calib matrix 'final'


%set subplot axes
axTS(1)=subplot(3,2,1); %for depth range 1
axTS(2)=subplot(3,2,3); %for depth range 2
axTS(3)=subplot(3,2,5); %for depth range 3
axTO(1)=subplot(3,2,2); %for depth range 1
axTO(2)=subplot(3,2,4); %for depth range 2
axTO(3)=subplot(3,2,6); %for depth range 3

% modif glazeo - 21/01/2015
read_moorspice_nc % this is to know what are the profiles (casts)

if p_onlyonecast %takes only the casts for which o2 samples were made
    if ~exist('giprof2cal')
        iprofmin=find(floor(gprof2get/100)==min(gistat2cal));
        iprofmax=find(floor(gprof2get/100)==max(gistat2cal));
        gprof2cal=gprof2get(iprofmin:iprofmax);
    end
    giprof=[];
    for iproff=1:length(gprof2cal)
        giprof=[giprof; find(actd.profileid==gprof2cal(iproff))];
    end
    giprof=giprof(:);
    profbeg=min(giprof);
    profend=max(giprof);
else %otherwise takes all casts
    gipbeg=find(actd.profileid==100*statbeg+castbeg);
    gipend=find(actd.profileid==100*statend+castend);
    giprof=gipbeg:gipend;
    if isempty(giprof),error('check castbeg and castend exist'),end
end

iswitch0=iswitch;
wswitch0=wswitch;
%
%ctd. variable contains current ctd profile while actd contails all profiles
%
varstr={'ctd.dpres' 'ctd.dtemp68' 'ctd.dsali' 'ctd.oxcal'};
for ipp=1:length(giprof)
    iprof=giprof(ipp);
    if iprof~=0
        %ctd. variable contains current ctd profile while actd contails all
        %profiles (see read_pandora_nc)
        disp(        '*********************************')
        disp(sprintf('**** STATION %03d, CAST %02d ****',Stat(iprof),Cast(iprof)))
        disp(        '*********************************')    
        %
        % loads the current ctd profiles
        %   (ctd iprof iswitch wswitch csuf)
        eval(sprintf(...
            'load Caliprofiles/ox_cal_%s_STN%03d_CAST%02d.mat ',...
            cname,actd.stat(iprof),actd.cast(iprof)))
        if ctd.iswitch~=iswitch0 || ctd.wswitch~=wswitch0
            error('wrong data set'),
        end
        for iii=1:size(gstatgroup,1)
            if ctd.stat>=gstatgroup(iii,1) && ctd.stat <= gstatgroup(iii,2)
                icolor=iii;
            end
        end
        %
        %plot bottles
        %
        % selects bottle data for graphics, later on
        %select bottles for that profile (for graphic purpose)
        gistatb=find((ismember(btl0.stat,actd.stat(iprof))==1)&(ismember(btl0.cast,actd.cast(iprof))==1));
        if isempty(gistatb)
            gistatb=find((ismember(btl0.stat,Stat(iprof))==1));
            %get bottles from other cast of that station
            p_othercast=1;
        else
            p_othercast=0;
        end
        %variable btl contains only bottles for current ctd profile
        btl=keep_data(btl0, gistatb);
        %find flagged bottles for graphics
        %good bottles
        gigoodbtl=find(btl.wnklo2_flag==2);
        %bad bottles
        gibadbtl=find(btl.wnklo2_flag==3);
        %measured salinit bottles:
        
        gigoodbtls=find(btl.bsali~=35); %35 exact means no data
        btl.theta=sw_ptmp(btl.usali,btl.utemp68,btl.upres,1000);
        ctd.theta=sw_ptmp(ctd.dsali,ctd.dtemp68,ctd.dpres,1000);
                            
        %%%%%%%%%%%%%%%
        % plot CTD data
        %%%%%%%%%%%%%%%%
        for idepr=1:length(gdeprange)-1
            gidep=find(ctd.dpres>gdeprange(idepr) & ctd.dpres<=gdeprange(idepr+1));% & oxi_qc==2);
            %TS
            axes(axTS(idepr))
            plot(ctd.dsali(gidep),ctd.theta(gidep),[gicolor{icolor} 'o'],'markersize',2); hold on; grid on
            ylabel(sprintf('Theta-1000, %5d-%5d m',gdeprange(idepr),gdeprange(idepr+1)))
            title('Salinity')
            %O2S
            axes(axTO(idepr))
            plot(ctd.oxcal(gidep),ctd.theta(gidep),[gicolor{icolor} 'o'],'markersize',2); hold on; grid on
            %if icolor==4
             %   keyboard
            %end
        end
        
        %%%%%%%%%%%%%%%%%%
        % now plot bottles
        %%%%%%%%%%%%%%%%%%%
        for idepr=1:length(gdeprange)-1
            % TS
            gidep=find(btl.upres(gigoodbtls)>gdeprange(idepr) & btl.upres(gigoodbtls)<=gdeprange(idepr+1));
            axes(axTS(idepr))
            ylabel(sprintf('Theta-1000, %5d-%5d m',gdeprange(idepr),gdeprange(idepr+1)))
            title('Salinity')
            plot(btl.usali(gigoodbtls(gidep)),btl.theta(gigoodbtls(gidep)),[gicolor{icolor} 'o'],'markersize',2); hold on; grid on
            plot(btl.usali(gigoodbtls(gidep)),btl.theta(gigoodbtls(gidep)),['ko'],'markersize',8)
            %O2-S
            gidep=find(btl.upres(gigoodbtl)>gdeprange(idepr) & btl.upres(gigoodbtl)<=gdeprange(idepr+1));
            axes(axTO(idepr))
            title('Oxygen (QC=2 for bottles) / Theta-1000')
            % modif glazeo - definition of btl.wnklo2, 21/01/2015
            btl.wnklo2 = btl.wnklo2um ; 
            % end of modification
            pll=plot(btl.wnklo2(gigoodbtl(gidep)),btl.theta(gigoodbtl(gidep)),[gicolor{icolor} 'x'],'markersize',8); hold on; grid on
            plot(btl.wnklo2(gigoodbtl(gidep)),btl.theta(gigoodbtl(gidep)),['ko'],'markersize',8);
            if ~isempty(pll)
                hdlp(icolor)=pll; %save to put legend later
            end
            
        end
        
    end
    
end
legend(axTO(1),hdlp,gstatgroupstr,'location','northeast')

