function plot_cycle_slips( filePrefix, ObsRec )
% 
% plot_cycle_slips()
%
% Input:
%  filePrefix - String that will be at the beginning of each plot filename
%  ObsRec     - viewdat output of record 27 or 35:19
%
% viewdat -d27 -h -translate_rec35_sub19_to_rec27 --mat=Obs.mat Input.T04 
% or
% viewdat -d35:19 -h --mat=Obs.mat Input.T04 
%
%  
% In Matlab:-
% Obs = load_viewdat('Obs.mat');
% plot_cycle_slips('MyString',Obs);
%
% Useful utility to plot tracked data along with cycle slips for GPS,
% GLONASS ... observables.
%
%
% Developed by Andrew Cartmell
% Copyright Trimble Inc
%

% Rec35:19 has time in the first columnd, rec 27 has seconds. We are
% unlikely to have a data file that spans beyond one week rollover so
% check the range of the fist column and make an assumption based on
% this range
minCol1 = min(ObsRec(:,1));
maxCol1 = max(ObsRec(:,1));

Col1Deta = maxCol1 - minCol1;

if(Col1Deta > 1)
  RecType27 = 1;
  disp('Assume Rec 27')

  key.TIME       = 1;  % GPS Time [Sec]
  key.CLKBIAS    = 2;  % Clock Offset [ms]
  key.GPSGLNBIAS = 3;  % GPS-GLONASS clock offset [ms]
  key.RAIM       = 4;  % RAIM INFO
  key.EPOCH      = 5;  % Epoch Flags (Deleted, System offset present, RAIM info present)
  key.SV         = 6;  % SV ID
  key.FDMA       = 7;  % FDMA Channel (GLONASS) otherwise the hardware channel number
  key.SAT_TYPE   = 8;  % Satellite Type
  key.SV_FLAGS   = 9;  % Satellite Flags 
  key.AZ         = 10; % Azimuth [deg]
  key.EL         = 11; % Elevation [deg]
  key.FREQ       = 12; % Frequency Band
  key.TRACK      = 13; % Track Type
  key.PHASE      = 14; % Phase [Cycles]
  key.RANGE      = 15; % Pseudorange [m]
  key.DOPP       = 16; % Doppler [Hz]
  key.CNO        = 17; % C/No [dBHz]
  key.SLIP_FLG   = 18; % Cycle Slip Flag
  key.SLIP_CNT   = 19; % Cycle Slip Counter
  key.MEAS1      = 20; % Measurement Flags 
  key.ANTENNA    = 21; % Antenna #
  key.MEAS2      = 22; % Measurement Flags 2
  key.MEAS3      = 23; % Measurement Flags 3
else
  RecType27 = 0;
  disp('Assume Rec 35:19')

  % Record 35:19 fields
  key.WN       = 1;  % GPS Week Number
  key.TIME     = 2;  % GPS Time [Sec]
  key.CLKBIAS  = 3;  % Clock Offset [ms]
  key.RAIM     = 4;  % RAIM INFO
  key.SV       = 5;  % SV ID
  key.FDMA     = 6;  % FDMA Channel (GLONASS)
  key.channel  = 7;  % Receiver Channel
  key.SAT_TYPE = 8;  % Satellite Type
  key.SV_FLAGS = 9;  % Satellite Flags
  key.AZ       = 10; % Azimuth [deg]
  key.EL       = 11; % Elevation [deg]
  key.FREQ     = 12; % blockType = Frequency Band
  key.TRACK    = 13; % trackType = Track Type
  key.PHASE    = 14; % Phase [Cycles]
  key.RANGE    = 15; % Pseudorange [m]
  key.DOPP     = 16; % Doppler [Hz]
  key.CNO      = 17; % C/No [dBHz]
  key.SLIP_CNT = 18; % Cycle Slip Counter
  key.MEAS     = 19; % Measurement Flags
  key.ANTENNA  = 20; % Antenna #
  key.CSTATE   = 21; % CSTATE in rec35:19 - not st_iface*.h value

  flag.dMEAS_PHASE_VALID    = 1;
  flag.dMEAS_DOPP_VALID     = 2;
  flag.dMEAS_RANGE_VALID    = 4;
  flag.dMEAS_LOCK_POINT     = 8;
  flag.dMEAS_POS_LOCK_POINT = 16;
  flag.dMEAS_XCORR_WAIT     = 32;
  flag.dMEAS_UNHEALTHY      = 64;
  flag.dMEAS_MULTIPATH      = 128;
  flag.dMEAS_SLIP           = 256;
  flag.dMEAS_MASTER_SUB     = 512;
  flag.dMEAS_REC27          = 1024;
  flag.dMEAS_JAMMER         = 2048;
  flag.dMEAS_CODE_PULLIN    = 4096;
  flag.dMEAS_SCINT_LO       = 8192;
  flag.dMEAS_SCINT_HI       = 16384;
