P<X> := PolynomialRing(Rationals());
K<phi> := NumberField(X^2-X-1);
OK := RingOfIntegers(K);

E := EllipticCurve([ 0, phi - 1, phi + 1, -2*phi, 0 ]);
E1 := EllipticCurve([ 1, -phi - 1, phi + 1, 
                    3904*phi - 6577, 53806*phi - 88477 ]);
E2 := EllipticCurve([ phi, phi, phi + 1, 
          -5179442*phi - 3201071, -6788348449*phi - 4195430069 ]);

assert Conductor(E) eq ideal<OK|3*phi-21>;
assert Conductor(E1) eq ideal<OK|-6*phi+42>;
assert Conductor(E2) eq ideal<OK|-6*phi+42>;

  c4,c6 := Explode(cInvariants(E));
  a := -27*c4;
  b := -54*c6;

RR<x,y,z> := PolynomialRing(K,3);

// Formula for X_E(7) due to Halberstadt and Kraus
F := a*x^4 + 7*b*x^3*z + 3*x^2*y^2 - 3*a^2*x^2*z^2 - 6*b*x*y*z^2 
       - 5*a*b*x*z^3 + 2*y^3*z + 3*a*y^2*z^2 + 2*a^2*y*z^3 - 4*b^2*z^4;

// Some minimising and (ad hoc) reducing suggested the following 
// change of coordinates.

T := Matrix(K,3,3,[ -36*phi + 12, 48*phi + 24, -84*phi + 24, 
  1224*phi - 1440, 720*phi - 1008, 2880*phi - 4752, 
		       -2*phi, 1, phi - 1 ]);
mu := 14693280768*phi + 14693280768;
F1 := (1/mu)*F^T;

C := Curve(Proj(RR),F);
C1 := Curve(Proj(RR),F1);

// I did a PointSearch on the restriction of scalars, to find 
// these (rather small) points
pts1 := [[phi - 1,phi,1 ], [4*phi - 6, phi, 1 ], [ phi - 1, -phi + 3, 1 ]];
pts1 := [C1!pt : pt in pts1];

// Back to the original curve
pts := [ [ 48*phi + 24, 720*phi + 936, 1 ], 
	 [ 12*phi + 24, 720*phi - 576, 1 ], [ 0, 1, 0 ]];
pts := [C!pt : pt in pts];

assert pts eq [C!Eltseq(Vector(Eltseq(pt))*Transpose(T)): pt in pts1];

Deriv := Derivative;

function Invariant(F)
  R9<x1,y1,z1,x2,y2,z2,x3,y3,z3> := PolynomialRing(K,9);
  poly := Evaluate(F,[x1,y1,z1])*Evaluate(F,[x2,y2,z2])*Evaluate(F,[x3,y3,z3]);
  G := SymmetricGroup(3);
  for i in [1..4] do
    poly := &+[Sign(g)*Deriv(Deriv(Deriv(poly,1^g),3+(2^g)),6+(3^g)): g in G];
  end for;
  return (1/5184)*K!poly;
end function;

function Covariants(F)
  Deriv := Derivative;
  Hmat := Matrix(3,3,[Deriv(Deriv(F,i),j): i,j in [1..3]]);
  Psi_6 := (-1/54)*Determinant(Hmat);
  vv := Matrix(3,1,[Deriv(Psi_6,i): i in [1..3]]);
  HHmat := HorizontalJoin(Hmat,vv);
  HHmat := VerticalJoin(HHmat,Matrix(1,4,Eltseq(vv) cat [0]));
  Psi_14 := (1/9)*Determinant(HHmat);
  Jmat := Matrix(3,3,[Deriv(f,i): i in [1..3],f in [F,Psi_6,Psi_14]]);
  Psi_21 := (1/14)*Determinant(Jmat);
  return Psi_6,Psi_14,Psi_21;
end function;

cubics := [
      2*(a*x^2 + 3*b*x*z + 3*y^2 + 2*a*y*z)*z,
      2*(-b*x^3 - 2*a*x^2*y - 6*b*x*y*z - a*b*x*z^2 - 8*y^3 
         - 8*a*y^2*z - 2*a^2*y*z^2 + (4*a^3 + 28*b^2)*z^3),
      2*(-b*x^3 - a^2*x^2*z - 12*b*x*y*z - 8*a*b*x*z^2 - 2*y^3 
         + a*y^2*z + 4*a^2*y*z^2 + (4*a^3 + 16*b^2)*z^3),
      2*(3*a*b*x^3 - 2*a^2*x^2*y + (a^3 + 21*b^2)*x^2*z - 2*a*b*x*y*z 
	  - 2*a^2*b*x*z^2 - a^2*y^2*z + 6*b^2*y*z^2 + a*b^2*z^3) 
];

curvelist := {};

for pt in pts do
  assert Evaluate(F,Eltseq(pt)) eq 0;
  Psi_0 := Invariant(F);
  Psi_6,Psi_14,Psi_21 := Covariants(F);
  cc4 := Evaluate(Psi_14,Eltseq(pt));
  cc6 := Evaluate(Psi_21,Eltseq(pt));
  dd := [Evaluate(d,Eltseq(pt)): d in cubics];
  assert exists(d){d : d in dd | d ne 0};
  assert Evaluate(Psi_6,Eltseq(pt)) ne 0; // i.e. not a cusp
  E0 := MinimalModel(EllipticCurve([-27*cc4*d^2,54*cc6*d^3]));
  Include(~curvelist,E0);
end for;

assert curvelist eq {E,E1,E2};