function [out_vec]=LoveNonLinearFit(x_Love,y_Love, in_vec)
% Version: 1.0, 19 March 2008
%
% Fits perturbed contour to two Love modes, minimizing sum
% of the squared distances between perturbed contour points and
% the nearest point of the unperturbed ellipse surface
persistent fit_results % previous fit results, for estimates
persistent r_nn % radius used by nearestneighbors.m, wish to use smallest possible
persistent pi_add % to keep rotation angle monotonically increasing
persistent S_Phi_old
% ******************* LOAD VALUES AND DEFINE NEEDED QUANTITIES *************
b = in_vec(1);
m_1 = in_vec(2);
m_2 = in_vec(3);
m_3 = in_vec(4); % flag only, if negative allow a/b to vary
a = in_vec(5);
n_ellipse_points = in_vec(6);
i_force_perp_matrix_creation = in_vec(7); % not used here
i_plot_perp_matrix = in_vec(8); % not used here
i_fit_analysis_prev_use = in_vec(9);
i_plot_fit = in_vec(10);
c = (a^2 - b^2)^0.5;
u_0 = atanh(b/a);
s_1 = linspace(0,2*pi,n_ellipse_points+1); s_1 = s_1(1:n_ellipse_points);
mult_s = 5;
s_2 = linspace(0,2*pi,n_ellipse_points*mult_s+1); s_2 = s_2(1:n_ellipse_points*mult_s);
h_0_squared = 1. ./ (b^2*cos(s_2).^2 + a^2.*sin(s_2).^2);
if i_fit_analysis_prev_use == 0; r_nn = 0.005; end % initial value, increased if too small
r_vec = sqrt(x_Love.^2 + y_Love.^2);
min_r = min(r_vec);
max_r = max(r_vec);
% ******************* DETERMINE ELLIPSE ANGLE PHI *************
S = ellipsefit(x_Love,y_Love);
if i_fit_analysis_prev_use == 0
ellipse_phase = S.Phi;
pi_add = 0;
else
if S.Phi < S_Phi_old; pi_add = pi_add + pi; end % ensure angle monotonically increasing
ellipse_phase = S.Phi + pi_add;
end
S_Phi_old = S.Phi;
% ******************* LOAD ESTIMATES AND FIT *************
if i_fit_analysis_prev_use == 0
out_vec = LoveLinearFit(x_Love,y_Love, [in_vec(1:6), [1 0 0 0]]);
if out_vec(1) < 1
alpha_m1 = out_vec(1); beta_m1 = out_vec(2); % LoveLinearFit used for estimates
alpha_m2 = out_vec(3); beta_m2 = out_vec(4);
else
alpha_m1 = 0.02; beta_m1 = 0.2; % Use this init if LoveLinearFit failed
alpha_m2 = 0.005; beta_m2 = 0.2;
end
else
if fit_results(2) == 0; fit_results(2) = -0.01; end % fitter can get stuck at angle beta = 0
alpha_m1 = fit_results(1); beta_m1 = fit_results(2); % from previous fit
alpha_m2 = fit_results(3); beta_m2 = fit_results(4);
if m_3 < 0; a = fit_results(5); b = fit_results(6); end % a and b not varied unless m_3 negative
end
if m_3 >= 0
x_in = [alpha_m1 beta_m1 alpha_m2 beta_m2 ellipse_phase];
options = optimset('Largescale','off', 'TolFun', 1E-10, 'TolX', 1e-10, 'Display', 'off', ...
'MaxIter', 5, 'MaxFunEvals', 40);
[x_out,resnorm, residual, exitflag, output] = lsqnonlin(@CalcDispConstAB,x_in, [], [], options); % Invoke optimizer
else
x_in = [alpha_m1 beta_m1 alpha_m2 beta_m2 a b ellipse_phase];
options = optimset('Largescale','off', 'TolFun', 1E-10, 'TolX', 1e-10, 'Display', 'off', ...
'MaxIter', 7, 'MaxFunEvals', 50);
[x_out,resnorm, residual, exitflag, output] = lsqnonlin(@CalcDispVaryAB,x_in, [], [], options); % Invoke optimizer
end
% ******************* LOAD FIT RESULTS *************
% alpha_m1 beta_m1 alpha_m2 beta_m2 alpha_m3 beta_m3 phi min_r max_r chi when m_3 positive, m_3 mode not fit
% alpha_m1 beta_m1 alpha_m2 beta_m2 a b phi min_r max_r chi when m_3 negative, now a b are fit
% 1 2 3 4 5 6 7 8 9 10
if m_3 >= 0
if m_2 > 0
fit_results = [x_out(1:4) a b x_out(5)];
else
fit_results = [x_out(1:2) 0 0 a b x_out(5)];
end
else
if m_2 > 0
fit_results = x_out(1:7);
else
fit_results = [x_out(1:2) 0 0 x_out(5:7)];
end
end
chi = sum(residual.^2);
out_vec = [fit_results, min_r, max_r, chi];
if i_plot_fit
cax = newplot;
plot(x_Love_p,y_Love_p,'r'); % rotated data, from last function eval
hold on
plot(x_pert_ellipse, y_pert_ellipse, '-', 'color', 'g'); % perturbed, from last function eval
plot(fit_results(5) * cos(s_1), fit_results(6) * sin(s_1), 'k'); % unperturbed
title(num2str([x_out(1), x_out(3), chi]));
title('Nonlinear fit. Red: data Green: fit Black: unperturbed');
if max(r_vec) > 1
xlim(max(r_vec) * [-1.02 1.02]);
else
xlim([-1.02 1.02]);
end
hold off
axis equal
drawnow
end
% *** DEFINE TWO FUNCTIONS, RETURN DISPLACEMENT FROM UNPERTURBED ELLIPSE AT EACH POINT **********
function F = CalcDispConstAB(x_in) % aspect ratio fixed
alpha_m1_0 = x_in(1);
beta_m1_0 = x_in(2);
alpha_m2_0 = x_in(3);
beta_m2_0 = x_in(4);
ellipse_phase_0 = x_in(5);
i_flag_m_nn = 0;
u_pert_1 = alpha_m1_0 * cos(m_1*(s_2-beta_m1_0)) .* h_0_squared;
if m_2 > 0
u_pert_2 = alpha_m2_0 * cos(m_2*(s_2-beta_m2_0)) .* h_0_squared;
else
u_pert_2 = 0;
end
x_pert_ellipse = c * cosh(u_0 + u_pert_1 + u_pert_2) .* cos(s_2);
y_pert_ellipse = c * sinh(u_0 + u_pert_1 + u_pert_2) .* sin(s_2);
x_Love_p = x_Love * cos(ellipse_phase_0) + y_Love * sin(ellipse_phase_0);
y_Love_p = y_Love * cos(ellipse_phase_0) - x_Love * sin(ellipse_phase_0);
x = [x_Love_p.' ; y_Love_p.' ]; % unperturbed contour
p = [x_pert_ellipse ; y_pert_ellipse]; % perturbed contour
idx = nearestneighbour(x, p, 'n', 1, 'r', r_nn); % nearest point in set of p to each point in x
while min(idx) == 0 % increase r_nn until every point has nearest neighbor
r_nn = r_nn * 2;
idx = nearestneighbour(x, p, 'n', 1, 'r', r_nn);
if r_nn > 2; break, end
end
% Displacement of perturbation: altitude of triangle of x point and nearest 2 points in p
% X point coords: (x_x, y_x) p points coords: (x_p0, y_p0) and (x_p2, y_p2)
for i_ss = 1:n_ellipse_points
x_x = x_Love_p(i_ss); y_x = y_Love_p(i_ss);
switch idx(i_ss)
case 0
if i_flag_m_nn == 0
disp('missing nearest neighbors!');
end
i_flag_m_nn = 1;
x_p0 = x_pert_ellipse(mult_s * n_ellipse_points); y_p0 = y_pert_ellipse(mult_s * n_ellipse_points);
x_p1 = 0; y_p1 = 0; x_p2 = 0; y_p2 = 0;
case 1
x_p0 = x_pert_ellipse(mult_s * n_ellipse_points); y_p0 = y_pert_ellipse(mult_s * n_ellipse_points);
x_p1 = x_pert_ellipse(idx(i_ss)); y_p1 = y_pert_ellipse(idx(i_ss));
x_p2 = x_pert_ellipse(idx(i_ss)+1); y_p2 = y_pert_ellipse(idx(i_ss)+1);
case mult_s * n_ellipse_points % account for endpoint
x_p0 = x_pert_ellipse(idx(i_ss)-1); y_p0 = y_pert_ellipse(idx(i_ss)-1);
x_p1 = x_pert_ellipse(idx(i_ss)); y_p1 = y_pert_ellipse(idx(i_ss));
x_p2 = x_pert_ellipse(1); y_p2 = y_pert_ellipse(1);
otherwise
x_p0 = x_pert_ellipse(idx(i_ss)-1); y_p0 = y_pert_ellipse(idx(i_ss)-1);
x_p1 = x_pert_ellipse(idx(i_ss)); y_p1 = y_pert_ellipse(idx(i_ss));
x_p2 = x_pert_ellipse(idx(i_ss)+1); y_p2 = y_pert_ellipse(idx(i_ss)+1);
end
a_t(i_ss) = sqrt( (x_x - x_p0)^2 + (y_x - y_p0)^2);
b_t(i_ss) = sqrt( (x_x - x_p2)^2 + (y_x - y_p2)^2);
c_t(i_ss) = sqrt( (x_p0 - x_p2)^2 + (y_p0 - y_p2)^2);
end
F = sqrt((a_t + b_t + c_t) .* (b_t + c_t - a_t) .* (c_t + a_t - b_t) .* (a_t + b_t - c_t) ./ (4 * c_t.^2));
end % Function CalcDispConstAB end
% ***************** THIS FUNCTION ALLOWS ASPECT RATIO TO VARY *********************
function F = CalcDispVaryAB(x_in)
alpha_m1_0 = x_in(1);
beta_m1_0 = x_in(2);
alpha_m2_0 = x_in(3);
beta_m2_0 = x_in(4);
a_0 = x_in(5);
b_0 = x_in(6);
ellipse_phase_0 = x_in(7);
i_flag_m_nn = 0;
c_0 = (a_0^2 - b_0^2)^0.5;
u_0_0 = atanh(b_0/a_0);
h_0_squared_0 = 1. ./ (b_0^2*cos(s_2).^2 + a_0^2.*sin(s_2).^2);
u_pert_1 = alpha_m1_0 * cos(m_1*(s_2-beta_m1_0)) .* h_0_squared_0;
if m_2 > 0
u_pert_2 = alpha_m2_0 * cos(m_2*(s_2-beta_m2_0)) .* h_0_squared_0;
else
u_pert_2 = 0;
end
x_pert_ellipse = c_0 * cosh(u_0_0 + u_pert_1 + u_pert_2) .* cos(s_2);
y_pert_ellipse = c_0 * sinh(u_0_0 + u_pert_1 + u_pert_2) .* sin(s_2);
x_Love_p = x_Love * cos(ellipse_phase_0) + y_Love * sin(ellipse_phase_0);
y_Love_p = y_Love * cos(ellipse_phase_0) - x_Love * sin(ellipse_phase_0);
x = [x_Love_p.' ; y_Love_p.' ]; % unperturbed contour
p = [x_pert_ellipse ; y_pert_ellipse]; % perturbed contour
idx = nearestneighbour(x, p, 'n', 1, 'r', r_nn); % nearest point in set of p to each point in x
while min(idx) == 0 % increase r_nn until every point has nearest neighbor
r_nn = r_nn * 2;
idx = nearestneighbour(x, p, 'n', 1, 'r', r_nn);
if r_nn > 2; break, end
end
% Displacement of perturbation: altitude of triangle of x point and nearest 2 points in p
% X point coords: (x_x, y_x) p points coords: (x_p0, y_p0) and (x_p2, y_p2)
for i_ss = 1:n_ellipse_points
x_x = x_Love_p(i_ss); y_x = y_Love_p(i_ss);
switch idx(i_ss)
case 0
if i_flag_m_nn == 0
disp('missing nearest neighbors!');
end
i_flag_m_nn = 1;
x_p0 = x_pert_ellipse(mult_s * n_ellipse_points); y_p0 = y_pert_ellipse(mult_s * n_ellipse_points);
x_p1 = 0; y_p1 = 0; x_p2 = 0; y_p2 = 0;
case 1
x_p0 = x_pert_ellipse(mult_s * n_ellipse_points); y_p0 = y_pert_ellipse(mult_s * n_ellipse_points);
x_p1 = x_pert_ellipse(idx(i_ss)); y_p1 = y_pert_ellipse(idx(i_ss));
x_p2 = x_pert_ellipse(idx(i_ss)+1); y_p2 = y_pert_ellipse(idx(i_ss)+1);
case mult_s * n_ellipse_points % account for endpoint
x_p0 = x_pert_ellipse(idx(i_ss)-1); y_p0 = y_pert_ellipse(idx(i_ss)-1);
x_p1 = x_pert_ellipse(idx(i_ss)); y_p1 = y_pert_ellipse(idx(i_ss));
x_p2 = x_pert_ellipse(1); y_p2 = y_pert_ellipse(1);
otherwise
x_p0 = x_pert_ellipse(idx(i_ss)-1); y_p0 = y_pert_ellipse(idx(i_ss)-1);
x_p1 = x_pert_ellipse(idx(i_ss)); y_p1 = y_pert_ellipse(idx(i_ss));
x_p2 = x_pert_ellipse(idx(i_ss)+1); y_p2 = y_pert_ellipse(idx(i_ss)+1);
end
a_t(i_ss) = sqrt( (x_x - x_p0)^2 + (y_x - y_p0)^2);
b_t(i_ss) = sqrt( (x_x - x_p2)^2 + (y_x - y_p2)^2);
c_t(i_ss) = sqrt( (x_p0 - x_p2)^2 + (y_p0 - y_p2)^2);
end
F = sqrt((a_t + b_t + c_t) .* (b_t + c_t - a_t) .* (c_t + a_t - b_t) .* (a_t + b_t - c_t) ./ (4 * c_t.^2));
F = real(F); % remove epsilon imaginary part which (rarely) appears erroneously
end % Function CalcDispVaryAB end
end % Function LoveNonLinearFit end
