function fig9_4

figure(2)
clf
% get(gcf)
set(gcf,'Position', [2 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)

%%% the following code changes coordinates to align with elliptic contours;
%%% it is longer than it should be because it is written for the general
%%% quadratic form E=ax^2+bxy+cy^2+dx+ey+f
a11=3; a12=-2; a22=4; b1=1; b2=-1;
a=a11; b=a12/2; c=a22; d=-b1; e=-b2; f=0;

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=30;
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

Emax=fb+A*B*rM^2
Emin=fb

surfc(A1,A2,E)
view(60,20)

xlabel('x')
ylabel('y')
zlabel('F(x,y)')
set(gca,'ytick',[-0.5 0 0.5])
set(gca,'xtick',[-0.4 0 0.4 0.8])
set(gca,'FontSize',16,'FontWeight','bold')
box on
H=gca;
H.LineWidth=1;

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

figure(1)
clf
% get(gcf)
set(gcf,'Position', [4 414 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)

%contour(A1,A2,E,[Emin+1 50 60 80 100 140 180 220 300 380 500 600],'b','LineWidth',1.3)
contour(A1,A2,E,[-0.11 0.035 0.3 0.55 0.8 1.1],'b','LineWidth',1.3)
%contour(A1,A2,E,'LineWidth',1.3)

xlabel('x-axis','FontSize',16,'FontWeight','bold')
ylabel('y-axis','FontSize',16,'FontWeight','bold')
set(gca,'ytick',[-0.5 0 0.5])
set(gca,'xtick',[-0.4 0 0.4 0.8])
set(gca,'FontSize',16,'FontWeight','bold')

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














