c-----------------------------------------------------------------------
c     these are the subroutines for the
c
c        Ring Modulator (ODE case) (cs = 2.0d-12)
c        ODE of dimension 15
c
c-----------------------------------------------------------------------
      subroutine prob(problm,neqn,tbegin,tend,ijac,mljac,mujac)
      character*(*) problm
      integer neqn,ijac,mljac,mujac
      double precision tbegin,tend

      problm = 'hardring'
      neqn   = 15
      tbegin = 0d0
      tend   = 1d-3
      ijac   = 1
      mljac  = neqn
      mujac  = neqn

      return
      end
c-----------------------------------------------------------------------
      subroutine init(neqn,y)
      integer neqn
      double precision y(neqn)

      integer i

      do 10 i=1,neqn
         y(i) = 0d0
   10 continue

      return
      end
c-----------------------------------------------------------------------
      subroutine feval(neqn,t,y,dy)
      integer neqn
      double precision t,y(neqn),dy(neqn)

      double precision c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc,
     +                 gamma,delta,pi,
     +                 uin1,uin2,ud1,ud2,ud3,ud4,qud1,qud2,qud3,qud4
      parameter (c=1.6d-8,cs=2.0d-12,cp=1d-8,r=25d3,rp=50d0,
     +           lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4,
     +           rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2,
     +           gamma=40.67286402d-9,delta=17.7493332d0,
     +           pi=3.141592653589793238462643383d0)

      uin1   = 0.5d0*dsin(2d3*pi*t)
      uin2   = 2d0*dsin(2d4*pi*t)
      ud1    = +y(3)-y(5)-y(7)-uin2
      ud2    = -y(4)+y(6)-y(7)-uin2
      ud3    = +y(4)+y(5)+y(7)+uin2
      ud4    = -y(3)-y(6)+y(7)+uin2
      qud1   = gamma*(dexp(delta*ud1)-1d0)
      qud2   = gamma*(dexp(delta*ud2)-1d0)
      qud3   = gamma*(dexp(delta*ud3)-1d0)
      qud4   = gamma*(dexp(delta*ud4)-1d0)

      dy(1)  = (y(8)-0.5d0*y(10)+0.5d0*y(11)+y(14)-y(1)/r)/c
      dy(2)  = (y(9)-0.5d0*y(12)+0.5d0*y(13)+y(15)-y(2)/r)/c
      dy(3)  = (y(10)-qud1+qud4)/cs
      dy(4)  = (-y(11)+qud2-qud3)/cs
      dy(5)  = (y(12)+qud1-qud3)/cs
      dy(6)  = (-y(13)-qud2+qud4)/cs
      dy(7)  = (-y(7)/rp+qud1+qud2-qud3-qud4)/cp
      dy(8)  = -y(1)/lh
      dy(9)  = -y(2)/lh
      dy(10) = (0.5d0*y(1)-y(3)-rg2*y(10))/ls2
      dy(11) = (-0.5d0*y(1)+y(4)-rg3*y(11))/ls3
      dy(12) = (0.5d0*y(2)-y(5)-rg2*y(12))/ls2
      dy(13) = (-0.5d0*y(2)+y(6)-rg3*y(13))/ls3
      dy(14) = (-y(1)+uin1-(ri+rg1)*y(14))/ls1
      dy(15) = (-y(2)-(rc+rg1)*y(15))/ls1

      return
      end
