function [M,S, R_mean] = Apical_fit(Xo, Yo)


N = length(Xo);
Ro = sqrt(Xo.^2 + Yo.^2); 

Theta = atan2(Yo, Xo);
%Theta = linspace(-pi/2, pi/2, N)' + pi/2;

X = linspace(-5, 5, N)';
R_mean = mean(Ro(1:20));

R = @(P) ones(N, 1)*R_mean - P(1)*exp(-((X - P(3))/P(2)).^2);
F_fun = @(P) sum((Xo - R(P).*cos(Theta)).^2 + (Yo - R(P).*sin(Theta)).^2);
P_fit = fminsearch(F_fun, [5 1 0]);


M = P_fit(1);
S = P_fit(2);
plot(R(P_fit).*cos(Theta), R(P_fit).*sin(Theta), 'r');

end