function feq=HillRadINC(yi,C,in)
global GEO wt mu
U=pot(0,yi,0,GEO,mu);

Vi=abs(cos(in)*sqrt(mu/yi)-wt*yi);
vx=Vi;
vz=sqrt(mu/yi)*sin(in);

feq=1/2*wt^2*yi^2+U-C-1/2*(vx^2+vz^2);
end