function fig9_18

global thetad xd yd

% generate some data points
n=50;
thetad=2*pi*rand(n,1);
xd=zeros(n,1);  yd=xd;
for i=1:n
    a(i)=2+0.4*(rand-0.5); b(i)=1+0.4*(rand-0.5); phi(i)=pi*(1+0.3*(rand-0.5))/3;
    xd(i)=a(i)*cos(phi(i))*cos(thetad(i))-b(i)*sin(phi(i))*sin(thetad(i));
    yd(i)=b(i)*cos(phi(i))*sin(thetad(i))+a(i)*sin(phi(i))*cos(thetad(i));
end

% circle_data=[cxd;yd; a ;b ;phi];
% save('/Users/mark/Desktop/cicrle_data.txt','circle_data','-ascii');

% circle_data=importdata('cicrle_data.txt');
% xd=circle_data(1,:);
% yd=circle_data(2,:);
% a=circle_data(3,1); b=circle_data(4,1); phi=circle_data(5,1);

x=[1.0, 0.5, 0]';
tol=1e-6;
xpoints=x(1);
ypoints=x(2);
f=F(x);
ff=norm(f,2)^2;
D=zeros(3,3);
mu=10;
alpha=3;
beta=2;
counter=1;
err=10*tol;
while err>tol %& counter<200
    counter=counter+1;
    xx=x;
    ff0=ff;

    J=Jac(x);
    JJ=J'*J;
    gk=2*J'*f;
    D(1,1)=2*JJ(1,1);  D(2,2)=2*JJ(2,2);  D(3,3)=2*JJ(3,3);
    S=2*JJ+mu*D;
    dx=-inv(S)*gk;
    x=xx+dx;
    f=F(x);
    ff=norm(f,2)^2;

    if ff<ff0
        mu=mu/beta;
        xpoints=[xpoints, x(1)];
        ypoints=[ypoints, x(2)];
        err=norm(-inv(S)*gk./x,inf);
    else
        mu=alpha*mu; x=xx;  ff=ff0;
    end

end
counter

clf
hold on
% get(gcf)
set(gcf,'Position', [12 950 421 395])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

plot(xd,yd,'bo','MarkerSize',8,'LineWidth',2)
axis square
box on
nt=100;
t=linspace(0,2*pi,nt);
a=x(1)
b=x(2)
phi=x(3)
for i=1:nt
    xc(i)=a*cos(phi)*cos(t(i))-b*sin(phi)*sin(t(i));
    yc(i)=b*cos(phi)*sin(t(i))+a*sin(phi)*cos(t(i));
end
plot(xc,yc,'r','LineWidth',2)
plot([0 0],[-2 2],'k','LineWidth',0.5)
plot([-2 2],[0 0],'k','LineWidth',0.5)
xlabel('x-axis')
ylabel('y-axis')
set(gca,'FontSize',16,'FontWeight','bold')
hold off

%exportgraphics(gcf,'/Users/mark/Desktop/ellipse2.eps')

% F=dot(f,f)
function f=F(x)
global thetad xd yd
a=x(1); b=x(2); phi=x(3);
n=length(xd);
f=zeros(n,1);
for i=1:n
    xf=a*cos(phi)*cos(thetad(i))-b*sin(phi)*sin(thetad(i)) -xd(i);
    yf=b*cos(phi)*sin(thetad(i))+a*sin(phi)*cos(thetad(i)) -yd(i);
    f(i)=xf^2+yf^2;
end

% jacobian
function g=Jac(x)
global thetad xd yd
a=x(1); b=x(2); phi=x(3);
n=length(xd);
g=zeros(n,3);
for i=1:n
    xf=a*cos(phi)*cos(thetad(i))-b*sin(phi)*sin(thetad(i)) -xd(i);
    yf=b*cos(phi)*sin(thetad(i))+a*sin(phi)*cos(thetad(i)) -yd(i);
    g(i,1)=2*xf*cos(phi)*cos(thetad(i))+2*yf*sin(phi)*cos(thetad(i));
    g(i,2)=-2*xf*sin(phi)*sin(thetad(i))+2*yf*cos(phi)*sin(thetad(i));

    g31=2*xf*(-a*sin(phi)*cos(thetad(i))-b*cos(phi)*sin(thetad(i)));
    g31=g31+2*yf*(-b*sin(phi)*sin(thetad(i))+a*cos(phi)*cos(thetad(i)));
    g(i,3)=g31;
end


















