function fig8_3

global xd yd

% plot error functionS for least square error
% this takes maybe 15 sec to run

tic

%  generate data
% n=50;
% xd=rand(1,n);
% yd=zeros(1,n);
% a1=3; a2=1;
% for i=1:n
%     yd(i)=a1*xd(i)+a2+(rand(1)-0.5)*(a1*xd(i)+a2);
% end
% surf_data=[xd;yd];
% save('/Users/mark/Desktop/surf_data.txt','surf_data','-ascii');

surf_data=importdata('surf_data.txt');
xd=surf_data(1,:);
yd=surf_data(2,:);

figure(1)
clf
% get(gcf)
set(gcf,'Position', [5 925 560 420]);
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)

% E_infty surface
v2a=-2; v2b=8; v1a=-4; v1b=6;
%n0=4000;
n0=1100;
v2=linspace(v2a,v2b,n0);
v1=linspace(v1a,v1b,n0);
[A1,A2]=meshgrid(v1,v2);
for i=1:n0
    for j=1:n0
        E(i,j)=Einf(v1(j),v2(i));
    end
end
surfc(A1,A2,E)
shading interp
colormap copper
view(-42, 16)
xlabel('v_1')
ylabel('v_2')
zlabel('E_\infty')
set(gca,'FontSize',20,'FontWeight','bold')

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

%fminsearch(@EM, a, optimset('TolFun',eps,'TolX',eps))
vinfmim = fminsearch(@EEinf, [2, 2])

% E_infty contours
figure(2)
clf
% get(gcf)
set(gcf,'Position', [7 429 560 420])
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)

n0=2000;
v2=linspace(v2a,v2b,n0);
v1=linspace(v1a,v1b,n0);
[A1,A2]=meshgrid(v1,v2);
for i=1:n0
    for j=1:n0
        E(i,j)=Einf(v1(j),v2(i));
    end
end

contour(A1,A2,E,12,'b','LineWidth',1.6)
hold on
plot(vinfmim(1),vinfmim(2),'r*','MarkerSize',14,'LineWidth',2)
xlabel('v_1-axis')
ylabel('v_2-axis')
set(gca,'ytick',[-4 -2 0 2 4 6 8])
set(gca,'FontSize',20,'FontWeight','bold')

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

% E_2 surface
figure(3)
clf
% get(gcf)
set(gcf,'Position', [569 925 560 420]);
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)

n0=1100;
v2=linspace(v2a,v2b,n0);
v1=linspace(v1a,v1b,n0);
[A1,A2]=meshgrid(v1,v2);
for i=1:n0
    for j=1:n0
        E(i,j)=Einf(v1(j),v2(i));
    end
end

rM=0.7;
[A1,A2,EE2]=E2_surf(rM);

surfc(A1,A2,EE2)
shading interp
colormap copper
view(-42, 16)
xlabel('v_1')
ylabel('v_2')
zlabel('E_2')
set(gca,'FontSize',20,'FontWeight','bold')

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

% E_2 contours
figure(4)
clf
% get(gcf)
set(gcf,'Position', [568 429 560 420])
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)

v2a=-2; v2b=8; v1a=-4; v1b=6;
n0=100;
v2=linspace(v2a,v2b,n0);
v1=linspace(v1a,v1b,n0);
[A1,A2]=meshgrid(v1,v2);
for i=1:n0
    for j=1:n0
        EE2(i,j)=E2(v1(j),v2(i));
    end
end

%fminsearch(@EM, a, optimset('TolFun',eps,'TolX',eps))
v2mim = fminsearch(@E22, [2, 2])

contour(A1,A2,EE2,[30 60 115 190 320 480 700 900 1200 1550 1900 2400],'b','LineWidth',1.6)
hold on
plot(v2mim(1),v2mim(2),'r*','MarkerSize',14,'LineWidth',2)
xlabel('v_1-axis')
ylabel('v_2-axis')
set(gca,'ytick',[-4 -2 0 2 4 6 8])
set(gca,'FontSize',20,'FontWeight','bold')
%exportgraphics(gcf,'/Users/mark/Desktop/cont2.eps')
toc


function g=Einf(v1,v2)
global xd yd
g=norm(v1+v2*xd-yd,inf);

function g=EEinf(v)
global xd yd
g=norm(v(1)+v(2)*xd-yd,inf);


function g=E2(v1,v2)
global xd yd
g=norm(v1+v2*xd-yd,2)^2;

function g=E22(v)
global xd yd
g=norm(v(1)+v(2)*xd-yd,2)^2;

function [A1,A2,E]=E2_surf(rM)
global xd yd
a=length(xd); b=2*sum(xd); c=dot(xd,xd); d=-2*sum(yd); e=-2*dot(xd,yd); f=dot(yd,yd);
theta=0.5*abs(atan(b/(a-c)));
theta/pi
if b>=0 & a>=c
    ct=cos(theta);
    st=sin(theta);
elseif b>=0 & a < c
    ct=-cos(theta);
    st=sin(theta);
elseif b < 0 & a>=c
    ct=cos(theta);
    st=-sin(theta);
elseif b < 0 & a < c
    ct=-cos(theta);
    st=-sin(theta);
end

%%% transform to E=A*xp^2+B*yp^2+fb
%%% where x=x0+xp*ct-yp*ct, y=y0+xp*ct+yp*st
x0=(2*c*d-b*e)/(b^2-4*a*c);
y0=(2*a*e-b*d)/(b^2-4*a*c);
A=a*ct^2+b*ct*st+c*st^2;
B=a*st^2-b*ct*st+c*ct^2;
fb=f+a*x0^2+b*x0*y0+c*y0^2+d*x0+e*y0;
n0=100;
% rM=1/3;
r=linspace(0.0001,rM,n0);
ang=linspace(0,2*pi,n0);
for ir=1:n0
    for it=1:n0
        xp=sqrt(B)*r(ir)*cos(ang(it));
        yp=sqrt(A)*r(ir)*sin(ang(it));
        A1(ir,it)=x0+xp*ct-yp*ct;
        A2(ir,it)=y0+xp*ct+yp*st;
        %        E(ir,it)=a*A1(ir,it)^2+b*A1(ir,it)*A2(ir,it)+c*A2(ir,it)^2+d*A1(ir,it)+e*A2(ir,it)+f;
        E(ir,it)=A*xp^2+B*yp^2+fb;
    end
end














