C =============================================================== C Model of Example 4 of Birta et. al(1985) C Ref : L.G.Birta, T.I.Oren, and D.L.Kettenis, "A Robust C Procedure for Discontinuity Handling in Continuous C System Simulation", Tran. Soc. Computer Simulation, C Vol.2, No.3, p.189, 1985 C =============================================================== subroutine res(neq,t,w,wdot,delta,ires,ichvar,rpar,ipar) !!independent{ w,wdot } !!dependent { delta } implicit none integer neq,ires,ichvar,ipar(1) double precision t,w(neq),wdot(neq),delta(neq),rpar(1) double precision x1,x2,x3,x1dot,x2dot,v1,v2,v3,z3,z4 double precision a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12 parameter(a1=12.0d0,a2=0.24d0,a3=12.0d0,a4=0.24d0) parameter(a5=13.64d0,a6=-11.64d0,a7=-50.0d0,a8=0.0d0) parameter(a9=-11.64d0,a10=13.64d0,a11=0.0d0,a12=-50.0d0) double precision r3,l3 parameter(r3=10.0d0,l3=0.2d0) double precision cos,atan intrinsic cos,atan c x1dot=wdot(1) x2dot=wdot(2) x1=w(1) x2=w(2) x3=w(3) v1=w(4) v2=w(5) v3=w(6) z3=w(7) z4=w(8) c delta(1)=x1+x2-x3 delta(2)=r3*(x1+x2)+l3*(x1dot+x2dot)-v3 delta(3)=100.0d0*cos(100.0d0*atan(1.0d0)*4.0d0*t)-v1 delta(4)=-v1-v2 delta(5)=v1-v3-z3 delta(6)=v2-v3-z4 if((x1.gt.0.0d0.or.v1.gt.v3).and. $ (x2.le.0.0d0.and.v2.lt.v3)) then delta(7)=(v1-a1*x1)/a2-x1dot delta(8)=0.0d0-x2dot else if((x2.gt.0.0d0.or.v2.gt.v3).and. $ (x1.lt.0.0d0.and.v1.lt.v3)) then delta(7)=0.0d0-x1dot delta(8)=(v2-a3*x2)/a4-x2dot else delta(7)=a5*v1+a6*v2+a7*x1+a8*x2-x1dot delta(8)=a9*v1+a10*v2+a11*x1+a12*x2-x2dot end if end if c return end