[Z0_Love_out]=FitLoveEnd(Z0_Love,m_1, trim_ang, trim_ang_2, i_plot)

function [Z0_Love_out]=FitLoveEnd(Z0_Love,m_1, trim_ang, trim_ang_2, i_plot)
%  Version: 1.0, 19 March 2008
%
%trim_ang = pi/8;
%trim_ang_2 = pi/4;
%i_plot = 0;
x_trim_bndy = 0.7;
n_ellipse_points = length(Z0_Love);
s = linspace(0,2*pi,n_ellipse_points+1);
s = s(1:n_ellipse_points);
x_Love = real(Z0_Love); y_Love = imag(Z0_Love);
      
i_trim = [find(abs(s) > trim_ang & abs(s) < pi - trim_ang) ...
    find(abs(s) > pi + trim_ang & abs(s) < 2 * pi - trim_ang )];
x_trim = x_Love(i_trim); y_trim = y_Love(i_trim);
       
i_fit = [find(abs(s) > trim_ang & abs(s) < trim_ang_2) ...
    find(abs(s) > pi - trim_ang_2 & abs(s) < pi - trim_ang) ...
    find(abs(s) > pi + trim_ang & abs(s) < pi + trim_ang_2 ) ...
    find(abs(s) > 2 * pi - trim_ang_2 & abs(s) < 2 * pi - trim_ang)];
                                              % add tip if they are unique - no knots
                                              
[x_0_intercept, y_0_intercept] = intersections([-5 5], [0 0], [x_Love.',x_Love(1)],[y_Love.',y_Love(1)],'false');
if length(x_0_intercept) == 2
    i_fit = [i_fit 1]; 
    i_fit = [i_fit length(Z0_Love)/2+1];
end
        
x_fit = x_Love(i_fit); y_fit = y_Love(i_fit);
t_fit = atan2(y_fit, x_fit);
i_rem = find(t_fit < -pi/2);
t_fit(i_rem) = t_fit(i_rem) + 2*pi;         
        
[t_fit_sort i_sort] = sort(t_fit);
x_fit_sort = x_fit(i_sort);
y_fit_sort = y_fit(i_sort);
x_new = x_trim.' ;
y_new = y_trim.' ;
        
for i_add = 1:2
    if i_add == 1;  i_add = find(x_fit_sort > 0); end;
    if i_add == 2;  i_add = find(x_fit_sort < 0); end;
    pp = csapi(y_fit_sort(i_add), x_fit_sort(i_add));
       
    y_vec = linspace(-1*min(abs(y_trim)), min(abs(y_trim)), (n_ellipse_points - length(x_trim))/2 + 2);
    y_vec = y_vec(2:end-1);
    yy = ppval(pp, y_vec);        
    axis equal
    fprime_1 = fnder(pp,1);
    fprime_2 = fnder(pp,2);
    fprime_3 = fnder(pp,3);
    x_new = [x_new yy];   y_new = [y_new y_vec];
end
        
theta_new = atan2(y_new, x_new);
[theta_sort i_theta_sort] = sort(theta_new);
x_new = x_new(i_theta_sort);
y_new = y_new(i_theta_sort);
if i_plot == 1
    plot(x_Love,y_Love, 'b');
    xlim([-1.1 1.1]);
    hold on            
    plot(x_trim, y_trim, 'o', 'MarkerSize', 5);
    plot(x_fit_sort, y_fit_sort, 'r');
    plot(x_new, y_new, '-+', 'color', 'k');
    axis equal
    hold off                       
end
       
Z0_Love_out = (x_new + i * y_new).';