C###  DAEPACK v0.3 - Copyright (C) M.I.T.
C###  - DERIVATIVE COMPUTATION -
      subroutine resad(neq,t,w,wdot,delta,ires,ichvar,rpar,ipar,
     $ zzzderiv,zzzne,zzzirn,zzzjcn,zzziw)
!!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###  Additional arguments for partial derivative computation
C###     zzzderiv - partial derivative array stored in sparse matrix
C###                format, pattern contained in zzzirn and zzzjcn.
C###     zzzne    - number of entries in zzzderiv.
C###     zzzirn   - row numbers for entries in zzzderiv.
C###     zzzjcn   - column numbers for entries in zzzderiv.
C###     zzziw    - integer workspace array.
      double precision zzzderiv(*)
      integer zzzne,zzzirn(zzzne),zzzjcn(zzzne),zzziw(*)

C###  Define elemenary variables and adjoints
      double precision zzzv1,zzzvbar1
      double precision zzzv2,zzzvbar2
      double precision zzzv3,zzzvbar3
      double precision zzzv4,zzzvbar4
      double precision zzzv5,zzzvbar5
      double precision zzzv6,zzzvbar6
      double precision zzzv7,zzzvbar7
      double precision zzzv8,zzzvbar8
      double precision zzzv9,zzzvbar9
      double precision zzzv10,zzzvbar10
      double precision zzzv11,zzzvbar11
      double precision zzzv12,zzzvbar12
      double precision zzzv13,zzzvbar13
      double precision zzzv14,zzzvbar14
      double precision zzzv15,zzzvbar15
      double precision zzzv16,zzzvbar16
      double precision zzzv17,zzzvbar17
C###  
C###  Active variable offsets and loop control variables
      integer deltaoft
      integer woft
      integer wdotoft
      integer x1dotoft
      integer x2dotoft
      integer x1oft
      integer x2oft
      integer x3oft
      integer v1oft
      integer v2oft
      integer v3oft
      integer z3oft
      integer z4oft
      integer zzzn
      integer zzzm
      integer zzzi
      integer zzzindx
C###  Variable offsets for independent, dependent, and 
C###  non-common block local active variables.
      deltaoft=0
      woft=deltaoft+neq
      wdotoft=woft+neq
      x1dotoft=wdotoft+neq
      x2dotoft=x1dotoft+1
      x1oft=x2dotoft+1
      x2oft=x1oft+1
      x3oft=x2oft+1
      v1oft=x3oft+1
      v2oft=v1oft+1
      v3oft=v2oft+1
      z3oft=v3oft+1
      z4oft=z3oft+1
C###  Number of independent and dependent variables.
      zzzn=neq+neq
      zzzm=neq
C###  Create SVM for functions/subroutines.
      call CREATESV(1,zzzn)
C###  Construct mapping between independent variable offsets and indices
      zzzindx=0
      do zzzi=1,neq
         zzzindx=zzzindx+1
         zzziw(zzzindx)=woft+zzzi
      end do
      do zzzi=1,neq
         zzzindx=zzzindx+1
         zzziw(zzzindx)=wdotoft+zzzi
      end do
C###  Set independent variable numbers
      call DSETIVCV(1,zzzn,zzziw)
C###  
C###  Compute elementary variables
      zzzv1=wdot(1)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      x1dot=zzzv1
C###  
      call DSVM1(x1dotoft+1,1,
     $           zzzvbar1,wdotoft+1,1)
C###  
C###  Compute elementary variables
      zzzv1=wdot(2)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      x2dot=zzzv1
C###  
      call DSVM1(x2dotoft+1,1,
     $           zzzvbar1,wdotoft+2,1)
C###  
C###  Compute elementary variables
      zzzv1=w(1)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      x1=zzzv1
C###  
      call DSVM1(x1oft+1,1,
     $           zzzvbar1,woft+1,1)
C###  
C###  Compute elementary variables
      zzzv1=w(2)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      x2=zzzv1
C###  
      call DSVM1(x2oft+1,1,
     $           zzzvbar1,woft+2,1)
C###  
C###  Compute elementary variables
      zzzv1=w(3)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      x3=zzzv1
C###  
      call DSVM1(x3oft+1,1,
     $           zzzvbar1,woft+3,1)
C###  
C###  Compute elementary variables
      zzzv1=w(4)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      v1=zzzv1
