function [out_vec]=LoveLinearFit(x_Love,y_Love, in_vec)
% Version: 1.0, 19 March 2008
%
% Fits perturbed contour to 1, 2 or 3 Love modes
% The perturbation (q_vec) is determined from the deviations
% along local perpendiculars to the unperturbed ellipse surface.
% It is then fit to sinusoids with m specified, see Eq. (15b) of paper
%
% May need to adjust intercept_lin_in and _out parameters from default
% values of 0.5 depending on aspect ratio of base ellipse. Use perp_matrix
% plot to identify problem
persistent perp_matrix
persistent fit_results % previous fit results, for estimates
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); % if non-zero, fit to two modes
m_3 = in_vec(4); % if non-zero, fit to three modes
a = in_vec(5); % a=1, ignored here (used in LoveNonlinearFit)
n_ellipse_points = in_vec(6);
i_force_perp_matrix_creation = in_vec(7);
i_plot_perp_matrix = in_vec(8);
i_fit_analysis_prev_use = in_vec(9);
i_plot_fit = in_vec(10);
c = (a^2 - b^2)^0.5;
u_0 = acosh(a/c);
s = linspace(0,2*pi,n_ellipse_points+1);
s = s(1:n_ellipse_points);
h_0_squared = 1. ./ (b^2*cos(s).^2 + a^2.*sin(s).^2);
intercept_len_in = 0.5 * b ; % perturbation boundaries as u + delta u
intercept_len_out = 0.5 ; % value of 0.5 may need to be adjusted by hand
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; % save angle for next run
x_Love_p = x_Love * cos(ellipse_phase) + y_Love * sin(ellipse_phase);
y_Love_p = y_Love * cos(ellipse_phase) - x_Love * sin(ellipse_phase);
x_ellipse = a * cos(s);
y_ellipse = b * sin(s);
x_pert = zeros(1,n_ellipse_points); y_pert = zeros(1,n_ellipse_points);
u_pert = zeros(1,n_ellipse_points); u_pert_y = zeros(1,n_ellipse_points);
x_0_intercept = zeros(1); y_0_intercept = zeros(1);
q_vec = zeros(1,n_ellipse_points);
% ******************* CREATE PERP MATRIX *************
if exist('perp_matrix') == 0 | i_force_perp_matrix_creation == 1 % slightly time consuming to make
scale_vec = zeros(n_ellipse_points); % so only constructed once
perp_matrix = zeros(2,3,n_ellipse_points);
x_tan_vec = -1 * a * sin(s) ./ sqrt(b^2 * cos(s).^2 + a^2 * sin(s).^2);
y_tan_vec = b * cos(s) ./ sqrt(b^2 * cos(s).^2 + a^2 * sin(s).^2);
theta_tan_vec = atan2(y_tan_vec, x_tan_vec);
theta_normal_vec = theta_tan_vec - pi / 2.;
perp_matrix(1,2,:) = x_ellipse; % perp matrix gives 3 point line intercept at each circumference point
perp_matrix(2,2,:) = y_ellipse; % only the two end points used with intersections routine
perp_matrix(1,3,:) = perp_matrix(1,2,:) + reshape(cos(theta_normal_vec),1,1,n_ellipse_points) * intercept_len_out;
perp_matrix(1,1,:) = perp_matrix(1,2,:) - reshape(cos(theta_normal_vec),1,1,n_ellipse_points) * intercept_len_in;
perp_matrix(2,3,:) = perp_matrix(2,2,:) + reshape(sin(theta_normal_vec),1,1,n_ellipse_points) * intercept_len_out;
perp_matrix(2,1,:) = perp_matrix(2,2,:) - reshape(sin(theta_normal_vec),1,1,n_ellipse_points) * intercept_len_in;
end
if i_plot_perp_matrix ==1 % diagnostic output, plot of contours, perp lines, intersection points
plot(x_ellipse, y_ellipse, 'b'); % can make sure intercept_len_out and _in are correct lengths
hold on
plot(x_Love_p, y_Love_p, 'g');
end
% ******************* DETERMINE INTERCEPTS *************
i_flag_no_intercept = 0;
for i_p = 1:n_ellipse_points
[x_0_intercept, y_0_intercept] = intersections(perp_matrix(1,[1,3],i_p),perp_matrix(2,[1,3],i_p), ...
[x_Love_p.',x_Love_p(1)],[y_Love_p.',y_Love_p(1)],'false');
if isempty([x_0_intercept, y_0_intercept]) == 0 & length(x_0_intercept) == 1 % intercepts OK
x_pert(i_p) = x_0_intercept;
y_pert(i_p) = y_0_intercept;
else
i_flag_no_intercept = 1; % problem with intercepts - mode amplitude may be too large
break;
end
u_pert(i_p) = acosh(x_pert(i_p) / c / cos(s(i_p)));
u_pert_y(i_p) = asinh(y_pert(i_p) / c / sin(s(i_p)));
u_pert(find(abs(x_pert) < 0.2)) = u_pert_y(find(abs(x_pert) < 0.2)); % avoid singularity when x near zero
if i_plot_perp_matrix == 1
plot(perp_matrix(1,:,i_p), perp_matrix(2,:,i_p), 'g');
plot(x_pert(i_p), y_pert(i_p), 'o');
if i_p == n_ellipse_points
hold off
axis equal
end;
end;
end;
% ******************* FIT MODES TO PERTURBED CONTOUR *************
if i_flag_no_intercept == 0
q_vec = real((u_pert - u_0) ./ h_0_squared); % See Eq. (14) of paper
q_vec = q_vec - mean(q_vec); % Remove m=0 component, if any
if i_fit_analysis_prev_use == 1 % use previous fit results as estimate if available and option is selected
alpha_m1 = fit_results(1);
beta_m1 = fit_results(2);
alpha_m2 = fit_results(3);
beta_m2 = fit_results(4);
alpha_m3 = fit_results(5);
beta_m3 = fit_results(6);
else % create rough estimates when first fitting
alpha_m1 = max(q_vec);
max_first_peak = max(q_vec(1:int16(n_ellipse_points/m_1)));
peak_index = find(q_vec(1:int16(n_ellipse_points/m_1)) == max_first_peak);
beta_m1 = peak_index * 2 * pi / n_ellipse_points;
if beta_m1 > pi/m_1; beta_m1 = beta_m1 - 2*pi/m_1; end
alpha_m2 = alpha_m1 * 0.1;
beta_m2 = beta_m1 * 0.1;
alpha_m3 = alpha_m1 * 0.1;
beta_m3 = beta_m1 * 0.1;
end
options = optimset('Largescale','on', 'TolFun', 1E-10, 'TolX', 1e-10, 'Display', 'off', ...
'LevenbergMarquardt', 'off', 'MaxIter', 4, 'MaxFunEvals', 30);
if m_2 == 0 % one mode
fit_function = @(x) (x(1) * cos(m_1 .* (s - x(2)))) - q_vec;
x_in = [alpha_m1 beta_m1];
lb = [-1 -2*pi/m_1]; ub = [1 4*pi/m_1];
else
if m_3 == 0 % two modes
fit_function = @(x) (x(1) * cos(m_1 .* (s - x(2)))) + (x(3) * cos(m_2 .* (s - x(4)))) - q_vec;
x_in = [alpha_m1 beta_m1 alpha_m2 beta_m2];
lb = [-1 -2*pi/m_1 -1 -2*pi/m_2]; ub = [1 2*pi/m_1 1 2*pi/m_2];
else % three modes
fit_function = @(x) (x(1) * cos(m_1 .* (s - x(2)))) + (x(3) * cos(m_2 .* (s - x(4)))) + (x(5) * cos(m_3 .* (s - x(6)))) - q_vec;
x_in = [alpha_m1 beta_m1 alpha_m2 beta_m2 alpha_m3 beta_m3];
lb = [-1 -2*pi/m_1 -1 -2*pi/m_2 -1 -2*pi/m_3]; ub = [1 2*pi/m_1 1 2*pi/m_2 1 2*pi/m_3];
end
end
[x_out, resnorm, residual, exitflag, output] = lsqnonlin(fit_function,x_in, lb, ub, options); % Invoke optimizer
% ******************* LOAD FIT RESULTS *************
fit_results = x_out;
if length(x_out) < 6
fit_results = [fit_results zeros(1,6 - length(x_out))];
end
chi = sum(residual.^2);
else
fit_results = [1 0 0 0 0 0]; % failed fit
chi = 0.;
end
out_vec = [fit_results, ellipse_phase, min_r, max_r, chi];
% ******************* DIAGNOSTIC PLOT OF FIT AND RESIDUALS *************
if i_plot_fit & i_flag_no_intercept == 0
plot(s, q_vec,'r')
hold on
q_fit = (fit_results(1) * cos(m_1 .* (s - fit_results(2)))) ...
+ (fit_results(3) * cos(m_2 .* (s - fit_results(4)))) ...
+ (fit_results(5) * cos(m_3 .* (s - fit_results(6))));
plot(s, q_fit, 'g');
plot(s, residual, 'k');
hold off
title('Linear fit to q(\eta). Red: data Green: fit Black: residuals');
xlabel('\eta');
ylabel('q');
drawnow
end
