####################################################################
#
#   The copyright to this file, and the images it draws are owned by 
# Frank Sottile.  We allow the free use of this, but you must first
# ask permission of either Frank Sottile or Adriana Ortiz, and 
# acknowledge our creation of these images.
#
####################################################################
interface(quiet=true):
with(plots):
plotsetup(x11,plotoptions=`width=6in,height=6in`):



f:=  -2*y^2  +2*x*y   +12*x^2
    +10*y^3  +3*x*y^2 -10*x^2*y  -13*x^3
    -11*y^4  +6*x*y^3  +9*x^2*y^2 -2*x^3*y -x^4:

#lprint(primpart(f));

H:=linalg[det](linalg[hessian](f,[x,y])):


#A:=implicitplot(H,x=-6.8..4.8,y=-2.8..3.5,grid=[400,400],axes=none):
#plotsetup(gif,plotoutput=sprintf("4.gif"),plotoptions=`height=400,width=400`):
#display(A);
#B:=implicitplot(0,x=-6.8..4.8,y=-2.8..3.5,grid=[200,200],axes=none,color=red):
A:=implicitplot(H,x=-6.8..4.8,y=-2.8..3.5,grid=[200,200],axes=none,color=blue,thickness=1):
#plotsetup(ps,plotoutput=`B.eps`,plotoptions=`color,portrait,width=2in,height=2in,noborder`):
#
#   This draws the annotated picture
#
P:=pointplot([[0,1/5],[-2,0],[2,2],[2,-1]], color=red, symbol=POINT):
l1:=plot([-2/5,t,t=-2.8..3.5],color=magenta, numpoints=2):
l2:=plot([t,3/4-t/2,t=-6.8..4.8],color=magenta, numpoints=2):
l3:=plot([t,-1/4+t/2,t=-6.8..4.8],color=magenta, numpoints=2):
xax:=plot([t,0,t=-7...4.8],color=black, numpoints=2):
yax:=plot([0,t,t=-2.8..3.5],color=black, numpoints=2):
t1:=plot([-4,t,t=-1/8..1/8],color=black,numpoints=2):
t2:=plot([3,t,t=-1/8..1/8],color=black,numpoints=2):
t3:=plot([t,-2,t=-1/8..1/8],color=black,numpoints=2):
t4:=plot([t,2,t=-1/8..1/8],color=black,numpoints=2):
a1:=plot([-7.+abs(t),t,t=-1/5..1/5],color=black,numpoints=3):
a2:=plot([4.8-abs(t),t,t=-1/5..1/5],color=black,numpoints=3):
a3:=plot([t,-2.8+abs(t),t=-1/5..1/5],color=black,numpoints=3):
a4:=plot([t,3.5-abs(t),t=-1/5..1/5],color=black,numpoints=3):
Axes:=[xax,yax,t1,t2,t3,t4,a1,a2,a3,a4]:
display(A,P,l1,l2,l3,Axes[],view=[-7.05..4.85,-2.85..3.55],scaling=CONSTRAINED);

#interface(quiet=false):
#
#  Verifies that there are no points at infinity
#
#topH:= subs(z=0, simplify(z^4*subs(x=x/z,y=y/z,H)));
#A:=simplify(-topH/36);
#B:=simplify(A-(x^2+3*x*y-13/3*y^2)^2 - (sqrt(6)*x^2 - sqrt(506/9)*y^2)^2);
#evalf(B);

subs(x=0,y=1/5,H); # 5124/125
subs(x=-2,y=0,H);   # 7068
subs(x=2,y=2,H);    # 8508
subs(x=2,y=-1,H);   # 6828


A:=-625*subs(y=y/5,x=-2/5,H)/4;

simplify( A - 35000*(1 -y)^2 - 481*(3 -4*y + y^2)^2 - (19 -21*y + 64/42*y^2)^2);