!spin with constant pitch
 program spin_check

   use bmad
!   use muon_mod
!   use muon_interface
!   use parameters_bmad
!   use quad_scrape_parameters   
      type (lat_struct),target:: lat
      type (branch_struct),pointer:: branch
      type (ele_struct) ele 
      type (coord_struct), allocatable :: co(:)
      type (coord_struct) orb_at_s 

      real(rp) w,s1,s2, half
      real(rp) offset
      real(rp) dist
      real(rp) f_rev
      real(rp) delta_e, chrom_x, chrom_y
      real(rp) delta_B/0./
      real(rp) pangle
      real(rp) s
      real(rp) px,py,pz
      real(rp) betahat(1:3)
      real(rp) beta_dot_s
      real(rp) deltap
      real(rp) xoff
      real(rp) voff
      real(rp) xp
      real(rp) theta, beta(1:3),beta_cart(1:3), spin_cart(1:3), sCrossBeta(1:3), sCrossB(1:3)
      real(rp) B(1:3)/0.,1.,0./
<<<<<<< .mine
      real(rp) gamma
=======
      real(rp) radius/7.112/
>>>>>>> .r49426

      integer ix, i, j, m, ii, n_turns, index, ixx, k
      integer nargs,iargc
      integer datetime_values(8)
      integer n
      integer ib/0/, track_state
      integer lun
            
      character*16 date, time, zone
      character*100 lat_file, pitch_angle/'0.5'/
      character*100 delta_p/'0.118034'/
      character*100 x_offset/'0.0'/
      character*100 v_offset/'0.0'/
      character*100 x_angle/'0.0'/

      logical err_flag

    call date_and_time(date,time,zone,datetime_values)
    print '(a10,a10,a1,a10)',' date = ',date,' time = ', time
