program turn_by_turn
! turn by turn tracking with radiation excitation and damping
! and modulated kicker (shaker) in x,y, pr z plane
! Write out turn by turn position and phase word
! Write out phase and amplitude

   use bmad
   use cesr_basic_mod
   use rad_int_common
   use radiation_mod
   use normal_mode_decomposition_mod
   use nr
   use ran_state, only: ran_seed
 
     implicit none

   interface
     subroutine save_last_turns(n,m,end,saved_turns)
      use bmad
      implicit none
      type (coord_struct) end
      type (coord_struct) saved_turns(:)
      integer n, m
     end subroutine
   end interface

   interface
    subroutine fft_plane(n,co, Q, phase)
    use bmad
    implicit none
    type(coord_struct) co(:)
    integer n
    real(rp) Q(3), phase(3)
    end subroutine
   end interface

      type (lat_struct) lat
      type (cesr_struct) cesr
      type (ele_struct) ele 
      type (coord_struct), allocatable :: orbit(:), high(:), low(:), saved_turns(:)
      type (coord_struct) start, end, bpm_resolution, start_offset
      type (normal_modes_struct) mode
      type (coord_struct) traj_det(0:10000, -1:n_det_maxx)
      type phase_amp_struct
       real(rp) horz_sin_sum(-1:n_det_maxx), horz_cos_sum(-1:n_det_maxx) 
       real(rp) vert_sin_sum(-1:n_det_maxx), vert_cos_sum(-1:n_det_maxx) 
       real(rp)long_sin_sum(-1:n_det_maxx), long_cos_sum(-1:n_det_maxx)
       real(rp)energy_sin_sum(-1:n_det_maxx), energy_cos_sum(-1:n_det_maxx)
       real(rp)sin_phi, cos_phi
       real(rp) horz_delta_phi(-1:n_det_maxx),vert_delta_phi(-1:n_det_maxx),long_delta_phi(-1:n_det_maxx)
       real(rp) horz_amplitude(-1:n_det_maxx), vert_amplitude(-1:n_det_maxx), long_amplitude(-1:n_det_maxx)
       real(rp) alt_horz_delta_phi(-1:n_det_maxx),alt_vert_delta_phi(-1:n_det_maxx),alt_long_delta_phi(-1:n_det_maxx)
       real(rp) alt_horz_amplitude(-1:n_det_maxx), alt_vert_amplitude(-1:n_det_maxx), alt_long_amplitude(-1:n_det_maxx)
       real(rp) energy_delta_phi(-1:n_det_maxx), energy_amplitude(-1:n_det_maxx)
       real(rp) alt_energy_delta_phi(-1:n_det_maxx), alt_energy_amplitude(-1:n_det_maxx)
       real(rp) phi
       real(rp) sum

       integer unit

       character*19 file_name
       character*4 type
      end type

      type (phase_amp_struct) plane(3) 


      integer ix, i, j, m, ii, n_turns, index, ixx, k
      integer n_save
      integer ix_cache
      integer nargs,iargc
      integer num_det, ix_x(100),  jx
      integer, allocatable :: ix_ele(:)
      integer iu,iu1
      integer ix_beg, ix_end
      integer res,np
      integer ios
      integer np1
      integer time(8)
      integer phih, phiv, bpm
      integer iu2
      integer n_fft_turns, kick_plane_index
      integer ix_shake

      character*140 lat_file
      character*120 line, last_line, changes(1:3)/3*' '/
      character*7 det_name
      character*40 fake_data_params_file
      character*16 shake_ele/'RF_W1'/, kick_plane/'Z'/
      character*3 bpm_name

      real(rp) mat(6,6), W(6,6), Y(6,6), tune(3)
      real(rp), parameter :: c_q = 3.84e-13
      real(rp) energy, gamma2_factor, i2, i5b
      real(rp) delta_kick, kick/2.e-5/

      real(rp) A, B 
      real(rp) delta/1.e-4/
      real(rp) RadToDeg, MetersToMM
      real(rp) idum
      real(rp) s
      real(rp) bpm_res, synch_tune
      real(rp) sum_long_amp, sum_energy_amp
      real(rp) tune_fit, phase, amp, chi2
      real(rp) but(4),h_diff,v_diff
      real(rp) Q(3), ph(3)
      real(rp) tan_phi, phi_at_start(2)

      logical radiation

      RadToDeg = 360./twopi
      MetersToMM = 1000.

      plane(1)%file_name = 'phase_amp_x.fake'
      plane(2)%file_name = 'phase_amp_y.fake'
      plane(3)%file_name = 'phase_amp_z.fake'
      plane(1)%type = 'HORZ'
      plane(2)%type = 'VERT'
      plane(3)%type = 'LONG'

