! 2009.07.08 - j. shanks
! program to simulate dispersion data 

program sim_eta

  use bmad
  use bmadz_interface
  use random_mod
  use dr_misalign_mod
  use radiation_mod
  use cesr_basic_mod
  use sim_utils
  use nonlin_bpm_mod
  use sim_bpm_mod
  use sim_eta_mod

 implicit none


    type (lat_struct) ring
    type (coord_struct), allocatable :: orb(:), eta(:)
    type(normal_modes_struct) modes
    type(det_struct), allocatable :: bpm(:)
    type(det_error_struct) :: bpm_error_sigmas
    type (all_pointer_struct), allocatable, save :: vkick_pointer(:)

    logical err
    integer ios, version, arg_num, iargc, readstatus, rad_cache
    
    integer i, j, ix, kx, bpm_index, dummy, n_iter, jx, ix_lat
    integer :: n_bpms = 0
    real(rp) :: bpm_offset = 0, bpm_rotation = 0, bpm_noise = 0
    real(rp) :: harvest, x, y, mag_vkick = 0.! vertical kick at 10w to generate vertical dispersion

    character(100) lat_file, file_name, lat_file_name,path
    character(50) orbit_file
    integer :: ran_seed = 0.

    real(rp) :: shear_x = 0., current = 7500 ! amount of horizontal shearing on bpm's, in meters
    real(rp) :: butn_res_sigma = 0., gain_sigma = 0., timing_sigma = 0. ! butn_res_sigma = position uncertainty/resolution, in units of button signal
    real(rp) :: delta, etaysq(1:99)  ! for manual dispersion 'measurement'
    real(rp), allocatable :: etay_rms(:)
    character(40) :: bpm_mask = "^DET\_[0-9]{2}[ewEW]$" ! default to using CESR naming convention



  namelist /parameters/ lat_file, bpm_offset, bpm_rotation, bpm_noise, &
       shear_x, current, gain_sigma, delta, n_iter, mag_vkick, timing_sigma, &
       bpm_mask

  !initialize randomizer's seed (harvest)
  ran_seed = 0. ! ran_seed = 0 <--> seed randomizer w/ current CPU time
  call ran_seed_put(ran_seed)

  arg_num=iargc()
  if(arg_num==0) then
     file_name='sim_eta.in'
  else
     call getarg(1, file_name)
  end if
  call string_trim (file_name, file_name, ix)
  open (unit= 1, file = file_name, status = 'old', iostat = ios)

  if(ios.ne.0)then
     write(*,*)
     write(*,*) 'ERROR: CANNOT OPEN FILE: ', trim(file_name)
  endif

  ! read in the parameters
  rewind(1)
  read(1, nml = parameters,iostat=readstatus)
  if(readstatus > 0) then
     print *,"CAN NOT READ FROM ",file_name
     stop
  end if

  allocate(etay_rms(1:n_iter))

  ! load parameters into bpm_error_sigmas structure:
  bpm_error_sigmas%bpm_offset   = bpm_offset
  bpm_error_sigmas%bpm_rotation = bpm_rotation
  bpm_error_sigmas%shear_x      = shear_x
  bpm_error_sigmas%gain_sigma   = gain_sigma
  bpm_error_sigmas%timing_sigma = timing_sigma

  ! init lattice
  if(lat_file(1:8) == 'digested')then
     call read_digested_bmad_file (lat_file, ring, version)
  else
     call fullfilename(lat_file,lat_file)
     call bmad_parser (lat_file, ring)
  endif

  dummy = splitfilename(lat_file,path,lat_file_name)
  lat_file_name = lat_file_name(1:len_trim(lat_file_name)-4)

! turn off all kicks / separators; commented out 2016.01.08 JSh
!  ring%ele(:)%value(hkick$) = 0
!  ring%ele(:)%value(vkick$) = 0


  call set_on_off(rfcavity$, ring, off$)
  call twiss_at_start(ring)
  call twiss_propagate_all(ring)
  call closed_orbit_calc(ring, orb, 4)


  ! Locate BPMs we wish to study:
  call find_bpms(ring, bpm_mask, bpm)
  n_bpms = ubound(bpm,1)


! turn on one steering to generate vertical dispersion of a predetermined amplitude:
  call pointers_to_attribute(ring, ring%ele(885)%name, "VKICK", .true., vkick_pointer, err)
  vkick_pointer(1)%r = vkick_pointer(1)%r + mag_vkick ! mag_kick is user-defined in the input file
  call closed_orbit_calc(ring, orb, 4)

!=== Misalignment ====================================================
     
  
  call resolution_to_button_error(bpm(1), current, bpm_noise)
  bpm(:)%butn_res_sigma = bpm(1)%butn_res_sigma

!=====================================================================
!=====================================================================
!===== Begin averaging loop! =========================================
!=====================================================================
!=====================================================================

  open(unit=31, file='rms_etay.out', status='replace')
  open(unit=32, file='eta.out', status='replace')
  write(31,'(a6,4a15)') '!kx',  'bpm_noise', 'bpm_rotation', 'gain_sigma', 'etay_rms(kx)'
  do kx = 1,n_iter ! number of iterations to average over
 
     !write(*,*) "kx = ", kx
!=== Generate BPM errors =============================================
! Note that putting values here have no effect on tracking. They're
! only for storage


     open(unit=47, file='bpms.ma', status='replace')
     write(47, '(a7,7a14)') "!index", "x-offset", "y-offset", "tilt", "g(1:4)"
     write(47, '(a7, 3a14)') ""
     
     do i = 1, n_bpms
        call bpm_errors(bpm(i),bpm_error_sigmas)
        write(47, '(i7,7e14.4)') bpm(i)%ix_db, bpm(i)%x_offset, &
             bpm(i)%y_offset, bpm(i)%tilt, bpm(i)%gain(:)
     end do
     close(unit=47)



     call sim_eta_data(ring, bpm, delta, eta)

     do ix = 1, n_bpms
        ix_lat = bpm(ix)%ix_lat
        etaysq(ix) = (eta(ix_lat)%vec(3))**2
     enddo
     etay_rms(kx) = sqrt(sum(etaysq)/n_bpms) ! RMS
     write(31,'(i6,4es15.3)') kx, bpm_noise, bpm_rotation, gain_sigma, etay_rms(kx)

  enddo ! matches [do kx=1, n_iter]

  
  do ix = 1, n_bpms
     ix_lat = bpm(ix)%ix_lat
     write(32,'(i6,5es15.3)') ix, bpm_noise, bpm_rotation, gain_sigma, eta(ix_lat)%vec(1), eta(ix_lat)%vec(3)
  enddo

  close(31)
  close(32)
  
  
end program sim_eta