function [ change_point, odds_21, log_prob, log_prob_noseg ] = ...
    find_change( photon_times, t_0, t_n )
%
% Find most probable two-rate model for Poisson arrival time data, 
% based on Bayesian analysis.
%
% Input:  photon_times -- photon arrival times
%                         (Note: These must be microseconds, not seconds,
%                         because the time values correspond to the 
%                         clock rate at which the data are sampled.)
%         t_0          -- time just previous to photon_times(1)
%         t_n          -- time just after last time in photon_times
%
% Output: change_point -- index of "photon_times" which provides the maximum 
%                         likelihood segmented model (that is, with one   
%                         Poisson rate to the left of 
%                               photon_times(change_point) 
%                         and another to the right
%              odds_21 -- odds ratio: 2 unequal rates / 1 rate
%             log_prob -- log probability of segmented model, as a  
%                         function of changepoint 
%       log_prob_noseg -- log prob of nonsegmented model
%--------------------------------------------------------------------------------------

dt_half = 0.5 * diff( photon_times );
n_ph = length( photon_times ); % Number of photons

min_time = photon_times(    1 );
max_time = photon_times( n_ph );

t_left  = 0.5 * ( t_0 + min_time );
t_right = 0.5 * ( max_time + t_n );

% Number of microsecond "ticks" in the whole (extended) interval:
n_ticks = t_right - t_left + 1; 

%------------------------------------------------------------------
%    Evaluate log-probability of the unsegmented model:
%------------------------------------------------------------------

log_prob_noseg = gammaln( n_ph + 1  ) + ...
                 gammaln( n_ticks - n_ph + 1 ) - ...
                 gammaln( n_ticks + 2);

%------------------------------------------------------------------
% Evaluate the log-probability of the segmented model as a 
% function of changepoint; find optimum changepoint.
%------------------------------------------------------------------

% Array of trial changepoints:
n_1  = (1: n_ph - 1)';
n_2 = n_ph - n_1;

m_1 = photon_times( n_1 ) + dt_half( n_1 ) - t_left;
m_2 = n_ticks - m_1;

log_prob = - 1.e55 * ones( n_ph, 1 ); % mark all points as invalid
arg_1 = m_1 - n_1 + 1;
arg_2 = m_2 - n_2 + 1;
ii = find( arg_1 > 0 & arg_2 > 0 ); % indices of valid points
log_prob(ii) = gammaln( n_1(ii)+1 ) + gammaln( arg_1(ii) ) - gammaln( m_1(ii)+2 );
log_prob(ii) = log_prob(ii) + gammaln( n_2(ii)+1 ) + gammaln( arg_2(ii) ) - ...
              gammaln( m_2(ii) + 2 );
[ max_log, change_point ] = max( log_prob(ii) );

%------------------------------------------------------------------
% Compute odds ratio: prob(seg) / prob(no_seg) 
%------------------------------------------------------------------
odds_21 = sum( exp( log_prob - log_prob_noseg ) );

if ~isfinite( odds_21 )
   odds_21 = 1000000; 
end