/* damped, linear oscillation: equation theta'' + 2kt theta' + w^2 theta =0
   2k: damping factor
   w: pulse when k=0

   if k>w, decreasing exponential behaviour
   if k=w: limit
   if k<w: damped oscillation, with pulse sqrt(w^2-k^2)

   we take a sampling period of 1 sec
*/

var thetap,thetav;
initial linit;

linit -> l0 :
      guard thetap=8 and thetav=0 assign { };
linit -> l1 :
      guard thetap=8 and thetav=0 assign { };
linit -> l2 :
      guard thetap=8 and thetav=0 assign { };
linit -> l3 :
      guard thetap=8 and thetav=0 assign { };
linit -> l4 :
      guard thetap=8 and thetav=0 assign { };
linit -> l5 :
      guard thetap=8 and thetav=0 assign { };
linit -> l6 :
      guard thetap=8 and thetav=0 assign { };

l0 -> l0:
      /* we take w=1/8 and k=0
	 load("diag"); set_display('xml)$
	 m:matrix([1,1,0],[-1/64,1,0],[0,0,1]); j:jordan(m); mj:dispJordan(j); pinv:ModeMatrix(m,j);p:invert(pinv);pinv.mj.p;  set_display('none)$ float(rectform(pinv)); float(rectform(j)); float(rectform(p));
      */
      guard true assign { thetap := thetap + thetav; thetav := - 1/64 thetap + thetav; }
      jordan
      matrix([1.0,1.0,0.0],[-0.125*%i,0.125*%i,0.0],[0.0,0.0,1.0])
      [[1.0-0.125*%i,1.0],[0.125*%i+1.0,1.0],[1.0,1.0]]
      matrix([0.5,4.0*%i,0.0],[0.5,-4.0*%i,0.0],[0.0,0.0,1.0])
      ;

l1 -> l1:
      guard true
      /* we take w=1/8 and k=1/64
	 load("diag"); set_display('xml)$
	 m:matrix([1,1,0],[-1/64,62/64,0],[0,0,1]); j:jordan(m); mj:dispJordan(j); pinv:ModeMatrix(m,j);p:invert(pinv);pinv.mj.p;  set_display('none)$ float(rectform(pinv)); float(rectform(j)); float(rectform(p));
      */
      assign { thetap := thetap + thetav; thetav := - 1/64 thetap + 62/64 thetav; }
      jordan
      matrix([1.0,1.0,0.0],
	       [-.1240195927061527*%i-0.015625,.1240195927061527*%i-0.015625,
		0.0],[0.0,0.0,1.0])
      [[0.984375-.1240195927061527*%i,1.0],
	 [.1240195927061527*%i+0.984375,1.0],[1.0,1.0]]
      matrix([0.0629940788348712*%i+0.5,4.031621045431756*%i,0.0],
	       [0.5-0.0629940788348712*%i,-4.031621045431756*%i,0.0],
	       [0.0,0.0,1.0])
;

l2 -> l2:
      guard true
      /* we take w=1/8 and k=1/16
	 set_display('xml)$
	 m:matrix([1,1,0],[-1/64,14/16,0],[0,0,1]); j:jordan(m); mj:dispJordan(j); pinv:ModeMatrix(m,j);p:invert(pinv);pinv.mj.p;  set_display('none)$ float(rectform(pinv)); float(rectform(j)); float(rectform(p));
      */
      assign { thetap := thetap + thetav; thetav := - 1/64 thetap + 14/16 thetav; }
      jordan
      matrix([1.0,1.0,0.0],
	       [-.1082531754730548*%i-0.0625,.1082531754730548*%i-0.0625,0.0],
	       [0.0,0.0,1.0])
      [[0.9375-.1082531754730548*%i,1.0],[.1082531754730548*%i+0.9375,1.0],
	 [1.0,1.0]]
      matrix([.2886751345948129*%i+0.5,4.618802153517007*%i,0.0],
	       [0.5-.2886751345948129*%i,-4.618802153517007*%i,0.0],
	       [0.0,0.0,1.0])
;
l3 -> l3:
      guard true
      /* we take w=1/8 and k=1/8
	 set_display('xml)$
	 m:matrix([1,1,0],[-1/64,6/8,0],[0,0,1]); j:jordan(m); mj:dispJordan(j); pinv:ModeMatrix(m,j);p:invert(pinv);pinv.mj.p;  set_display('none)$ float(rectform(pinv)); float(rectform(j)); float(rectform(p));
      */
      assign { thetap := thetap + thetav; thetav := - 1/64 thetap + 6/8 thetav; }
      jordan
      matrix([0.125,1.0,0.0],[-0.015625,0.0,0.0],[0.0,0.0,1.0])
      [[0.875,2.0],[1.0,1.0]]
      matrix([0.0,-64.0,0.0],[1.0,8.0,0.0],[0.0,0.0,1.0])
;
l4 -> l4:
      /* we take w=1/8 and k=2/8
	 set_display('xml)$
	 m:matrix([1,1,0],[-1/64,4/8,0],[0,0,1]); j:jordan(m); mj:dispJordan(j); pinv:ModeMatrix(m,j);p:invert(pinv);pinv.mj.p;  set_display('none)$ float(rectform(pinv)); float(rectform(j)); float(rectform(p));
      */
      guard true assign { thetap := thetap + thetav; thetav := - 1/64 thetap + 4/8 thetav; }
      jordan
      matrix([1.0,1.0,0.0],[-.4665063509461096,-.03349364905389035,0.0],
	       [0.0,0.0,1.0])
      [[.5334936490538904,1.0],[.9665063509461096,1.0],[1.0,1.0]]
      matrix([-0.0773502691896258,-2.309401076758503,0.0],
	       [1.077350269189626,2.309401076758503,0.0],[0.0,0.0,1.0])
      ;

l5 -> l5:
      /* we take w=1/8 and k=1/2
	 set_display('xml)$
	 m:matrix([1,1,0],[-1/64,0,0],[0,0,1]); j:jordan(m); mj:dispJordan(j); pinv:ModeMatrix(m,j);p:invert(pinv);pinv.mj.p;  set_display('none)$ float(rectform(pinv)); float(rectform(j)); float(rectform(p));
      */
      guard true assign { thetap := thetap + thetav; thetav := - 1/64 thetap; }
      jordan
 matrix([1.0,1.0,0.0],[-.9841229182759271,-.01587708172407287,0.0],
	       [0.0,0.0,1.0])
 [[.01587708172407287,1.0],[.9841229182759271,1.0],[1.0,1.0]]
 matrix([-.01639777949432223,-1.032795558988644,0.0],
	       [1.016397779494322,1.032795558988644,0.0],[0.0,0.0,1.0])
;

l6 -> l6:
      /* we take w=1/8 and k=3/4
	 set_display('xml)$
	 m:matrix([1,1,0],[-1/64,-1/2,0],[0,0,1]); j:jordan(m); mj:dispJordan(j); pinv:ModeMatrix(m,j);p:invert(pinv);pinv.mj.p;  set_display('none)$ float(rectform(pinv)); float(rectform(j)); float(rectform(p));
      */
      guard true assign { thetap := thetap + thetav; thetav := - 1/64 thetap - 1/2 thetav; }
      jordan
matrix([1.0,1.0,0.0],[-1.489509972887452,-.01049002711254798,0.0],
	       [0.0,0.0,1.0])
[[-0.489509972887452,1.0],[0.989509972887452,1.0],[1.0,1.0]]
matrix([-0.00709255283710994,-.6761234037828133,0.0],
	       [1.00709255283711,.6761234037828133,0.0],[0.0,0.0,1.0]);