!      nargs = cesr_iargc()
!      if(nargs == 1)then
!         call cesr_getarg(1,lat_file)
!         print *, 'Using ', trim(lat_file)
!       else
!      lat_file = 'bmad.'
!      type '(a,$)',' Lattice file name ? (default= bmad.) '
!      read(5,'(a)') line
!       call string_trim(line, line, ix)
!       lat_file = line
!       if(ix == 0) lat_file = 'bmad.'
!       type *, ' lat_file = ', lat_file
!      endif



   call get_fake_data_params(n_turns,n_save, kick,bpm_res,synch_tune, &
            idum, radiation, changes, shake_ele, start_offset, n_fft_turns, kick_plane, lat_file)

       call bmad_parser (lat_file, lat)


   if(kick_plane == 'X')kick_plane_index = 1
   if(kick_plane == 'Y')kick_plane_index = 2
   if(kick_plane == 'Z')kick_plane_index = 3
   print '(1x,a11,1x,i6,1x,a11,1x,i6)',' n_turns = ',n_turns,'  n_save = ', n_save
   print '(1x,a11,e12.4,a14,e12.4)', ' bpm_res = ', bpm_res,' synch_tune = ',synch_tune
   print '(1x,a,a,a,e12.4)',' Amplitude of modulated ',kick_plane,'-plane kick = ', kick
   print '(1x,a36,l)',' Radiation damping and excitation = ', radiation

     if (idum < 1) then
        call date_and_time(values=time)
        idum = time(8)
     end if
     call ran_seed(time(8))
     print *,' idum = ', idum


      line = 'ALL H_SEP 0'
      call find_change( line, lat)
      do i = 1, 3
       if(changes(i) /= ' ')call find_change(changes(i), lat)
      end do
      do while(line(1:1) /= 'G')
10    print '(a, $)', '  element change or GO> '
      read(5, '(a)',err=10) line
     
  ix = index(line, '!')
  if (ix /= 0) line = line(:ix-1)        ! strip off comments

  call str_upcase(line, line)
  call string_trim(line, line, ix)

  if (ix == 0) then       ! nothing typed. do the same thing
      line = last_line
  endif

   last_line = line

   call str_upcase(line,line)
   if(line(1:1) /= 'G')call find_change( line, lat)
   end do


       call reallocate_coord(orbit,lat%n_ele_max)
       call reallocate_coord(high,lat%n_ele_max)
       call reallocate_coord(low,lat%n_ele_max)
       call twiss_at_start(lat)

       call twiss_propagate_all(lat)
       call bmad_to_cesr(lat, cesr)

       lat%z%tune = synch_tune*twopi
       call set_z_tune(lat)
       call calc_z_tune(lat)
       call transfer_matrix_calc(lat, .true., mat)
       call normal_mode_decomposition(mat, tune, Y,W, .true.)

       print *,' tune from calc_z = ',lat%z%tune,'  tune from normal_mode_decomposition = ',tune(3)
       print '(a10,e12.4,a10,e12.4)',' a_tune = ',lat%a%tune/twopi,' b_tune = ', lat%b%tune/twopi

        ix_beg = 0
        if(lat_file(1:4) == 'bmad')ix_beg = index(lat_file, '_')
        ix_end = index(lat_file, '.')   

   call element_locator(shake_ele, lat, ix_shake)
   if(ix_shake /= -1)then
        print '(a17,a16,a8,e12.4)',' Shake element = ', shake_ele, '  eta = ', lat%ele(ix_shake)%x%eta
    else
        print '(3a)',' Shake element ', trim(shake_ele), ' not found '
   endif

      do k=1,3