C###  
      call DSVM1(v1oft+1,1,
     $           zzzvbar1,woft+4,1)
C###  
C###  Compute elementary variables
      zzzv1=w(5)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      v2=zzzv1
C###  
      call DSVM1(v2oft+1,1,
     $           zzzvbar1,woft+5,1)
C###  
C###  Compute elementary variables
      zzzv1=w(6)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      v3=zzzv1
C###  
      call DSVM1(v3oft+1,1,
     $           zzzvbar1,woft+6,1)
C###  
C###  Compute elementary variables
      zzzv1=w(7)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      z3=zzzv1
C###  
      call DSVM1(z3oft+1,1,
     $           zzzvbar1,woft+7,1)
C###  
C###  Compute elementary variables
      zzzv1=w(8)
C###  Compute elementary partial derivatives
      zzzvbar1=1.0d0
C###  
      z4=zzzv1
C###  
      call DSVM1(z4oft+1,1,
     $           zzzvbar1,woft+8,1)
C###  
C###  Compute elementary variables
      zzzv1=x1
      zzzv2=x2
      zzzv3=zzzv1+zzzv2
      zzzv4=x3
      zzzv5=zzzv3-zzzv4
C###  Compute elementary partial derivatives
      zzzvbar5=1.0d0
      zzzvbar4=-zzzvbar5
      zzzvbar3=zzzvbar5
      zzzvbar2=zzzvbar3
      zzzvbar1=zzzvbar3
C###  
      delta(1)=zzzv5
C###  
      call DSVM3(deltaoft+1,1,
     $           zzzvbar4,x3oft+1,1,
     $           zzzvbar2,x2oft+1,1,
     $           zzzvbar1,x1oft+1,1)
C###  
C###  Compute elementary variables
      zzzv1=r3
      zzzv2=x1
      zzzv3=x2
      zzzv4=zzzv2+zzzv3
      zzzv5=zzzv1*zzzv4
      zzzv6=l3
      zzzv7=x1dot
      zzzv8=x2dot
      zzzv9=zzzv7+zzzv8
      zzzv10=zzzv6*zzzv9
      zzzv11=zzzv5+zzzv10
      zzzv12=v3
      zzzv13=zzzv11-zzzv12
C###  Compute elementary partial derivatives
      zzzvbar13=1.0d0
      zzzvbar12=-zzzvbar13
      zzzvbar11=zzzvbar13
      zzzvbar10=zzzvbar11
      zzzvbar9=zzzvbar10*zzzv6
      zzzvbar8=zzzvbar9
      zzzvbar7=zzzvbar9
      zzzvbar5=zzzvbar11
      zzzvbar4=zzzvbar5*zzzv1
      zzzvbar3=zzzvbar4
      zzzvbar2=zzzvbar4
C###  
      delta(2)=zzzv13
C###  
      call DSVM5(deltaoft+2,1,
     $           zzzvbar12,v3oft+1,1,
     $           zzzvbar8,x2dotoft+1,1,
     $           zzzvbar7,x1dotoft+1,1,
     $           zzzvbar3,x2oft+1,1,
     $           zzzvbar2,x1oft+1,1)
C###  
C###  Compute elementary variables
      zzzv1=100.0d0
      zzzv2=100.0d0
      zzzv3=1.0d0
      zzzv4=atan(zzzv3)
      zzzv5=zzzv2*zzzv4
      zzzv6=4.0d0
      zzzv7=zzzv5*zzzv6
      zzzv8=t
      zzzv9=zzzv7*zzzv8
      zzzv10=cos(zzzv9)
      zzzv11=zzzv1*zzzv10
      zzzv12=v1
      zzzv13=zzzv11-zzzv12
C###  Compute elementary partial derivatives
      zzzvbar13=1.0d0
      zzzvbar12=-zzzvbar13
C###  
      delta(3)=zzzv13
C###  
      call DSVM1(deltaoft+3,1,
     $           zzzvbar12,v1oft+1,1)
C###  
C###  Compute elementary variables
      zzzv1=v1
      zzzv2=-zzzv1
      zzzv3=v2
      zzzv4=zzzv2-zzzv3
C###  Compute elementary partial derivatives
      zzzvbar4=1.0d0
      zzzvbar3=-zzzvbar4
      zzzvbar2=zzzvbar4
      zzzvbar1=-zzzvbar2
