out_vec=PiecewiseIntegration(in_vec)

function out_vec=PiecewiseIntegration(in_vec)
%  Version: 1.0, 19 March 2008
%
%  Integrates in time steps, allowing observation of evolution and controlled
%  halt if run is about to hang

RelTolValue             = in_vec(1); 
AbsTolValue             = in_vec(2);
a                       = in_vec(3);
b                       = in_vec(4);
zeta                    = in_vec(5);
m_1                     = in_vec(6);
alpha_m                 = in_vec(7);
beta_m                  = in_vec(8);
sim_end_time            = in_vec(9);
sim_steps               = in_vec(10);
frame_steps             = in_vec(11);
safe_max_del_r_mult     = in_vec(12);
safe_max_min_r_mult     = in_vec(13);
safe_max_max_r          = in_vec(14);
safe_mult_time          = in_vec(15);
safe_min_r              = in_vec(16);
i_resample              = in_vec(17);
t_resample              = in_vec(18);
t_resample_start        = in_vec(19);
i_plot_safe_calc        = in_vec(20);
n_ellipse_points        = in_vec(21);
n_time_points_sim       = in_vec(22);
Z0_start                = in_vec(23:end);
Z0_start                = Z0_start';
        
%   ************   FLAG INITIALIZATIONS   *********************
i_flag_end              = 0;    % use to end when safe value exceeded
pi_add                  = 0;    % keeps ellipse angle monotonically increasing as ellipse rotation in CCW direction
ellipse_fit_old         = 0;    % previous value of ellipse angle

s = linspace(0,2*pi,n_ellipse_points+1);  s = s(1:n_ellipse_points);
options = odeset('RelTol',RelTolValue,'AbsTol',AbsTolValue);
safe_max_min_r = safe_max_min_r_mult * b;

for i_sim_step = 1:double(int16(sim_end_time / sim_steps - 0.01))
    
    time_start  = (i_sim_step-1) * sim_steps;                       % determine start and end time
    time_end    = i_sim_step * sim_steps;
    if time_end > sim_end_time || sim_end_time - time_end < frame_steps
        time_end = sim_end_time;
        i_flag_end = 1;
    end
    t_0 = cputime;                                                  % perform integration over time step
    [T_step Y_Love_step] = ode113(@(t,Z) dZdt(Z,s,zeta),time_start:frame_steps:time_end,Z0_start,options);  
    time_needed = cputime - t_0;
           
    if i_sim_step == 1                                          % if first step, define baseline quantities to check for halt conditions
        T = T_step;
        Y_Love = Y_Love_step;
        standard_time = time_needed;
        del_r_vec = sqrt( (real(Y_Love(end,1:end-1)) - real(Y_Love(end,2:end)) ).^2 + (imag(Y_Love(end,1:end-1)) - imag(Y_Love(end,2:end)) ).^2);
        del_r_mean = sum(del_r_vec) / length(del_r_vec);
        max_spacing_init = del_r_mean;
        min_radius_init = b;  
        disp(['standard time: ', num2str(standard_time)]);
    else                                                        % add new step onto existing ones
        T      = [T.', T_step(2:end).'].';
        Y_Love  = [Y_Love.', Y_Love_step(2:end,:).'].';
    end
                                                                % quantities for various tests of nodes
    del_r_vec = sqrt( (real(Y_Love(end,1:end-1)) - real(Y_Love(end,2:end)) ).^2 + (imag(Y_Love(end,1:end-1)) - imag(Y_Love(end,2:end)) ).^2);
    del_r_mean = sum(del_r_vec) / length(del_r_vec);
    r_vec = sqrt( (real(Y_Love(end,:))).^2 + (imag(Y_Love(end,1:end))).^2);
            
    if (i_resample == 1 && rem(T(end),t_resample) == 0 && T(end) >= t_resample_start) 
        resam_str = ' resample';
    else
        resam_str = '';
    end
            
    disp(['i_sim_step: ',num2str(sim_steps*i_sim_step, '%3.2f'), '  time: ',num2str(time_needed / standard_time, '%4.3f'), ...
    '  mean del r: ', num2str(del_r_mean / max_spacing_init, '%4.3f'), ' min r: ',  num2str(min(r_vec) , '%4.3f'), ' max r: ',  ...
    num2str(max(r_vec), '%4.3f' ), resam_str  ]);
    if del_r_mean > max_spacing_init * safe_max_del_r_mult || min(r_vec) < safe_min_r || ...
        time_needed > standard_time * safe_mult_time || min(r_vec) > safe_max_min_r || max(r_vec) > safe_max_max_r;
        i_flag_end = 1;
    end
    if i_flag_end, break, end;

    if (i_resample == 1 && rem(T(end),t_resample) == 0 && T(end) >= t_resample_start) 
        x = real(Y_Love(end,:));  x = [x,x(1)];
        y = imag(Y_Love(end,:));  y = [y,y(1)];
        N = n_ellipse_points + 1;
        p = [x',y'];
        q = curvspace(p,N);
        x = q(1:end-1,1); y = q(1:end-1,2);
        Z0_start = x + i * y;
    else
        Z0_start = Y_Love(end,:);
    end   
    
    if i_plot_safe_calc == 1
        hold off
        x_Love = real(Z0_start); y_Love = imag(Z0_start);
        S=ellipsefit(x_Love,y_Love);  
        
        if i_sim_step > 1 && S.Phi < ellipse_fit_old
            pi_add = pi_add + pi;
        end
        
        x_Love_p = x_Love * cos(S.Phi + pi_add) + y_Love * sin(S.Phi + pi_add);
        y_Love_p = y_Love * cos(S.Phi + pi_add) - x_Love * sin(S.Phi + pi_add);
        ellipse_fit_old = S.Phi;
        x_ellipse  = a * cos(s);    y_ellipse  = b * sin(s);
        area_init = polyarea(x_ellipse, y_ellipse);
        area_now = polyarea(x_Love, y_Love);
        plot(x_Love_p, y_Love_p, '-+','color', 'b', 'MarkerSize', 4);
        hold on
        plot(real(Y_Love(1,:)),imag(Y_Love(1,:)),'-', 'color','k');
        plot(x_ellipse,y_ellipse,'r');
        if b<0.5 ; axis equal; end;
        if max(r_vec) > 1
            xlim(max(r_vec) * [-1.02 1.02]);
        else
            xlim([-1.02 1.02]);
        end
                
        ylim([min(y_Love_p)*1.05 1.05*max(y_Love_p)]);
        title(['t= ',num2str(T(end)), '  m: ', num2str(m_1), '  a/b: ', num2str(1/b), '  \alpha_m^0: ', ...
            num2str(alpha_m), '  \beta_m^0: ', num2str(beta_m), '   del area: ', num2str(area_now / area_init)]);
        drawnow
        hold off
    end % i_plot_safe_calc == 1
end
        safe_time_end = time_end;
        disp(['final time end: ', num2str(safe_time_end)]);
        n_time_points_safe = length(T);
    
        if safe_time_end ~= sim_end_time || n_time_points_safe ~= n_time_points_sim
            disp('simulation ended early');
        end
     out_vec = [T Y_Love];