! 2009.05.13 - JSh
! updated 2009.10.14 - modularized


program sim_gain_map

  use bmad
  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_decay_mod

  implicit none
  
  type (lat_struct) ring
  type(det_struct), allocatable :: bpm(:)
  type(det_error_struct) :: bpm_error_sigmas
  type(meas_beam_struct) :: beam
  type(gain_map_data_struct), allocatable :: gain_map_data(:)
  character(40) :: bpm_mask = "^DET\_[0-9]{2}[ewEW]$" ! default to using CESR naming convention

  real(rp) ::  dx = 1.e-3, dy = 1.e-3, b(1:4)
  integer :: nx = 3, ny = 3 !nx, ny = number of x,y data points; dx,dy = grid spacing
  logical err
  integer ios, version, arg_num, iargc, readstatus
  
  integer i, j, ix, jx, bpm_index, dummy
  real(rp) :: bpm_noise = 0, button_offset(4) = 0., tau = -1 ! tau = lifetime
  real(rp) harvest
  
  character(100) lat_file, file_name, lat_file_name, path, output_file_name
!  integer time(8)
  integer :: ran_seed = 0., debug_bpm = 15, n_bpms ! debug_bpm = single BPM index
  
  real(rp) :: current = 7500, shear_x = 0. ! amount of horizontal shearing on bpm's, in meters
  real(rp) :: gain_sigma = 0., bpm_offset=0, bpm_rotation=0, t, timing_sigma=0.
  real(rp) :: fixed_gain(1:4) = 1.
  

  namelist /parameters/ lat_file, dx, dy, nx, ny, bpm_noise, shear_x, current, &
       gain_sigma, button_offset, debug_bpm, tau, fixed_gain, 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_gain_map.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

  ! 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

  if (debug_bpm .eq. -2) then
     ! 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
     
     ! Locate BPMs we wish to study:
     call find_bpms(ring, bpm_mask, bpm)
     n_bpms = ubound(bpm,1)
     
  else
     allocate(bpm(1))
     bpm(1)%ix_db = debug_bpm
  endif

  do j=1,size(bpm)
     bpm(j)%gain(1:4) = fixed_gain(:)
  enddo
  
  dummy = splitfilename(lat_file,path,lat_file_name)
  lat_file_name = lat_file_name(1:len_trim(lat_file_name)-4)

  beam%current = current
  beam%I0 = current
  beam%tau = tau


  call string_trim (file_name, file_name, ix)
  file_name = file_name(:ix-3)
  write(output_file_name,'(3a)') 'bpms_', trim(file_name), '.ma'
  open(unit=47, file=trim(output_file_name), status='replace')
  write(47, '(a7,7a14)') "!index", "x-offset", "y-offset", "tilt", "g1", "g2", "g3", "g4"
  write(47, '(a7, 7a14)') ""
  do jx = 1, size(bpm)
     if ((fixed_gain(1) .eq. 1.) .and. (fixed_gain(2) .eq. 1.) .and. (fixed_gain(3) .eq. 1.) &
          .and. (fixed_gain(4) .eq. 1.)) then
        call bpm_errors(bpm(jx), bpm_error_sigmas)
     endif
     write(47, '(i7,7e14.4)') bpm(jx)%ix_db, bpm(jx)%x_offset, bpm(jx)%y_offset, bpm(jx)%tilt, bpm(jx)%gain(:)
  enddo
  close(unit=47)


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


  write(output_file_name,'(2a)') trim(file_name), '.input'
  open(unit=23,file=trim(output_file_name), status='replace')

  write(23,'(a)') "&bpm_gain_anal_params"
  write(23,'(a)') "    read_current_from_butns_file = .false."
  write(23,'(a)') "    read_fake_data = .true."
  write(23,'(a)') "    sigma_cutoff = 1000"
  write(23,'(a)') "    chisq_cutoff = 1000"

  if (debug_bpm .eq. -2) then
     write(23,'(a)') "    fit_bpm_idx(0) = -2" ! fit all BPMs
  else
     write(23,'(a,i3)') "    fit_bpm_idx(0) = ", bpm(1)%ix_db
  endif

  if (.not. allocated(gain_map_data)) allocate(gain_map_data(nx*ny))

  do jx = 1, size(bpm)

     call sim_gain_map_data(bpm(jx), beam, nx, ny, dx, dy, gain_map_data)
     
     do j=1, nx*ny
        write(23, '(a,i2,a,i2)') "    orbit(",j,")%butns_index =", j
        t = (30.*60/(nx*ny))*j  ! assume 30 mins for measurement
        call sim_decay(beam, t) 
        write(23, '(a,i2,a,es12.3)') "    orbit(",j,")%current = ", beam%current
        write(23,'(a,i2,a,i2,a,4f12.3)') "    orbit(",j,")%butns%det(", bpm(jx)%ix_db,")%amp = ", &
             gain_map_data(j)%amp(1:4)
        write(23,'(a,i2,a,i2,a,es12.3)') "    orbit(",j,")%butns%det(",bpm(jx)%ix_db,")%vec(1) = ", &
             gain_map_data(j)%vec(1)
        write(23,'(a,i2,a,i2,a,es12.3)') "    orbit(",j,")%butns%det(",bpm(jx)%ix_db,")%vec(3) = ", &
             gain_map_data(j)%vec(3)
     enddo
     
  enddo
 
  write(23,'(a)') "/"
  close(23)
  


  open(unit=47, file='fort.47', status='replace')
  write(47, '(11a12,2a30)') '!Orbit#', 'x (m)', 'y (m)', 'g1', 'g2', 'g3', 'g4', 'b1', 'b2', 'b3', 'b4', &
       'b1-b2-b3+b4', '(b1-b2+b3-b4)(b1+b2-b3-b4)'

  do j=1, nx*ny
     b(:) = gain_map_data(j)%amp(1:4)
     write(47, '(i12,10es12.3, 2es30.3)') j, gain_map_data(j)%vec(1), gain_map_data(j)%vec(3), &
          bpm(1)%gain(:), b(:), (b(1)-b(2)-b(3)+b(4)), (b(1)-b(2)+b(3)-b(4))*(b(1)+b(2)-b(3)-b(4))
  enddo
  close(47)
   



end program sim_gain_map