end


% Week extend the data - useful if we want to look at data that crosses
% the week boundary
if(RecType27 == 0)
  % For rec 35:19 we have the week number
  if(ObsRec(end,key.WN) ~= ObsRec(1,key.WN))
    ObsRec(:,key.TIME) = ObsRec(:,key.TIME) + ((ObsRec(:,key.WN) - ObsRec(1,key.WN)) * 7 * 86400);
  end
else
  if(ObsRec(end,key.TIME) < ObsRec(1,key.TIME))
    % Time went backwards
    delta = diff(ObsRec(:,key.TIME));
    % Find where it jumps backwards
    i = find(delta < -3.5*86400);
    % Will only handle 1 week roll
    ObsRec((i(1)+1):end,1) = ObsRec((i(1)+1):end,1) + 86400 * 7;
  end
end

%
% A few options (could make parameters on the input)
%

% Show all bands, if 0 just L1 C/A
show_all_bands = 1;
% Plot the X-axis in hours (0) or seconds (1)
plot_sec_axis  = 1;
% Only valid when the x-axis is in hours, force the X-axis to 24 hours
force_to_one_day = 0;

% Take a guess at the data rate
drate = 1.0;
stime = ObsRec(floor(length(ObsRec)/2),key.TIME);
for tmpx=floor(length(ObsRec)/2):floor(length(ObsRec)/2)+100
  if ( stime ~= ObsRec(tmpx,key.TIME) )
    drate = ObsRec(tmpx,key.TIME) - stime;
    break;
  end
end
fprintf('Guessed data is rate %.2f\n',drate);

% Align in days and convert to seconds
day= floor(ObsRec(1,key.TIME)/86400);


if(plot_sec_axis)
  Scale = 1;
else
  Scale = 3600;
end
drate = drate/Scale;
ObsRec(:,key.TIME) = (ObsRec(:,key.TIME) - day*86400) / Scale;
taxis = round( (ObsRec(1,key.TIME):drate:ObsRec(end,key.TIME)) * 1000 * Scale);

