!........................................................................
!
! Subroutine : GET_INIT_VEC(LAT, DEL_ACTION, NTRACKS, VEC_INIT)
!
! Description: subroutine to calculate the initial set of vectors to be 
!              propagated to compute jacobian. For fixed energy there are 
!              5 vectors. If there is RF and synchrotron oscillations there 
!              are 7 vectors
!
! Arguments  : 
!       Input:
!         LAT       -- lat_struct : Lat
!         DEL_ACTION -- Coord_struct : displacement from closed orbit in 
!                       action Units. With no coupling and alpha = 0 for 
!                       example:
!                        delta_x  = sqrt(del_action.pos) * sqrt(beta)
!                        delta_xp = sqrt(del_action.vel) / sqrt(beta)
!         NTRACKS    -- Integer: number of tracks used to compute jacobian
!
!       Output:
!         VEC_INIT  -- Select_mat_struct = Array of 5 or 7 starting vectors
!
! Mod/Commons:    
!
! Calls      :      
!
! Author     : 
!
! Modified   :     
!             
!........................................................................
!
!
! $Log$
! Revision 1.7  2007/01/30 16:15:13  dcs
! merged with branch_bmad_1.
!
! Revision 1.3  2003/06/05 18:33:27  cesrulib
! synch with bmad union removal
!
! 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 get_init_vec(lat, del_action, ntracks, vec_init)
  use bmad


  use nonlin_mod

  implicit none

  type (lat_struct)  lat
  type (coord_struct)  del_action
  type (select_mat)  mat(7)

  real(rp) g_mat(4,4), g_inv_mat(4,4), v_mat(4,4), v_inv_mat(4,4), gv_mat(6,6)
  real(rp) vec_init(7,6), product(6,6)
  integer ntracks
  integer i,j,k

  real(rp) matr(4,5) / 1.0, 0.0, 1.0, 0.0, &
           0.309017,  0.951057,  0.309017,  0.951057, &
          -0.809017,  0.587785, -0.809016, -0.587785, &
          -0.809017, -0.587785,  0.309017, -0.951056, &
           0.309017, -0.951056, -0.809017,  0.587785 /

!

  if (del_action%vec(1) == 0.) del_action%vec(1) = 2.e-7   ! default sigma
  if (del_action%vec(2) == 0.) del_action%vec(2) = del_action%vec(1)
  if (del_action%vec(3) == 0.) del_action%vec(3) = del_action%vec(1)
  if (del_action%vec(4) == 0.) del_action%vec(4) = del_action%vec(1)
  if (del_action%vec(5) == 0.) del_action%vec(5) = del_action%vec(1)
  if (del_action%vec(6) == 0.) del_action%vec(6) = del_action%vec(1)

  call make_v_mats (lat%ele(0), v_mat, v_inv_mat)
  call make_g_mats (lat%ele(0), g_mat, g_inv_mat)
  gv_mat(1:4,1:4) = matmul (v_mat, g_inv_mat)
  gv_mat(5,5) = 1
  gv_mat(6,6) = 1

  vec_init = 0

  do k=1,ntracks

    mat(k)%select(:,:) = 0
    do i=1,ntracks-1
      mat(k)%select(i,i) = matr(i,k)
    end do

    product = matmul (gv_mat, mat(k)%select)

    do i=1,ntracks-1
      do j=1,ntracks-1
        vec_init(k,i)= product(i,j)*sqrt(del_action%vec(j)) + vec_init(k,i)
      end do
    end do

  end do

  return

end subroutine