xbttsg_biaspertrip.m

Contenu du fichier

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT GENERAL BIAS PER TRIP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for iship =1: length(thxbt)
    disp(['Treating ' thxbt(iship).shipid])
    thg=thxbt(iship).thg;
    xbt=thxbt(iship).xbt;
    clear xbtmean
    for ix=1:length(xbt)
      xbtmean(ix)=mean(xbt(ix).xsst);
    end
     %find indices of first element of each trip;
    thtrip=cat(1,thg.thcode);
    if exist('badtrips')
      %REMOVE FLAGGED TRIPS
      if ~strcmp(badtrips(iship).ship,thxbt(iship).shipid)
        error('ship ID error in badtrips:must have same indices')
      end
      gibadtrips=[];
      for ibt=1:length(badtrips(iship).trips)
        gitrip0=find(thtrip==badtrips(iship).trips(ibt));
        gibadtrips=[gibadtrips;gitrip0];
        disp(sprintf('REMOVING TRIP %10d', ...
                     badtrips(iship).trips(ibt)))
        if isempty(gibadtrips)
          error(sprintf('TRIP %10d NOT FOUND',...
                        badtrips(iship).trips(ibt)))
        end
        if exist('pltbadtrips') & pltbadtrips
          figure(15);
          plot(cat(1,thg(gitrip0).thsst)-xbtmean(gitrip0)')
          nx=length(gitrip0);
          xlabel('XBT Shot');tirids=cat(1,xbt(gitrip0).xtirid);
          ginx=1:5:nx; %where to put tirid on x axis
          set(gca,'xlim',[1 nx+1],'xtick',ginx,'xticklabel',[]);
          yl=get(gca,'ylim');
          tx=text(ginx,yl(1)*ones(length(ginx),1),...
                  num2str(tirids(ginx)),'Rotation',40);
          ylabel('THG - XBT mean'); 
          title(sprintf('SHIP %s TRIP %10d',thxbt(iship).shipid,...
                        badtrips(iship).trips(ibt)))   
          grid on;
          %print -Phpecopsingle
        end
      end %ibt
    end
    thg(gibadtrips)=[];
    xbt(gibadtrips)=[];
    thxbt(iship).thg(gibadtrips)=[];
    thxbt(iship).xbt(gibadtrips)=[];
    thtrip=cat(1,thg.thcode);
    githtrip=[1;find(diff(thtrip))+1];
    ginpoints=[diff([githtrip;length(xbt)+1])];
    ntrip=length(githtrip);
    nxbt=length(xbt);
    tlat=cat(1,thg.thlat);
    tlon=cat(1,thg.thlon);
    clear xbtmean
    for ix=1:nxbt
      xbtmean(ix)=mean(xbt(ix).xsst);
    end
    for itrip=1:ntrip
      gitrip=githtrip(itrip):githtrip(itrip)+ginpoints(itrip)-1;
      if any(diff(cat(1,thg(gitrip).thcode)))
        error('pb with trip id')
      end
      %Mean Bias for trip 1 point minimum
      if length(gitrip)
        xdev=cat(1,thg(gitrip).thsst)-xbtmean(gitrip)';
      else
        xdev=NaN;
      end
      thxbt(iship).trip(itrip).dev=mmean(xdev);
      thxbt(iship).trip(itrip).std=mstd(xdev);
      thxbt(iship).trip(itrip).tripid=thtrip(githtrip(itrip));
      thxbt(iship).trip(itrip).gitrip=gitrip;
    end %itrip
    xdev=cat(1,thg.thsst)-xbtmean';
    thxbt(iship).dev=mmean(xdev);
    thxbt(iship).std=mstd(xdev);
    
    figure(iship);clf
    devt=cat(1,thxbt(iship).trip(:).dev);
    plot(devt,'ro');
    title(sprintf('%s mean bias, %5.2f +/- %5.2f deg',...
                  thxbt(iship).shipid,thxbt(iship).dev,...
                  thxbt(iship).std))
    stdt=cat(1,thxbt(iship).trip(:).std);
    hold on;grid on;
    ntrip=length(stdt);
    plot(1:ntrip,stdt,'b-+',1:ntrip,-stdt,'b-+');
    set(gca,'xlim',[1 ntrip+1],'xtick',1:ntrip,'xticklabel',[]);
    set(gca,'ylim',[-5 5]);
    yl=get(gca,'ylim');
    tx=text(1:ntrip,yl(1)*ones(ntrip,1),...
         num2str(cat(1,thxbt(iship).trip(:).tripid)),...
         'Rotation',40);
    ylabel('degres');xlabel('trip')
    legend('mean(tsg-xbt)','sdt dev (tsg-xbt)')
    drawnow
    %eval(sprintf('print -depsc figbiaspertrip_%s.eps',...
    %             thxbt(iship).shipid))
    eval(sprintf('print -djpeg XBT-TSG-results/figbiaspertrip-c_%s.jpg',...
                 thxbt(iship).shipid))
  end %on iship