"""Fresh rational verification of the paper's nominal reactor ellipsoid."""
import sympy as s
R=s.Rational
P=s.Matrix([[R(95,4),R(4075,142),R(3025,142)],[R(4075,142),R(57115,1207),R(1735,71)],[R(3025,142),R(1735,71),R(1920,71)]])

class InvariantEllipsoid:
    def __init__(self,model,center,level=R(1,5)):
        self.model=model;self.center=s.Matrix(center);self.level=R(level)
    def verify(self):
        if len(self.model.variables)!=3 or self.level<=0:raise ValueError('Three-dimensional model and positive level required.')
        x=self.model.variables;sub=dict(zip(x,self.center));f=s.Matrix(self.model.field);J=self.model.J.subs(sub)
        if any(s.simplify(v)!=0 for v in f.subs(sub)) or J.T*P+P*J!=-s.eye(3):return dict(status='not_applicable_to_this_operating_point')
        u=s.Matrix(s.symbols('ua ub uc'));shift=f.subs({v:c+w for v,c,w in zip(x,self.center,u)},simultaneous=True);remainder=s.simplify(shift-J*u);cubic=s.expand(2*(P*u).dot(remainder))
        # Recover the quadratic M from cubic = ua/10 * u^T M u.
        quotient=s.cancel(10*cubic/u[0]);M=s.hessian(quotient,list(u))/2
        if s.expand(quotient-(u.T*M*u)[0])!=0:return dict(status='unsupported_nonlinear_remainder')
        minors=[P[:i,:i].det() for i in [1,2,3]];frob=sum(v*v for v in M);bounds=[R(1,4),R(127,1000),R(161,1000)];Pinv=P.inv()
        passed=min(minors)>0 and frob<400 and all(Pinv[i,i]*self.level<bounds[i]**2 for i in range(3)) and all(self.center[i]>bounds[i] for i in range(3))
        if not passed:return dict(status='not_certified_at_requested_level',frobenius_squared=str(frob))
        return dict(status='certified_forward_invariant',P=P.tolist(),principal_minors=minors,lyapunov_identity=True,cubic_matrix=M.tolist(),frobenius_squared=frob,coordinate_radius_upper=bounds,level=self.level,multiplier_b_floor=self.center[1]-bounds[1],V_exponential_time_constant=2*s.trace(P),euclidean_ball_radius=R(45,1000),ball_contained=bool(s.trace(P)*R(45,1000)**2<=self.level),scope='Fixed rate constants and load, with initial state inside this ellipsoid; no entry or time-varying-load guarantee.')
    def contains(self,state):
        u=s.Matrix(state)-self.center;return bool((u.T*P*u)[0]<=self.level)