!      iu=10+k
      iu=lunget()
      open(unit=iu,file=plane(k)%file_name)
      print '(a,1x,i,a,a)',' Write phase and amplitude data for plane ',k,' to file = ',plane(k)%file_name
        write(iu,'(a16)')'  &LATTICE_NAME '
        write(iu,'(a23,a)')'     phase_data_file = ', plane(k)%file_name
        write(iu, *)'     lattice = ',"'",lat_file(ix_beg+1:ix_end-1),"'"
        write(iu,'(a19,f12.5)')' horiz_beta_freq = ',abs(lat%z%tune)/twopi*390.1
        write(iu,'(a19,f12.5)')' vert_beta_freq = ',abs(lat%b%tune)/twopi*390.1
!        do i = 1,3
!           if(changes(i) /= ' ')write(iu,'(a9,i1,a4,a)')' changes(',i,') = ',changes(i)
!        end do
        write(iu,'(a5)')' /END'
        write(iu,*)
          line = ''

          fake_data_params_file = "fake_data_params.dat"
          print *,' Params file = ', fake_data_params_file
          iu1 = lunget()
          open (iu1, file=fake_data_params_file, type='old', readonly, shared, iostat = ios)

          do while(index(line,'/END') == 0)
            read(iu1,'(a72)')line
            write(iu,'(a72)')line(1:72)
          end do
          rewind(iu1)
        close(iu1)          
        write(iu,*)
        write(iu,'(a17)')'  &PHASE_AMP_DATA'
        plane(k)%unit = iu
!        print *,' k, plane(k)%type, plane(k)%unit ', k,plane(k)%type, plane(k)%unit
      end do

  ! Find the indicies of all of the BPMs ix(1:100)

      allocate(ix_ele(0:lat%n_ele_max))
      jx=0
      do i = 0,n_det_maxx
       ix = cesr%det(i)%ix_lat
       call element_locator(cesr%det(i)%name,lat, ix)
        if(ix <= 0)cycle
!       print '(a25,2i,1x,a12,1x,a12)',  ' i, ix, , cesr%det(i)%name,lat%ele(ix)%name ', i,ix,cesr%det(i)%name,lat%ele(ix)%name
         ix_ele(ix) = i
      end do

      high(0)%vec(6)=delta
      call closed_orbit_calc(lat,high,4,1)
      low(0)%vec(6)=-delta
      call closed_orbit_calc(lat,low,4,1)
      call closed_orbit_calc(lat, orbit, 6, 1)
      print '(a11,6e12.4)', ' Orbit%vec ', orbit(0)%vec(1:6)
      print '(1x,a19,e)', 'Synchrotron tune = ', lat%z%tune/twopi
      if(radiation)then
         bmad_com%radiation_damping_on = .true.
         bmad_com%radiation_fluctuations_on = .true.
      endif
      start%vec(1:6) = orbit(0)%vec(1:6) + start_offset%vec(1:6)
     print '(a,6e12.4)',' start%vec ', start%vec(1:6)
!      if(.not. radiation)start%vec(6) = kick
      ii=0

     s=0
     np=0

     tune(3) = lat%z%tune
     tune(1) = lat%a%tune
     tune(2) = lat%b%tune
     call reallocate_coord(saved_turns,n_fft_turns + 100)

     j=0
      do i =1,n_turns ! damping time
       j = j+1
       do m = 1, lat%n_ele_track
        if(m==1 .and. (i/1000)*1000 == i)print '(1x,a7,1x,i,1x,a16,6e12.4)',' start ',i, lat%ele(m)%name, start%vec(1:6)
        call track1(start, lat%ele(m), lat%param, end)
         if(lat%param%lost)then
          if(lat%param%plane_lost_at == x_plane$) print '(a)',' Particle lost in x_plane '
          if(lat%param%plane_lost_at == y_plane$) print '(a)',' Particle lost in y_plane '
          stop
         endif
          if(lat%ele(m)%name == shake_ele)then
           call save_last_turns(n_fft_turns,i,end,saved_turns)
           if(radiation)then
            tan_phi = -lat%ele(ix_shake)%a%alpha - lat%ele(ix_shake)%a%beta &
                  * end%vec(2)/end%vec(1)
            phi_at_start(1) = atan(tan_phi)
            tan_phi = -lat%ele(ix_shake)%b%alpha - lat%ele(ix_shake)%b%beta & 
                  * end%vec(4)/end%vec(3)
            phi_at_start(2) = atan(tan_phi)
