anchor.diff:    alpha%gr = alpha%J*alpha%omega_dot
   
   end subroutine create_body
  +
  +
  +
  +subroutine disp(t,y)
  +
  +  use constants,only:dp,sp
  +  implicit none
  +
  +  real(dp)    :: t, y(12)
  +  integer(sp) :: i
  +
  +  write(*,'(13f13.2)') t, (y(i),i=1,12)
  +
  +
  +end subroutine disp
  +
  +
  +subroutine disp_mat(t, A,m,n)
  +
  +  use constants, only:dp,sp
  +  implicit none
  +
  +  real(dp), intent(in) :: t
  +  integer(sp)          :: i, j, m, n
  +  real(dp)             :: A(m,n)
  +
  +
  +  do i = 1, n
  +     write(*,'(13f13.2)') t, (A(j,i),j=1,m)
  +  end do
  +
  +
  +end subroutine disp_mat
  
positive.code:   use constants
    use utils,only:get_vector_elements
    use body_mod,only:alpha, R_rigid
  !  use jacobian_mod, only: res_rigid
    implicit none
  
    real(dp)  :: t
  
negative.code: 
    !?? get the residual
    delta = get_vector_elements(R_rigid(alpha),4)
  !  print*, "residual vector = ", delta
  
    call disp(t,delta)
  
  end subroutine res
  
  !********************************
  ! implementation of the jacobian
  !********************************
  subroutine jac(t,y,yprime,pd,cj,rpar,ipar)
  !!$C  to define the matrix of partial derivatives
  !!$C             PD=DG/DY+CJ*DG/DYPRIME
  !!$C         CJ is a scalar which is input to JAC.
  !!$C         For the given values of T,Y,YPRIME, the
  !!$C         subroutine must evaluate the non-zero partial
  !!$C         derivatives for each equation and each solution
  !!$C         component, and store these values in the
  !!$C         matrix PD.  The elements of PD are set to zero
  !!$C         before each call to JAC so only non-zero elements
  !!$C         need to be defined.
  !!$C
  !!$C         Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ.
  !!$C         You must declare the name JAC in an EXTERNAL statement in
  !!$C         your program that calls DDASSL.  You must dimension Y,
  !!$C         YPRIME and PD in JAC.
  
    use constants,only:dp,sp
    use body_mod,only:alpha
    use jacobian_mod,only:jac_rigid
    implicit none
    real(dp)     :: cj
    real(dp)     :: t
    real(dp)     :: y(12)
    real(dp)     :: yprime(12)
    integer(sp)  :: ires
    real(dp)     :: rpar(*)
    integer(sp)  :: ipar(*)
    real(dp)     :: pd(12,12)
    integer(sp)  :: i,j,k
  
    ! JACOBIAN  = jac_rigid(alpha,a) !?? implement jacobian
   ! print*,cj 
  !  stop
    cj = 0.1
    pd = jac_rigid (alpha, cj)
  
    print*,  shape(pd)
    
    call disp_mat(t, pd, 12, 12)
  
    stop
  
    !  stop"where is the impl"
  
  end subroutine jac
  
  