C###  
      delta(4)=zzzv4
C###  
      call DSVM2(deltaoft+4,1,
     $           zzzvbar3,v2oft+1,1,
     $           zzzvbar1,v1oft+1,1)
C###  
C###  Compute elementary variables
      zzzv1=v1
      zzzv2=v3
      zzzv3=zzzv1-zzzv2
      zzzv4=z3
      zzzv5=zzzv3-zzzv4
C###  Compute elementary partial derivatives
      zzzvbar5=1.0d0
      zzzvbar4=-zzzvbar5
      zzzvbar3=zzzvbar5
      zzzvbar2=-zzzvbar3
      zzzvbar1=zzzvbar3
C###  
      delta(5)=zzzv5
C###  
      call DSVM3(deltaoft+5,1,
     $           zzzvbar4,z3oft+1,1,
     $           zzzvbar2,v3oft+1,1,
     $           zzzvbar1,v1oft+1,1)
C###  
C###  Compute elementary variables
      zzzv1=v2
      zzzv2=v3
      zzzv3=zzzv1-zzzv2
      zzzv4=z4
      zzzv5=zzzv3-zzzv4
C###  Compute elementary partial derivatives
      zzzvbar5=1.0d0
      zzzvbar4=-zzzvbar5
      zzzvbar3=zzzvbar5
      zzzvbar2=-zzzvbar3
      zzzvbar1=zzzvbar3
C###  
      delta(6)=zzzv5
C###  
      call DSVM3(deltaoft+6,1,
     $           zzzvbar4,z4oft+1,1,
     $           zzzvbar2,v3oft+1,1,
     $           zzzvbar1,v2oft+1,1)
C###  
      if((x1.gt.0.0d0.or.v1.gt.v3).and.(x2.le.0.0d0.and.v2.lt.v3))then
C###  Compute elementary variables
         zzzv1=v1
         zzzv2=a1
         zzzv3=x1
         zzzv4=zzzv2*zzzv3
         zzzv5=zzzv1-zzzv4
         zzzv6=a2
         zzzv7=zzzv5/zzzv6
         zzzv8=x1dot
         zzzv9=zzzv7-zzzv8
C###  Compute elementary partial derivatives
         zzzvbar9=1.0d0
         zzzvbar8=-zzzvbar9
         zzzvbar7=zzzvbar9
         zzzvbar5=zzzvbar7/zzzv6
         zzzvbar4=-zzzvbar5
         zzzvbar3=zzzvbar4*zzzv2
         zzzvbar1=zzzvbar5
C###  
         delta(7)=zzzv9
C###  
         call DSVM3(deltaoft+7,1,
     $              zzzvbar8,x1dotoft+1,1,
     $              zzzvbar3,x1oft+1,1,
     $              zzzvbar1,v1oft+1,1)
C###  
C###  Compute elementary variables
         zzzv1=0.0d0
         zzzv2=x2dot
         zzzv3=zzzv1-zzzv2
C###  Compute elementary partial derivatives
         zzzvbar3=1.0d0
         zzzvbar2=-zzzvbar3
C###  
         delta(8)=zzzv3
C###  
         call DSVM1(deltaoft+8,1,
     $              zzzvbar2,x2dotoft+1,1)
C###  
      else
         if((x2.gt.0.0d0.or.v2.gt.v3).and.(x1.lt.0.0d0.and.v1.lt.
     $      v3))then
C###  Compute elementary variables
            zzzv1=0.0d0
            zzzv2=x1dot
            zzzv3=zzzv1-zzzv2
C###  Compute elementary partial derivatives
            zzzvbar3=1.0d0
            zzzvbar2=-zzzvbar3
C###  
            delta(7)=zzzv3
C###  
            call DSVM1(deltaoft+7,1,
     $                 zzzvbar2,x1dotoft+1,1)
C###  
C###  Compute elementary variables
            zzzv1=v2
            zzzv2=a3
            zzzv3=x2
            zzzv4=zzzv2*zzzv3
            zzzv5=zzzv1-zzzv4
            zzzv6=a4
            zzzv7=zzzv5/zzzv6
            zzzv8=x2dot
            zzzv9=zzzv7-zzzv8
C###  Compute elementary partial derivatives
            zzzvbar9=1.0d0
            zzzvbar8=-zzzvbar9
            zzzvbar7=zzzvbar9
            zzzvbar5=zzzvbar7/zzzv6
            zzzvbar4=-zzzvbar5
            zzzvbar3=zzzvbar4*zzzv2
            zzzvbar1=zzzvbar5
