H0=-y*D[f,y]+x0*D[f,x0];
H1=-3/2*f-y*D[f,y]-2*x0*D[f,x0];
b0=y*D[f,x0];
g0=x0*I*f;
bm0=-x0*D[f,y];
gm0=-y*D[D[I*f,x0],y]-x0*D[D[I*f,x0],x0]-3/2*D[I*f,x0];
mom=-y*D[D[I*f,y],y]-3/2*D[I*f,y]-x0*D[D[I*f,x0],y];