c-----------------------------------------------------------------------
      subroutine jeval(neqn,t,y,jac,ldim)
      integer neqn,ldim
      double precision t,y(neqn),jac(ldim,neqn)

      integer i,j
      double precision c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc,
     +                 gamma,delta,pi,
     +                 uin2,ud1,ud2,ud3,ud4,qpud1,qpud2,qpud3,qpud4
      parameter (c=1.6d-8,cs=2.0d-12,cp=1d-8,r=25d3,rp=50d0,
     +           lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4,
     +           rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2,
     +           gamma=40.67286402d-9,delta=17.7493332d0,
     +           pi=3.141592653589793238462643383d0)

      do 20 j=1,neqn
         do 10 i=1,neqn
            jac(i,j) = 0d0
   10    continue
   20 continue

      uin2  = 2d0*dsin(2d4*pi*t)
      ud1   = +y(3)-y(5)-y(7)-uin2
      ud2   = -y(4)+y(6)-y(7)-uin2
      ud3   = +y(4)+y(5)+y(7)+uin2
      ud4   = -y(3)-y(6)+y(7)+uin2
      qpud1 = gamma*delta*dexp(delta*ud1)
      qpud2 = gamma*delta*dexp(delta*ud2)
      qpud3 = gamma*delta*dexp(delta*ud3)
      qpud4 = gamma*delta*dexp(delta*ud4)

      jac(1,1)   = -1d0/(c*r)
      jac(1,8)   = 1d0/c
      jac(1,10)  = -0.5d0/c
      jac(1,11)  = -jac(1,10)
      jac(1,14)  = jac(1,8)
      jac(2,2)   = jac(1,1)
      jac(2,9)   = jac(1,8)
      jac(2,12)  = jac(1,10)
      jac(2,13)  = jac(1,11)
      jac(2,15)  = jac(1,14)
      jac(3,3)   = (-qpud1-qpud4)/cs
      jac(3,5)   = qpud1/cs
      jac(3,6)   = -qpud4/cs
      jac(3,7)   = (qpud1+qpud4)/cs
      jac(3,10)  = 1d0/cs
      jac(4,4)   = (-qpud2-qpud3)/cs
      jac(4,5)   = -qpud3/cs
      jac(4,6)   = qpud2/cs
      jac(4,7)   = (-qpud2-qpud3)/cs
      jac(4,11)  = -1d0/cs
      jac(5,3)   = qpud1/cs
      jac(5,4)   = -qpud3/cs
      jac(5,5)   = (-qpud1-qpud3)/cs
      jac(5,7)   = (-qpud1-qpud3)/cs
      jac(5,12)  = 1d0/cs
      jac(6,3)   = -qpud4/cs
      jac(6,4)   = qpud2/cs
      jac(6,6)   = (-qpud2-qpud4)/cs
      jac(6,7)   = (qpud2+qpud4)/cs
      jac(6,13)  = -1d0/cs
      jac(7,3)   = (qpud1+qpud4)/cp
      jac(7,4)   = (-qpud2-qpud3)/cp
      jac(7,5)   = (-qpud1-qpud3)/cp
      jac(7,6)   = (qpud2+qpud4)/cp
      jac(7,7)   = (-qpud1-qpud2-qpud3-qpud4-1d0/rp)/cp
      jac(8,1)   = -1d0/lh
      jac(9,2)   = jac(8,1)
      jac(10,1)  = 0.5d0/ls2
      jac(10,3)  = -1d0/ls2
      jac(10,10) = -rg2/ls2
      jac(11,1)  = -0.5d0/ls3
      jac(11,4)  = 1d0/ls3
      jac(11,11) = -rg3/ls3
      jac(12,2)  = jac(10,1)
      jac(12,5)  = jac(10,3)
      jac(12,12) = jac(10,10)
      jac(13,2)  = jac(11,1)
      jac(13,6)  = jac(11,4)
      jac(13,13) = jac(11,11)
      jac(14,1)  = -1d0/ls1
      jac(14,14) = -(ri+rg1)/ls1
      jac(15,2)  = jac(14,1)
      jac(15,15) = -(rc+rg1)/ls1

      return
      end
c-----------------------------------------------------------------------
      subroutine solut(neqn,yend)
      integer neqn
      double precision yend(neqn)

       yend(1 ) = -0.233905735842070174d-1 
       yend( 2) = -0.736748548598215210d-2 
       yend( 3) =  0.258295670891733276d+0     
       yend( 4) = -0.406446572164108511d+0     
       yend( 5) = -0.403945566551075719d+0     
       yend( 6) =  0.260796676505321678d+0     
       yend( 7) =  0.110676186126208484d+0     
       yend( 8) =  0.293990434238965346d-6 
       yend( 9) = -0.284002993248608379d-7 
       yend(10) =  0.726719826722729682d-3 
       yend(11) =  0.792948719688146885d-3 
       yend(12) = -0.725528349561929123d-3 
       yend(13) = -0.794140196849026266d-3 
       yend(14) =  0.708849541688199994d-4 
       yend(15) =  0.239005907525775679d-4 
      return
      end