C###  
            delta(8)=zzzv9
C###  
            call DSVM3(deltaoft+8,1,
     $                 zzzvbar8,x2dotoft+1,1,
     $                 zzzvbar3,x2oft+1,1,
     $                 zzzvbar1,v2oft+1,1)
C###  
         else
C###  Compute elementary variables
            zzzv1=a5
            zzzv2=v1
            zzzv3=zzzv1*zzzv2
            zzzv4=a6
            zzzv5=v2
            zzzv6=zzzv4*zzzv5
            zzzv7=zzzv3+zzzv6
            zzzv8=a7
            zzzv9=x1
            zzzv10=zzzv8*zzzv9
            zzzv11=zzzv7+zzzv10
            zzzv12=a8
            zzzv13=x2
            zzzv14=zzzv12*zzzv13
            zzzv15=zzzv11+zzzv14
            zzzv16=x1dot
            zzzv17=zzzv15-zzzv16
C###  Compute elementary partial derivatives
            zzzvbar17=1.0d0
            zzzvbar16=-zzzvbar17
            zzzvbar15=zzzvbar17
            zzzvbar14=zzzvbar15
            zzzvbar13=zzzvbar14*zzzv12
            zzzvbar11=zzzvbar15
            zzzvbar10=zzzvbar11
            zzzvbar9=zzzvbar10*zzzv8
            zzzvbar7=zzzvbar11
            zzzvbar6=zzzvbar7
            zzzvbar5=zzzvbar6*zzzv4
            zzzvbar3=zzzvbar7
            zzzvbar2=zzzvbar3*zzzv1
C###  
            delta(7)=zzzv17
C###  
            call DSVM5(deltaoft+7,1,
     $                 zzzvbar16,x1dotoft+1,1,
     $                 zzzvbar13,x2oft+1,1,
     $                 zzzvbar9,x1oft+1,1,
     $                 zzzvbar5,v2oft+1,1,
     $                 zzzvbar2,v1oft+1,1)
C###  
C###  Compute elementary variables
            zzzv1=a9
            zzzv2=v1
            zzzv3=zzzv1*zzzv2
            zzzv4=a10
            zzzv5=v2
            zzzv6=zzzv4*zzzv5
            zzzv7=zzzv3+zzzv6
            zzzv8=a11
            zzzv9=x1
            zzzv10=zzzv8*zzzv9
            zzzv11=zzzv7+zzzv10
            zzzv12=a12
            zzzv13=x2
            zzzv14=zzzv12*zzzv13
            zzzv15=zzzv11+zzzv14
            zzzv16=x2dot
            zzzv17=zzzv15-zzzv16
C###  Compute elementary partial derivatives
            zzzvbar17=1.0d0
            zzzvbar16=-zzzvbar17
            zzzvbar15=zzzvbar17
            zzzvbar14=zzzvbar15
            zzzvbar13=zzzvbar14*zzzv12
            zzzvbar11=zzzvbar15
            zzzvbar10=zzzvbar11
            zzzvbar9=zzzvbar10*zzzv8
            zzzvbar7=zzzvbar11
            zzzvbar6=zzzvbar7
            zzzvbar5=zzzvbar6*zzzv4
            zzzvbar3=zzzvbar7
            zzzvbar2=zzzvbar3*zzzv1
C###  
            delta(8)=zzzv17
C###  
            call DSVM5(deltaoft+8,1,
     $                 zzzvbar16,x2dotoft+1,1,
     $                 zzzvbar13,x2oft+1,1,
     $                 zzzvbar9,x1oft+1,1,
     $                 zzzvbar5,v2oft+1,1,
     $                 zzzvbar2,v1oft+1,1)
C###  
         end if
      end if

C###  Construct derivative matrix and sparsity pattern before returning

C###  Construct mapping between dependent variable offsets and indices
      zzzindx=0
      do zzzi=1,neq
         zzzindx=zzzindx+1
         zzziw(zzzindx)=deltaoft+zzzi
      end do
      call DCPRVS(1,zzzm,zzziw,zzzderiv,zzzne,zzzirn,zzzjcn)
C###  Delete SVM for functions/subroutines.
      call DELETESV(1)
C###  Done
      return
      end