!........................................................................
!
! Subroutine : FIND_JACOBIAN(T, NTRACKS, I_ELE, Y)
!
! Description:
!
! Arguments  : 
!          Input:
!             T       -- Bmadz_track_struct: Orbit of all of the tracks
!             NTRACKS -- Integer: number of tracks (5 or 7)
!             I_ELE   -- Integer: element number
!
!          Output:
!             Y       -- Real(rp): Place to stash the jacobian and 
!                        closed orbit information at element I_ELE         
!
! Mod/Commons:    
!
! Calls      :      
!
! Author     : 
!
! Modified   :     
!             
!........................................................................
!
!
! $Log$
! Revision 1.6  2007/01/30 16:15:13  dcs
! merged with branch_bmad_1.
!
! Revision 1.3  2004/05/03 17:58:24  dcs
! strack_struct -> bmadz_track_struct to avoid conflict.
!
! Revision 1.2  2003/04/30 17:14:50  cesrulib
! dlr's changes since last import
!
! Revision 1.1.1.1  2002/12/13 19:23:28  cesrulib
! import bmadz
!
!
!........................................................................
!


subroutine find_jacobian(t, ntracks, i_ele, y)

  use bmad


  use nonlin_mod

  implicit none

  type (bmadz_track_struct)  t(7)

  integer ntracks, i_ele, k, j, n
  real(rp) xx_in(7,7), xx_out(7,7), xx_in_inverse(7,7), y(7,7), temp(7,7)

  do j=1,ntracks
    xx_in(1:6,j) = t(j)%orb(0)%vec
    xx_in(ntracks,j) = 1.
    xx_out(1:6,j) = t(j)%orb(i_ele)%vec
    xx_out(ntracks,j) = 1.
  end do

!print *
!print *,' XX_IN'
!do j=1,ntracks
!  print '(7es12.4)',xx_in(j,1:7)
!end do


  n = ntracks
  call mat_inverse(xx_in(1:n,1:n), xx_in_inverse(1:n,1:n))
  y(1:n,1:n) = matmul (xx_out(1:n,1:n), xx_in_inverse(1:n,1:n))

temp(1:n,1:n) = matmul(xx_in_inverse(1:n,1:n),xx_in(1:n,1:n))
!print *
!print *,' XX * XX^{-1}'
!do j=1,ntracks
 ! print '(7es12.4)',temp(j,1:7)
!end do
!print *,' determinant (xx_in(1:5,1:5)) = ', determinant(xx_in(1:5,1:5))
!pause


end subroutine