!             if(100*(j/100) == j)print '(3e12.4)',phi_at_start(1), sin(phi_at_start(1)), end%vec(2)/end%vec(1)
             if(kick_plane == 'Z')end%vec(6) = kick*sin(-tune(3) * j + ph(3)) + end%vec(6)
             if(kick_plane == 'X')end%vec(2) = kick*sin(tune(1) * i + ph(1)) + end%vec(2)
             if(kick_plane == 'Y')end%vec(4) = kick*sin(tune(2) * i + ph(2)) + end%vec(4)
!             if(kick_plane == 'X')end%vec(2) = kick*sin(phi_at_start(1)) + end%vec(2)
!             if(kick_plane == 'Y')end%vec(4) = kick*sin(phi_at_start(2)) + end%vec(4)
            write(13,'(3e12.4)')end%vec(1:2), sin(tune(1) * j + ph(1))
          endif

        endif

        if(i >= n_turns - n_save .and.  any(cesr%det(0:n_det_maxx)%ix_lat == m))then
         ii=i-n_turns+n_save


            do res = 1,6
             if(bpm_res /= 0.)call gasdev(idum)
              bpm_resolution%vec(res) = idum*bpm_res
             if(np < 10)print '(1x,e12.4)',bpm_resolution%vec(res)
             np = np+1
             s = s+ bpm_resolution%vec(res)**2
            end do
          if(ii > 10000)then
            print *,' Cannot save more than 10000 turns. TRAJ_DET array limited to 10000 '
            stop
          endif
          traj_det(ii,ix_ele(m))%vec(1:6) = end%vec(1:6) + bpm_resolution%vec(1:6)
        endif

        start%vec(1:6) = end%vec(1:6)
           if(ix_shake <0 .and. m == lat%n_ele_track)call save_last_turns(n_fft_turns,i,end,saved_turns)
       end do


        if((i/100)*100 == i .and. i > n_fft_turns)then
            call fft_plane(n_fft_turns,saved_turns, Q, ph)
            tune(kick_plane_index) = Q(kick_plane_index) * twopi
            j = n_fft_turns
         endif
      if((i/1000)*1000 == i .and. i > n_fft_turns) print *,' Q = ', Q(kick_plane_index)
      end do

      plane(1:3)%sum=0.
      do j =0,n_det_maxx
!       if(cesr%det(j)%ix_lat == 0)cycle
        do k=1,3
         plane(k)%horz_cos_sum(j)=0.
         plane(k)%horz_sin_sum(j)=0.
         plane(k)%vert_cos_sum(j)=0.
         plane(k)%vert_sin_sum(j)=0.
         plane(k)%long_cos_sum(j)=0.
         plane(k)%long_sin_sum(j)=0.
         plane(k)%energy_cos_sum(j)=0.
         plane(k)%energy_sin_sum(j)=0.
        end do
      end do
      do i = 1,ii
       plane(1)%phi = (i+n_turns-n_save)*tune(1)
       plane(2)%phi = (i+n_turns-n_save)*tune(2)
       plane(3)%phi = (i+n_turns-n_save)*tune(3)
       do k=1,3
         plane(k)%cos_phi = cos(plane(k)%phi)
         plane(k)%sin_phi = sin(plane(k)%phi)
         plane(k)%sum = plane(k)%sum + plane(k)%cos_phi**2 + plane(k)%sin_phi**2
       end do
 
       do j = 0, n_det_maxx
        do k=1,3
         plane(k)%horz_cos_sum(j) = plane(k)%horz_cos_sum(j) + traj_det(i,j)%vec(1) * plane(k)%cos_phi
         plane(k)%horz_sin_sum(j) = plane(k)%horz_sin_sum(j) + traj_det(i,j)%vec(1) * plane(k)%sin_phi
         plane(k)%vert_cos_sum(j) = plane(k)%vert_cos_sum(j) + traj_det(i,j)%vec(3) * plane(k)%cos_phi
         plane(k)%vert_sin_sum(j) = plane(k)%vert_sin_sum(j) + traj_det(i,j)%vec(3) * plane(k)%sin_phi
         plane(k)%long_cos_sum(j) = plane(k)%long_cos_sum(j) + traj_det(i,j)%vec(5) * plane(k)%cos_phi
         plane(k)%long_sin_sum(j) = plane(k)%long_sin_sum(j) + traj_det(i,j)%vec(5) * plane(k)%sin_phi
         plane(k)%energy_cos_sum(j) = plane(k)%energy_cos_sum(j) + traj_det(i,j)%vec(6) * plane(k)%cos_phi
         plane(k)%energy_sin_sum(j) = plane(k)%energy_sin_sum(j) + traj_det(i,j)%vec(6) * plane(k)%sin_phi
        end do !k x,y,z

       end do !j  detectors
     end do !ii   turns

   
      print '(a6,3e12.4)',' sum= ',plane(1:3)%sum