%for sy = 0 : 10 % RT_SatType_e
%fprintf('Set to GPS/GLN only\n');
for sy = [0, 1, 2, 3, 4, 10, 9] % RT_SatType_e
  figure
  clear figLabel;
  clear pltHandle;
  figIndex = 1;
  % Following gotXXX used for the legend generation
  gotL1CA = 0;
  gotL1C = 0;
  gotL1E = 0;
  gotL2C = 0;
  gotL2E = 0;
  gotL5  = 0;

  clf; hold on; grid on;
  yind = 1;
  ytlbl = char(1,4);

  if ( sy == 1 )
    sv_range=120:158;
    SVString = 'SBAS';
  elseif(sy  == 0)
    sv_range=1:32;
    SVString = 'GPS';
  elseif(sy  == 2)
    sv_range=1:24;
    SVString = 'GLN';
  elseif(sy  == 3)
    sv_range=1:36;
    SVString = 'GAL';
  elseif(sy  == 4)
    sv_range=193:197;
    SVString = 'QZSS';
  elseif(sy  == 10)
    sv_range=1:30;
    SVString = 'BDS';
  elseif(sy  == 9)
    sv_range=1:7;
    SVString = 'IRNSS';
  else
    sv_range=1:32;
    SVString = 'Unknown';
  end

  for sv = sv_range
    ind = find((ObsRec(:,key.SV)==sv)&(ObsRec(:,key.SAT_TYPE)==sy));
    if ( ~isempty(ind) )
      iL1CA=[]; iL1C=[]; iL1E=[]; iL2E=[]; iL2C=[]; iL5=[];

      if ( sy == 0 )
        % [iL1CA, iL1C, iL1E, iL2E, iL2C, iL5] = find_gps_ind( arr );
        [iL1CA, iL1C, iL1E, iL2E, iL2C, iL5] = find_gps_ind(ObsRec(ind,:));
        sychar = 'G';
      elseif ( sy == 1 )
        % [iL1CA, iL5] = find_sbas_ind( arr );
        [iL1CA, iL5] = find_sbas_ind(ObsRec(ind,:));
        sychar = 'S';
      elseif ( sy == 2 )
        % [iG1C, iG1P, iG2C, iG2P, iG3] = find_gln_ind( arr );
        [iL1CA, iL1E, iL2C, iL2E, iL5] = find_gln_ind(ObsRec(ind,:));
        sychar = 'R';
      elseif ( sy == 3 )
        % [iE1, iE5A, iE5B, iE5L, iE6] = find_gal_ind( arr );
        [iL1CA, iL1E, iL2E, iL5, ~] = find_gal_ind(ObsRec(ind,:));
        sychar = 'E';
      elseif ( sy == 4 )
        % [iL1CA, iL1C, iL1SAIF, iL2C, iL5, iLEX] = find_qzss_ind( arr );
        [iL1CA, iL1C, iL1E, iL2C, iL5, ~] = find_qzss_ind(ObsRec(ind,:));
        sychar = 'Q';
      elseif ( sy == 5 || sy == 7 || sy == 10 )
        % [iB1, iB2, iB3] = find_cmp_ind( arr );
        [iL1CA, iL5, iL2E] = find_cmp_ind(ObsRec(ind,:));
        sychar = 'B';
      elseif ( sy == 9 )
        % [iL5CA, iS1CA] = find_irss_ind( arr );
        [iL1CA, iL1C] = find_irnss_ind(ObsRec(ind,:));
        sychar = 'I';
      end

      % Blue    = L1C/A, G1C/A, E1, B1
      % Black   = L1C, S1CA
      % Red     = L1E, G1P, E5A, L1SAIF
      % Green   = L2C, G2C
      % Cyan    = L2E, G2P, E5B
      % Magenta = L5, G3, E5AltBOC
      yoff = .0;
      if (~isempty(iL1CA))
        ObsRecms=round(ObsRec(ind(iL1CA),key.TIME)*1000*Scale);
        [~,A,B]=intersect(taxis,ObsRecms);
        tmpx = NaN*taxis; tmpx(A)=ObsRec(ind(iL1CA(B)),key.TIME);
        plot(tmpx,ones(length(tmpx),1)*yind + yoff,'b');
        
        % Cycle slips
        if(RecType27 == 1)
          cind = find(ObsRec(ind(iL1CA),key.SLIP_FLG));
        else
          % USe the slip flag in the measurement flag
          cind = find(bitand(ObsRec(ind(iL1CA),key.MEAS),flag.dMEAS_SLIP));

          % Alternate way to find slips - look for changes in slip count
          %delta = diff(ObsRec(ind(iL1CA),key.SLIP_CNT));
          %cind = find(delta) + 1;
        end

        fig = plot(ObsRec(ind(iL1CA(cind)),key.TIME),ones(length(cind),1)*yind + yoff,'bo');
        
        if ( gotL1CA == 0 && ~isempty(fig) )
          pltHandle(figIndex) = fig;
          if ( sy == 5 || sy == 7 || sy == 10 )
            figLabel(figIndex) = {'B1'};
          elseif ( sy == 3 )
            figLabel(figIndex) = {'E1'};
          elseif ( sy == 9 )
            figLabel(figIndex) = {'L5 C/A'};
          else
            figLabel(figIndex) = {'L1 C/A'};
          end
          figIndex = figIndex + 1;
          gotL1CA = 1;
        end

        % SV, StartElev->EndElev, # Epochs, # Slip Epochs, % Slip Epochs
        fprintf( 'L1CA %c%3d %4.1f->%4.1f %6d %6d %6.2f\n', ...
                 sychar, sv, ...
                 ObsRec(ind(iL1CA(1)),key.EL), ObsRec(ind(iL1CA(end)),key.EL), ...
                 length(iL1CA), length(cind), 100.0*length(cind)/length(iL1CA) ...
               );
        if ( show_all_bands && ~isempty(iL1C) )
          yoff = yoff + .15;
          %plot(ObsRec(ind(iL1C),1),ones(length(iL1C),1)*yind + yoff,'k.');
          ObsRecms=round(ObsRec(ind(iL1C),key.TIME)*1000*Scale);
          [~,A,B]=intersect(taxis,ObsRecms);
          tmpx = NaN*taxis; tmpx(A)=ObsRec(ind(iL1C(B)),key.TIME);
          plot(tmpx,ones(length(tmpx),1)*yind + yoff,'k');

          if(RecType27 == 1)
            cind = find(ObsRec(ind(iL1C),key.SLIP_FLG));
          else
            cind = find(bitand(ObsRec(ind(iL1C),key.MEAS),flag.dMEAS_SLIP));
          end
          
          fig = plot(ObsRec(ind(iL1C(cind)),key.TIME),ones(length(cind),1)*yind + yoff,'ko');
          if ( gotL1C == 0 && ~isempty(fig) )
            pltHandle(figIndex) = fig;
            if ( sy == 9 )
              figLabel(figIndex) = {'S1 C/A'};
            else
              figLabel(figIndex) = {'L1 C'};
            end
            figIndex = figIndex + 1;
            gotL1C = 1;
          end
        end
        if ( show_all_bands && ~isempty(iL1E) )
          yoff = yoff + .15;
          %plot(ObsRec(ind(iL1E),1),ones(length(iL1E),1)*yind + yoff,'r.');
          ObsRecms=round(ObsRec(ind(iL1E),key.TIME)*1000*Scale);
          [~,A,B]=intersect(taxis,ObsRecms);
          tmpx = NaN*taxis; tmpx(A)=ObsRec(ind(iL1E(B)),key.TIME);
          plot(tmpx,ones(length(tmpx),1)*yind + yoff,'r');
          if(RecType27 == 1)
            cind = find(ObsRec(ind(iL1E),key.SLIP_FLG));
          else
            cind = find(bitand(ObsRec(ind(iL1E),key.MEAS),flag.dMEAS_SLIP));
          end

          fig = plot(ObsRec(ind(iL1E(cind)),key.TIME),ones(length(cind),1)*yind + yoff,'ro');
          if ( gotL1E == 0 && ~isempty(fig) )
            pltHandle(figIndex) = fig;
            if ( sy == 3 )
              figLabel(figIndex) = {'E5A'};
            elseif ( sy == 4 )
              figLabel(figIndex) = {'L1SAIF'};
            else
              figLabel(figIndex) = {'L1E'};
            end
            figIndex = figIndex + 1;
            gotL1E = 1;
          end
        end
        if ( show_all_bands && ~isempty(iL2C) )
          yoff = yoff + .15;
          %plot(ObsRec(ind(iL2C),1),ones(length(iL2C),1)*yind + yoff,'g.');
          ObsRecms=round(ObsRec(ind(iL2C),key.TIME)*1000*Scale);
          [~,A,B]=intersect(taxis,ObsRecms);
          tmpx = NaN*taxis; tmpx(A)=ObsRec(ind(iL2C(B)),key.TIME);
          plot(tmpx,ones(length(tmpx),1)*yind + yoff,'g');
          if(RecType27)
            cind = find(ObsRec(ind(iL2C),key.SLIP_FLG));
          else
            cind = find(bitand(ObsRec(ind(iL2C),key.MEAS),flag.dMEAS_SLIP));
          end
          fig = plot(ObsRec(ind(iL2C(cind)),key.TIME),ones(length(cind),1)*yind + yoff,'go');
          if ( gotL2C == 0 && ~isempty(fig) )
            pltHandle(figIndex) = fig;
            if ( sy == 2 )
              figLabel(figIndex) = {'L2 C/A'};
            else
              figLabel(figIndex) = {'L2C'};
            end
            figIndex = figIndex + 1;
            gotL2C = 1;
          end
        end
        if ( show_all_bands && ~isempty(iL2E) )
          yoff = yoff + .15;
          %plot(ObsRec(ind(iL2E),1),ones(length(iL2E),1)*yind + yoff,'c.');
          ObsRecms=round(ObsRec(ind(iL2E),key.TIME)*1000*Scale);
          [~,A,B]=intersect(taxis,ObsRecms);
          tmpx = NaN*taxis; tmpx(A)=ObsRec(ind(iL2E(B)),key.TIME);
          plot(tmpx,ones(length(tmpx),1)*yind + yoff,'c');
          if(RecType27)
            cind = find(ObsRec(ind(iL2E),key.SLIP_FLG));
          else
            cind = find(bitand(ObsRec(ind(iL2E),key.MEAS),flag.dMEAS_SLIP));
          end
          fig = plot(ObsRec(ind(iL2E(cind)),key.TIME),ones(length(cind),1)*yind + yoff,'co');
          if ( gotL2E == 0 && ~isempty(fig) )
            pltHandle(figIndex) = fig;
            if ( sy == 3 )
              figLabel(figIndex) = {'E5B'};
            else
              figLabel(figIndex) = {'L2E'};
            end
            figIndex = figIndex + 1;
            gotL2E = 1;
          end
        end
        if ( show_all_bands && ~isempty(iL5) )
          yoff = yoff + .15;
          %plot(ObsRec(ind(iL5),1),ones(length(iL5),1)*yind + yoff,'m.');
          ObsRecms=round(ObsRec(ind(iL5),key.TIME)*1000*Scale);
          [~,A,B]=intersect(taxis,ObsRecms);
          tmpx = NaN*taxis; tmpx(A)=ObsRec(ind(iL5(B)),key.TIME);
          plot(tmpx,ones(length(tmpx),1)*yind + yoff,'m');
          if(RecType27)
            cind = find(ObsRec(ind(iL5),key.SLIP_FLG));
          else
            cind = find(bitand(ObsRec(ind(iL5),key.MEAS),flag.dMEAS_SLIP));
          end
          fig = plot(ObsRec(ind(iL5(cind)),key.TIME),ones(length(cind),1)*yind + yoff,'mo');
          if ( gotL5 == 0 && ~isempty(fig) )
            pltHandle(figIndex) = fig;
            if ( sy == 5 || sy == 7 || sy == 10 )
              figLabel(figIndex) = {'B2'};
            elseif ( sy == 3 )
              figLabel(figIndex) = {'E5Alt'};
            else
              figLabel(figIndex) = {'L5'};
            end
            figIndex = figIndex + 1;
            gotL5 = 1;
          end
        end

        ytlbl(yind,1:4) = sprintf('%c%3d',sychar,sv);
        yind = yind + 1;
      end
    end
  end
  set(gca(),'YTick',1:yind-1);
  set(gca(),'YTickLabel',ytlbl);

  if(plot_sec_axis)
    xlabel('Time [GPS Secs]');
  else
    if(force_to_one_day)
      a = axis;
      a(1) = 0;
      a(2) = 24;
      a(3) = 0;
      a(4) = yind;
      axis(a);
    end
    xlabel('Time [Hour into GPS day]');
  end
 
  if(exist('pltHandle'))
    legend(pltHandle,figLabel,'location','best');
  end
  
  title([SVString,' Tracking']);

  if(exist('OCTAVE_VERSION', 'builtin') ~= 0)
    print('-dpng',[filePrefix,'-',SVString,'-Track.png'],'-S960,720');
  else
    print('-dpng',[filePrefix,'-',SVString,'-Track.png']);
  end
end



