///////////////////////////////////
// Projections to xz-coordinates //
///////////////////////////////////

LIB "all.lib";
LIB "poly.lib";
LIB "tropical.lib";
LIB "elim.lib";

ring rr = (0,b2,b3,b4,b5), (x,y,z),dp;
poly f=y^2-x*(x-b2^2)*(x-b3^2)*(x-b4^2)*(x+b5^2);
// f;
// -x^5+(b2^2+b3^2+b4^2-b5^2)*x^4+(-b2^2*b3^2-b2^2*b4^2+b2^2*b5^2-b3^2*b4^2+b3^2*b5^2+b4^2*b5^2)*x^3+(b2^2*b3^2*b4^2-b2^2*b3^2*b5^2-b2^2*b4^2*b5^2-b3^2*b4^2*b5^2)*x^2+y^2+(b2^2*b3^2*b4^2*b5^2)*x

ideal I = (f,z-y+b3*b4*b5*x-b5*x^2);

poly fnoy= eliminate(I,y)[1];


// Replace z by y so that the pictures are not flipped.
ring r2 = (0,b2, b3, b4, b5),(x,y),dp;
setring(r2);
map PP = rr,x,0,y;
poly newfnoy = PP(fnoy);
// newfnoy;
// x^5+(-b2^2-b3^2-b4^2)*x^4+(b2^2*b3^2+b2^2*b4^2-b2^2*b5^2+b3^2*b4^2-b3^2*b5^2+2*b3*b4*b5^2-b4^2*b5^2)*x^3+(2*b5)*x^2*y+(-b2^2*b3^2*b4^2+b2^2*b3^2*b5^2+b2^2*b4^2*b5^2)*x^2+(-2*b3*b4*b5)*x*y-y^2+(-b2^2*b3^2*b4^2*b5^2)*x

matrix m = coef(newfnoy,xy);
m;
m[1,1]=x^5
m[1,2]=x^4
m[1,3]=x^3
m[1,4]=x^2*y
m[1,5]=x^2
m[1,6]=x*y
m[1,7]=y^2
m[1,8]=x
m[2,1]=1
m[2,2]=(-b2^2-b3^2-b4^2)
m[2,3]=(b2^2*b3^2+b2^2*b4^2-b2^2*b5^2+b3^2*b4^2-b3^2*b5^2+2*b3*b4*b5^2-b4^2*b5^2)
m[2,4]=(2*b5)
m[2,5]=(-b2^2*b3^2*b4^2+b2^2*b3^2*b5^2+b2^2*b4^2*b5^2)
m[2,6]=(-2*b3*b4*b5)
m[2,7]=-1
m[2,8]=(-b2^2*b3^2*b4^2*b5^2)


poly mx4 = m[2,2];
poly mx3 = m[2,3];
poly mx2 = m[2,5];

// We write the output coefficients:
write(":w coefficientsProjectionMxz.txt", "MX4 =", mx4, "", "MX3 =", mx3, "", "MX2 =", mx2, "");