! WARNING: THIS CODE IS DEPRECATED AS OF 2016.01.08 --JSh

! 2008.10.20 - j. shanks
! program to simulate data for both ORM and gain map analysis
! using the same set of misalignments for both data sets

! 2009.09.24 - This code is HEAVILY outdated, and will not run as-is.
! Probably another 1-2hrs debugging necessary to revive this program. -JSh

! 2009.10.14 - Begin overhaul. -JSh

program sim_orm

  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_decay_mod

 implicit none

 

    ! variables used in ORM data generation:
    type (lat_struct) ring
    type (coord_struct), allocatable :: orb(:)
    type(normal_modes_struct) modes
    type (db_struct) cesr
    type(ma_struct) ma_params(20)
    type(det_struct), allocatable :: bpm(:)
    type(det_error_struct) :: bpm_error_sigmas
    type(meas_beam_struct) :: beam
    type(orm_data_struct), allocatable :: orm_data(:)
    character(40) :: bpm_mask = "^DET\_[0-9]{2}[ewEW]$" ! default to using CESR naming convention


    real(rp) ::  mag_hkick = 1.e-7, mag_vkick = 1.e-7 !magnitude of the kick; identical at each kicker
    integer, allocatable :: hkick_ix(:), vkick_ix(:) ! array of element numbers of the steering magnets
    logical err
    logical :: misalign_magnets = .false.

    integer ios, version, arg_num, iargc, readstatus, rad_cache
    
    integer i, j, jmax, k, ix, jx, kx, bpm_index, dummy
    integer :: steering_index = 1
    integer :: n_bpms = 0, n_steering_tot = 0, n_steering_h = 0, n_steering_v = 0 !break into horizontal and vertical, for tallying
    real(rp) :: bpm_offset = 0, bpm_rotation = 0, bpm_noise = 0, eta_error = 0, single_bpm_tilt = 0.
    real(rp) harvest, x, y

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

		! variables for gain map data generation:
    real(rp) :: shear_x = 0., current = 7500 ! amount of horizontal shearing on bpm's, in meters
    real(rp) b(4), d(4,0:2,0:2) ! four button responses, for each (x,y) combination
    real(rp) :: butn_res_sigma = 0., gain_sigma = 0., timing_sigma=0. ! butn_res_sigma = position uncertainty/resolution, in units of button signal


  namelist /parameters/ lat_file, output, bpm_offset, bpm_rotation, mag_hkick, mag_vkick, &
                        ma_params, bpm_noise, misalign_magnets, shear_x, current, &
                        gain_sigma, timing_sigma, bpm_mask

  write(*,'(a)') 'WARNING: THIS CODE IS DEPRECATED AS OF 2016.01.08 --JSh'

  !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_orm.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

  output = trim(output)
  if((output .ne. 'loco') .and. (output .ne. 'tao_cesr')) then
     write(*,'(a,a)') trim(output), " is not a valid output format. Use either 'loco' or 'tao_cesr'."
     stop
  endif

  ! 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
  
  ring%ele(:)%value(hkick$) = 0
  ring%ele(:)%value(vkick$) = 0 
 
  call set_on_off(rfcavity$, ring, off$)
  bmad_com%radiation_damping_on = .false.

  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)

  beam%current = current
  beam%I0 = current

  ! find all horiz, vert steerings:
  call find_steerings(cesr, hkick_ix, vkick_ix)
  n_steering_h = ubound(hkick_ix,1)
  n_steering_v = ubound(vkick_ix,1)
  n_steering_tot = n_steering_h + n_steering_v

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



  !=== 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", "g1", "g2", "g3", "g4"
  write(47, '(a7, 7a14)') ""
  
  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(:)
  enddo
  close(unit=47)
  

