function xdot = whaneq(t,x) %matrix formulation of Whan's equation
%constraint Ib as part of the vector
global  betac Ib betaL Rs Rsg
%
%
%
%The Matric, x(1)= phase, x(2)= voltage , x(3)= shunt current, x(4)=current
%
JJM= [ 0,                1,                 0,                 0
 (-sin(x(1))./(x(1)*betac)), -betaover,   -betaover,     1/betaL
            0                1/betaL      0              1/betaL
            0                 0           0                    0];
xdot=JJM*x;