!    directory = trim(date)//'_'//trim(time(1:6))
!    status= system('mkdir ' // directory)
!    if(status == 0)print '(a20,a)',' create directory = ', trim(directory)

!      OPEN (UNIT=5, FILE='quad_input.dat', STATUS='old', IOSTAT=ios)
!    READ (5, NML=quad_input, IOSTAT=ios)

!  print *, 'ios=', ios
!   rewind(unit=5)

!  CLOSE(5)
!  write(6,nml = quad_input)

!      quad_fringe_energy_change = .false.
!      quad_fringe_energy_change = .true.
      bmad_com%spin_tracking_on = .true.
      bmad_com%aperture_limit_on = .false.
      bmad_com%max_aperture_limit = 10000.
      nargs = cesr_iargc()
!      filename = 'test'
!      ix_branch = 0
      
      if(nargs >= 1)then
         call cesr_getarg(1,lat_file)
         print *, 'Using lattice ',trim(lat_file)
         if(nargs >= 2)then
            call cesr_getarg(2,x_offset)
            print *,'Using x offset ', trim(x_offset)
            read(x_offset,*)xoff
         endif
         if(nargs >= 3)then
            call cesr_getarg(3,x_angle)
            print *,'Using x angle ', trim(x_angle)
            read(x_angle,*)xp
         endif
         if(nargs >= 4)then
            call cesr_getarg(4,v_offset)
            print *,'Using v offset', trim(v_offset)
            read(v_offset,*)voff
         endif
         if(nargs >= 5)then
            call cesr_getarg(5,pitch_angle)
            print *,'Using pitch angle', trim(pitch_angle)
            read(pitch_angle,*)pangle
         endif
         if(nargs >= 6)then
            call cesr_getarg(6,delta_p)
            print *,'Using momentum offset ', trim(delta_p)
            read(delta_p,*)deltap
         endif
         
      else
                 
      lat_file = 'bmad.'
      read(x_offset,*)xoff
      read(pitch_angle,*)pangle
      read(x_angle,*)xp
      read(v_offset,*)voff
      read(delta_p,*)deltap
      print '(6a16)',lat_file,' x_offset = ',x_offset,' x_angle = ',x_angle,' v_offset = ', v_offset,' pitch_angle = ',pitch_angle,' delta_p = ',delta_p
   endif

       print '(2a)', ' lat_file = ', lat_file
       call bmad_parser (lat_file, lat)
       !       print *,' element offset = ', offset
       call reallocate_coord (co, lat%n_ele_track)
!       do ie=1,lat%branch(ix_branch)%n_ele_track
!          print *,ie,lat%branch(ix_branch)%ele(ie)%name
       !       end do

       call init_coord(co(0), particle=antimuon$)

       co(0)%vec=0
!       print '(a,es12.4)',' Pitch angle = ', pangle
       co(0)%vec(1) = xoff
       co(0)%vec(2) = xp
       co(0)%vec(4) = pangle
       co(0)%vec(6) = deltap
       co(0)%vec(3) = voff
       print '(a,6es12.4)',' vec = ', co(0)%vec

       n=lat%n_ele_track
       print '(a10, 6a12)','element','s','x','y','z','theta'
!       write(11, '(a10, 6a12)')'element','s','x','y','z','theta'
       do i=1,n
         print '(i10,6es12.4)',i,lat%ele(i)%s, lat%ele(i)%floor%r, lat%ele(i)%floor%theta
!         write(11, '(i10,6es12.4)')i,lat%ele(i)%s, lat%ele(i)%floor%r, lat%ele(i)%floor%theta
       end do
   print *,' lat%n_ele_track = ', lat%n_ele_track, n
   lun=lunget()
   open(unit=lun,file = 'tbt_'//trim(pitch_angle))
   print '(a,a)',' Write to file :','tbt_'//trim(pitch_angle)
   write(lun, '(2a10, 7a14,3a14, a14, a15,3a15, a15,9a14,3a14)') 'turn,','ele','s', 'x', 'xp', 'y', 'yp','z','zp', 's_x','s_y','s_z','beta_dot_s','time','beta_x','beta_y','beta_z','theta', &
        'beta_x cart','beta_y cart','beta_z cart','spin_x cart','spin_y cart','spin_z cart','sCrossBeta_x','sCrossBeta_y','sCrossBeta_z', &
        'spin X B x','spin X B y','spin X B z'
        co(0)%spin(:) = 0
        co(0)%spin(1) = 1. 
        print '(6es14.6,3es14.6)',co(0)%vec, co(0)%spin
   s=0
   do i=1,2000
        track_state = 0
        call track_all(lat,co, ib, track_state, err_flag)
        if(co(lat%n_ele_track)%state /= alive$)then
           print *,'Particle lost at element ', track_state
           print *,'Lost at ', coord_state_name(co(track_state)%state)
           print *,' co(track_state)%vec ', co(track_state)%vec
           exit
        else
        
          do j=1, lat%n_ele_track
            s = (i-1)*co(lat%n_ele_track)%s + co(j)%s
            px = co(j)%vec(2)
            py = co(j)%vec(4)
            pz = sqrt((co(j)%vec(6)+1)**2 - co(j)%vec(2)**2-co(j)%vec(4)**2)
            betahat(1:3)=(/px,py,pz/)/sqrt(px**2+py**2+pz**2)
            beta(1:3) = betahat(1:3) * co(j)%beta
            theta = co(j)%s/lat%ele(lat%n_ele_track)%s * twopi
            beta_cart = beta
            beta_cart(1) = beta(1)*cos(theta) - beta(3)*sin(theta)
            beta_cart(3) = beta(1)*sin(theta) +beta(3)*cos(theta)            
            beta_dot_s = dot_product(betahat,co(j)%spin)
            spin_cart  = co(j)%spin
            spin_cart(1) = co(j)%spin(1)*cos(theta) - co(j)%spin(3)*sin(theta)
            spin_cart(3) = co(j)%spin(1)*sin(theta) + co(j)%spin(3)*cos(theta)
            sCrossBeta(1) = spin_cart(2)*beta_cart(3) - spin_cart(3)*beta_cart(2)
            sCrossBeta(2) = spin_cart(3)*beta_cart(1) - spin_cart(1)*beta_cart(3)
            sCrossBeta(3) = spin_cart(1)*beta_cart(2) - spin_cart(2)*beta_cart(1)
            sCrossB(1) = -spin_cart(3)*B(2)
            sCrossB(2) = 0
            sCrossB(3) = spin_cart(1)*B(2)
            gamma = 1./sqrt(1-dot_product(beta_cart,beta_cart))
            
            
            write(lun, '(2i10,7es14.6,3es14.6,es14.6,es14.6,3es14.6, es14.6,9es14.6,3es14.6)')i,j,s,co(j)%vec, co(j)%spin, beta_dot_s, co(j)%t, betahat(1:3), theta, beta_cart, spin_cart, sCrossBeta, &
                 sCrossB
            write(11, '(i10,3es12.4)')j,co(j)%s,radius*cos(co(j)%s *twopi/lat%ele(lat%n_ele_track)%s),radius*sin(co(j)%s *twopi/lat%ele(lat%n_ele_track)%s)
          end do
          co(0) = co(lat%n_ele_track) 
        endif

   end do
   close(lun)


       print '(a,f12.4,a,f16.6)', ' The length is ',lat%ele(lat%n_ele_track)%s, ' gamma = ', gamma
end

!+
! Subroutine em_field_custom (ele, param, s_rel, orbit, local_ref_frame, field, calc_dfield, err_flag, &
!                                     calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles)
!
! Routine for handling custom (user supplied) EM fields.
! This routine is called when ele%field_calc = custom$ or when ele is a custom element (ele%key = custom$)
! In order to be used, this stub file must be modified appropriately. See the Bmad manual for more details.
!
! Note: Unlike all other elements, "s_rel" and "here" areguments for a patch element are with respect to 
! the exit reference frame of the element. See the Bmad manual for more details.
!
! Note: Fields should not have any unphysical discontinuities. 
! Discontinuities may cause Runge-Kutta integration to fail resulting in particles getting marked as "lost".
! The mode of failure here is that RK will try smaller and smaller steps to integrate through the 
! discontinuity until the step size gets lower than bmad_com%min_ds_adaptive_tracking. At this
! point the particle gets marked as lost.
!
! General rule: Your code may NOT modify any argument that is not listed as
! an output agument below.
!
! Input:
!   ele              -- Ele_struct: Custom element.
!   param            -- lat_param_struct: Lattice parameters.
!   s_rel            -- Real(rp): Longitudinal position relative to the start of the element.
!   orbit            -- Coord_struct: Coords with respect to the reference particle.
!   local_ref_frame  -- Logical, If True then take the 
!                          input coordinates and output fields as being with 
!                          respect to the frame of referene of the element. 
!   calc_dfield      -- Logical, optional: If present and True then the field 
!                          derivative matrix is wanted by the calling program.
!   use_overlap      -- logical, optional: Add in overlap fields from other elements? Default is True.
!   calc_potential   -- logical, optional: Calc electric and magnetic potentials? Default is false. 
!   grid_allow_s_out_of_bounds 
!                    -- logical, optional: For grids, allow s-coordinate to be grossly out of bounds 
!                          and return zero instead of an error? Default: False. Used internally for overlapping fields.
!   rf_time          -- real(rp), optional: Set the time relative to the RF clock. Normally this time is calculated using
!                          orbit%t or orbit%vec(5) but sometimes it is convenient to be able to override this.
!                          For example, time_runge_kutta uses this.
!   used_eles(:)     -- ele_pointer_struct, allocatable, optional: For use if this routine is called recursively. 
!                          This argument should be passed back to em_field_calc if em_field calc is called by this routine.
!
! Output:
!   field    -- Em_field_struct: Structure hoding the field values.
!   err_flag -- Logical, optional: Set true if there is an error. False otherwise.
!-

recursive subroutine em_field_custom (ele, param, s_rel, orbit, local_ref_frame, field, calc_dfield, err_flag, &
                                         calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles)

use em_field_mod, except_dummy => em_field_custom

implicit none

type (ele_struct) :: ele
type (lat_param_struct) param
type (coord_struct), intent(in) :: orbit
type (em_field_struct) :: field
type (ele_pointer_struct), allocatable, optional :: used_eles(:)

real(rp), intent(in) :: s_rel
real(rp), optional :: rf_time
real(rp) k
real(rp) x,z,length/44.686/
real(rp) Bl
real(rp) radius/7.112/

logical local_ref_frame
logical, optional :: calc_dfield, err_flag, calc_potential, grid_allow_s_out_of_bounds, use_overlap

character(*), parameter :: r_name = 'em_field_custom'

!
x=radius*cos(orbit%s *twopi/length)
z=radius*sin(orbit%s *twopi/length)
Bl=1.4513007*1.e-3
field%B= -z*Bl
field%B(2) = 1.4513007
field%E= x*Bl
!k= -1250000. (works for 0.1 radian)
k= -36000. !omega_a \sim omega_p
!k= -18000. !omega_a \sim 2*omega_p
k=0.
field%E(1) = -k * orbit%vec(1)
field%E(2) = k * orbit%vec(3)
!if (present(err_flag)) err_flag = .true.


end subroutine