! directly from ring_ma:

  write(*,*) "Introducing misalignments..."

  dr_misalign_params%sigma_cutoff = 3.
  dr_misalign_params%alignment_multiplier = 1
  if (misalign_magnets .eqv. .true.) then
     call dr_misalign(ring, ma_params)
  endif

  open(unit=31, file='quad.ma', status='replace')
  open(unit=32, file='sext.ma', status='replace')
  open(unit=33, file='bend.ma', status='replace')

  write (31, '(a7,a7,3a14)') "name", "index", "x-offset", "y-offset", "tilt"
  write (32, '(a7,a7,3a14)') "name", "index", "x-offset", "y-offset", "tilt"
  write (31, '(a7,a7,3a14)') ""
  write (32, '(a7,a7,3a14)') ""
  write (33, '(a7,a7,2a14)') "name", "index", "x-offset", "y-offset"
  write (33, '(a7,a7,3a14)') ""
  
   do i=0, 99 ! write quad misalignments to output
      if (cesr%quad(i)%ix_lat == 0) cycle
      write(31, '(a7,3e14.5)') ring%ele(cesr%quad(i)%ix_lat)%name, ring%ele(cesr%quad(i)%ix_lat)%value(x_offset$), &
           ring%ele(cesr%quad(i)%ix_lat)%value(y_offset$), ring%ele(cesr%quad(i)%ix_lat)%value(tilt$)
   enddo
   do i=0, 99 ! write sextupole misalignments to output
      if (cesr%sex(i)%ix_lat == 0) cycle
      write(32, '(a7,3e14.5)') ring%ele(cesr%sex(i)%ix_lat)%name, ring%ele(cesr%sex(i)%ix_lat)%value(x_offset$), &
           ring%ele(cesr%sex(i)%ix_lat)%value(y_offset$)
   enddo

   do i=1, ring%n_ele_max
      if (match_reg(ring%ele(i)%name, 'RAW_BET.*') .or. match_reg(ring%ele(i)%name, 'BUMP.*')) cycle
      if (match_reg(ring%ele(i)%name, '^B[0-9]*')) then
         write(33, '(a7,2e14.5)') ring%ele(i)%name, ring%ele(i)%value(x_offset$), ring%ele(i)%value(y_offset$)
      endif
   enddo

   close(unit=31)
   close(unit=32)
   close(unit=33)


   ! write mdsa_fakedata.in for tao_cesr
   if(output == 'tao_cesr') then
      open(unit=47,file='mdsa_fakedata.in', status='replace')
      write(47,'(a,a,a)') "&lattice file_name = '$BMAD_LAT/", trim(lat_file_name),".lat' /"
      do i=1, n_steering_tot
         write(47,'(a,i3.3,a,i3.3,a)') "&measurement who = orbit 'orbit_files/orbit.1",2*i,"' 'orbit_files/orbit.1",2*i-1, "' /"
      enddo
      close(unit=47)
   endif

   ! prep orm_data_struct for passing information to sim_orm_data:

   if (.not. allocated(orm_data)) allocate(orm_data(n_steering_tot*2))

   do jx = 1, ubound(hkick_ix,1)
      orm_data(2*jx-1:2*jx)%mag_hkick = mag_hkick
      orm_data(2*jx-1:2*jx)%mag_vkick = 0.
      orm_data(2*jx-1:2:jx)%hkick_ix = hkick_ix(jx)
      orm_data(2*jx-1:2:jx)%ele_name = ring%ele(hkick_ix(jx))%name
   enddo
   do jx = 1, ubound(vkick_ix,1)
      kx = ubound(hkick_ix,1) + jx
      orm_data(2*kx-1:2*kx)%mag_hkick = 0.
      orm_data(2*kx-1:2*kx)%mag_vkick = mag_vkick
      orm_data(2*kx-1:2:kx)%vkick_ix = vkick_ix(jx)
      orm_data(2*kx-1:2:kx)%ele_name = ring%ele(vkick_ix(jx))%name
   enddo
   
   ! simulate all ORM data-- save in orm_data_struct
   call sim_orm_data(ring, bpm, beam, orm_data)
   
   
   do i=1, 2*n_steering_tot
      write(orbit_file, fmt='(A,i3.3)') 'orbit_files/orbit.1', i
      open(unit=47, file=orbit_file, status='replace')
      
      if(output=='tao_cesr') then !write headers for orbit files
         write(47,'(a)') "&data_parameters      ! Fortran90 namelist."
         write(47,'(a)') "param%data_type = 'ORBIT'   ! Identifies file as an orbit file."
         write(47,'(a,a,a)') "param%lattice   = '", trim(lat_file_name), "'"
         !write(47,'(a)') "param%data_date = '1998-JAN-06 18:05:36'"
         !write(47,'(a)') "param%comment   = 'This optional comment is used for a title when plotting'"

         if(i .lt. n_steering_h*2+1) then
            write(47,'(a,a,a)') "param%var_attrib_name = '", "HKICK", "'"
            write(47,'(a,a,a)') "param%var_ele_name    = '", trim(orm_data(i)%ele_name),"'    !name of element being varied"
            write(47,'(a,es14.4,a)') "param%dvar            = ",  mag_hkick     , "        ! radians"
         else
            write(47,'(a,a,a)') "param%var_attrib_name = '", "VKICK", "'"
            write(47,'(a,a,a)') "param%var_ele_name    = '", trim(orm_data(i)%ele_name),"'    !name of element being varied"
            write(47,'(a,es14.4,a)') "param%dvar            = ",  mag_vkick     , "        ! radians"
         endif
         
         write(47,'(a)') "/"
         write(47,'(a)') ""
         write(47,'(a)') "&cesr_data"
      endif
      
      do j=1, n_bpms
         if(output=='loco') then ! output for LOCO is in meters
            write(47,'(i6, 2f12.7)') bpm(j)%ix_db, orm_data(i)%bpm(j)%vec(1), orm_data(i)%bpm(j)%vec(3)
         elseif(output=='tao_cesr') then ! output for tao_cesr must be in meters
            write(47,'(a,i2,a,2es18.5)') "   orbit(", bpm(j)%ix_db, ") =", orm_data(i)%bpm(j)%vec(1), orm_data(i)%bpm(j)%vec(3)
         endif
      enddo
      
      if(output=='tao_cesr') then
         write(47,'(a)') "/"
         write(47,'(a)') ""
      endif
      
      close(unit=47)
      
   enddo  ! writing ORM data to files.
   
 end program sim_orm