! write turn by turn data

    iu2 = lunget()
    open(unit=iu2, file = 'turn_by_turn.dat')
    do j=1,98
     bpm = j
     if(j > 49)bpm = 99-j
     bpm_name(1:3)='   '
     if(bpm > 9)write(bpm_name(1:2),'(i2)')bpm
     if(bpm <= 9)write(bpm_name(2:2),'(i1)')bpm
     if(j <=49)bpm_name(3:3)='W'
     if(j > 49)bpm_name(3:3) = 'E'
     write(iu2,'(a21,a3)')'# Location         : ',bpm_name
     do i=1,ii  !turns
      h_diff = tune(1)
      v_diff = tune(2)
      phih = mod((i+n_turns-n_save)*tune(1),twopi)/twopi * 512
      phiv = mod((i+n_turns-n_save)*tune(2),twopi)/twopi * 512
      but(1:4) = 100
      write(iu2,'(4f13.4,2x,2f11.6,i12,1x,2i8,2f13.6)')but(1:4),traj_det(i,j)%vec(1),traj_det(i,j)%vec(3),i,phih,phiv, h_diff, v_diff    
     end do
   end do

     sum_long_amp = 0
     sum_energy_amp = 0
     np1 = 0
     tune_fit = lat%z%tune/twopi

     do k=1,3
      print *,' k = ', k,'  plane(k)%unit = ', plane(k)%unit
      do j = 0, n_det_maxx

         A = plane(k)%horz_cos_sum(j)/plane(k)%sum
         B = plane(k)%horz_sin_sum(j)/plane(k)%sum
         plane(k)%horz_delta_phi(j) = -atan2(B,A)   *RadToDeg
         plane(k)%horz_amplitude(j) = sqrt(A**2+B**2) *MeterstoMM
         A = plane(k)%vert_cos_sum(j)/plane(k)%sum
         B = plane(k)%vert_sin_sum(j)/plane(k)%sum
         plane(k)%vert_delta_phi(j) = -atan2(B,A)   *RadToDeg
         plane(k)%vert_amplitude(j) = sqrt(A**2+B**2) *MetersToMM

         A = plane(k)%long_cos_sum(j)/plane(k)%sum
         B = plane(k)%long_sin_sum(j)/plane(k)%sum
         plane(k)%long_delta_phi(j) = -atan2(B,A)  *RadToDeg
         plane(k)%long_amplitude(j) = sqrt(A**2+B**2) *MetersToMM
         if(k == 3)sum_long_amp = sum_long_amp + plane(k)%long_amplitude(j)

         A = plane(k)%energy_cos_sum(j)/plane(k)%sum
         B = plane(k)%energy_sin_sum(j)/plane(k)%sum
         plane(k)%energy_delta_phi(j) = -atan2(B,A)  *RadToDeg
         plane(k)%energy_amplitude(j) = sqrt(A**2+B**2) *MetersToMM
         if(k == 3)then
            sum_energy_amp = sum_energy_amp + plane(k)%energy_amplitude(j)
            np1 = np1+1
         endif
         write(plane(k)%unit,'(1x,a4,a1,i3,a2,6f12.4)')plane(k)%type,'(',j,')=',&
                                           plane(k)%horz_amplitude(j),plane(k)%vert_amplitude(j), &
                                           plane(k)%horz_delta_phi(j),plane(k)%vert_delta_phi(j), &
                                           plane(k)%long_amplitude(j),plane(k)%long_delta_phi(j)
      end do ! j detectors
      end do ! k x,y,z
      do k =1,3
       write(plane(k)%unit,'(a5)')' /END'
       write(plane(k)%unit,'(a8)')'  &FACTS'
       write(plane(k)%unit, '(1x,a22,e12.4)')' standard_deviation = ',sqrt(s/np)
       write(plane(k)%unit, '(1x,a36,l)')' Radiation_damping_and_excitation = ',radiation
       write(plane(k)%unit, '(1x,a34,e12.4,1x,a5)')' Average_longitudinal_amplitude = ',sum_long_amp/np1, ' ! mm'
       write(plane(k)%unit, '(1x,a34,e12.4,1x,a20)')' Average_energy_amplitude = ',sum_energy_amp/np1,' ! Delta E/E X 1000 '
       write(plane(k)%unit,'(a5)')' /END'
     print *,' k, plane(k)%unit ', k, plane(k)%unit
       close(unit=plane(k)%unit)
      end do ! k x,y,z
      print '(1x,a22,e12.4)',' standard deviation = ',sqrt(s/np)
      print '(1x,a36,l)',' Radiation damping and excitation = ',radiation
      print '(1x,a34,e12.4,1x,a5)',' Average longitudinal amplitude = ',sum_long_amp/np1, ' ! mm'
      print '(1x,a34,e12.4,1x,a20)',' Average energy amplitude = ',sum_energy_amp/np1,' ! Delta E/E X 1000 '
      print '(1x,a17,a)',' Shake element = ', shake_ele

  end

 subroutine get_det_ref(lat,m,ixx)
   use bmad
   implicit none
   type (lat_struct) lat
   integer m,ixx, index

         ixx = -1
         if(lat%ele(m)%name(8:8) /= ' ')return
         read(lat%ele(m)%name(5:6),'(i2)')index
         if(lat%ele(m)%name(7:7) == 'W')ixx = index
         if(lat%ele(m)%name(7:7) == 'E')ixx = 99- index

   return
  end

    subroutine get_fake_data_params(n_turns,n_save, kick,bpm_res,synch_tune,  idum, radiation, changes, shake_ele, &
                                            start_offset, n_fft_turns, kick_plane, lat_file)
     use bmad
     implicit none
      integer  n_turns, n_save, iu, ios, n_fft_turns
      character*140 lat_file     
      character*120  changes(1:3)
      character*20 fake_data_params_file
      character*16 shake_ele, kick_plane
      type (coord_struct) start_offset

      real(rp) kick, bpm_res, synch_tune, idum

      logical radiation

      namelist /fake_data_params/n_turns,n_save, kick,bpm_res,synch_tune, &
             idum, radiation, changes, shake_ele, kick_plane, start_offset, n_fft_turns, lat_file

! default
   n_fft_turns = 1024
!

   iu = lunget()
   fake_data_params_file = "fake_data_params.dat"
   print *,' Params file = ', fake_data_params_file
   open (iu, file=fake_data_params_file, type='old', readonly, shared, iostat = ios)
   read(iu, nml = fake_data_params, iostat = ios)

   if(ios /= 0)then
     print *,' Warning: Error reading FAKE_DATA namelist'
     stop
     rewind(iu)
   endif
   rewind(iu)

   close(iu)
  return
  end

subroutine save_last_turns(n,m,end,saved_turns)
 use bmad
 implicit none

 type (coord_struct) end
 type (coord_struct) saved_turns(:)

 integer n, m
 integer i

 if(m <= n)then
  saved_turns(m) = end
 else  !shift all turns down one
  saved_turns(n+1) = end
  do i=1,n
   saved_turns(i) = saved_turns(i+1)
  end do
!  print '(a,11e12.4)',' Save_last ',saved_turns(n-10:n)%vec(1)
endif
 
 return
end