function   [ n_seg_vec, xx_vec ] = make_segments( photon_times )
% function [ n_seg_vec, xx_vec ] = make_segments( photon_times )
%
% Input:  photon_times -- photon arrival times, in microseconds
%
% Output: n_seg_vec  -- array of changepoint times
%         xx_vec     -- count rates (c/usec) in the corresponding segments
%
% Note: t_seg = photon_times( n_seg_vec ) generates the changepoint times
%
% ------------------------------------------------------------------------

global prior_ratio min_photons

n_times = length( photon_times );
  min_time = photon_times(       1 );
  max_time = photon_times( n_times );
    delta_t  = ( max_time - min_time ) / ( n_times - 1 );
    min_tick = floor( min_time - 0.5 * delta_t );
    max_tick =  ceil( max_time + 0.5 * delta_t );
    n_ticks = max_tick - min_tick + 1; % Number of microsecond "ticks"

nseg_1_vec = [ 1 ];
nseg_2_vec = [ n_times ];
nosubs_vec = [ 0 ];
xx_vec = [ n_times / n_ticks ];
no_seg_flag = 0;

while no_seg_flag == 0

   num_segments = length( nseg_1_vec );
   no_seg_flag = 1; % set escape unless do a sub-segmentation

   nseg_1_work = [];
   nseg_2_work = [];
   nosubs_work = [];
   xx_work  = [];

   for seg_id = 1: num_segments

      do_it = 1; % do this one, unless ...

      % ... this one has been done before!
      if nosubs_vec( seg_id ) == 1
         do_it = 0;
      end

      nseg_1 = nseg_1_vec( seg_id );
      nseg_2 = nseg_2_vec( seg_id );
      x_seg  = xx_vec( seg_id );

      times_this = photon_times( nseg_1: nseg_2 );
      nt_this = length( times_this );

      if do_it > 0

         % Determine previous time
         time_this_1 = times_this(1);

         if time_this_1 == photon_times(1);

            % Special handling for first point in full array,
            %   or if it is the second point, but the first two
            %   (or more) times are equal:
            ii = find( times_this > time_this_1 );
            if isempty(ii) 
               delt_t = 2; % Token value
            else
               delt_t = times_this(ii(1)) - time_this_1;
            end
            t_0 = time_this_1 - delt_t;

         else

            % t_0 is the time just previous to the sub-array
            t_0 = photon_times( nseg_1 - 1 );

         end

         % Determine subsequent time
         time_this_n = times_this(nt_this);

         if time_this_n == photon_times(n_times);

            % Special handling for last point in full array,
            %   or if it is the second-to-last point, but the 
            %   last two (or more) times are equal:
            ii = find( times_this < time_this_n );
            if isempty(ii) 
               delt_t = 2; % Token value
            else
               delt_t = time_this_n - times_this(ii(length(ii)));
            end
            t_n = time_this_n + delt_t;

         else

            % t_n is the time just after the sub-array
            t_n = photon_times( nseg_2 + 1 );

         end

         [ n_seg, odds_ratio, log_prob ] = find_change( times_this, t_0, t_n );

         % ... one of the proposed subsegments is too short:
         n_seg_right = nt_this - n_seg;
         if (n_seg <= min_photons) | (n_seg_right <= min_photons)
            do_it = 0; 
         end

         % ... the significance criterion not met:
         if odds_ratio < prior_ratio
            do_it = 0; 
         end

      end

      if do_it > 0

         % Subsegment this one; do not escape yet
         no_seg_flag = 0; 

         ntimes_1_left  = nseg_1;
         ntimes_1_right = nseg_1 + n_seg - 1;

         ntimes_2_left  = nseg_1 + n_seg;
         ntimes_2_right = nseg_2;

         n_ticks_left  = times_this( n_seg   ) - times_this(     1 ) + 1;
         n_ticks_right = times_this( nt_this ) - times_this( n_seg ) + 1;

         nn_left  = n_seg;
         nn_right = nt_this - n_seg;

         x_seg_left  = nn_left  / n_ticks_left;
         x_seg_right = nn_right / n_ticks_right;

         nseg_1_work = [ nseg_1_work ntimes_1_left ntimes_1_right ];
         nseg_2_work = [ nseg_2_work ntimes_2_left ntimes_2_right ]; 
         xx_work     = [ xx_work  x_seg_left  x_seg_right  ]; 
         nosubs_work = [ nosubs_work 0 0 ];

      else

         % No sub-segmenting of this segment;
         % so just stuff in the beginning, end, mark
         % as "nosubs" so that it will not be done again
         nseg_1_work = [ nseg_1_work nseg_1 ];
         nseg_2_work = [ nseg_2_work nseg_2 ]; 
         xx_work  = [ xx_work x_seg ]; 
         nosubs_work = [ nosubs_work 1 ];

      end

   end

   % Post the segmentations just done:
   nseg_1_vec = nseg_1_work;
   nseg_2_vec = nseg_2_work;
   xx_vec  = xx_work;
   nosubs_vec = nosubs_work;

end
 
n_seg_vec = nseg_2_vec;