program ring_ma2
  use bmad
  use correct_ring_mod2 ! modified version, using simulated data
  use dr_misalign_mod
  use ran_state, only: ran_seed
  use radiation_mod
  use cesr_basic_mod
  use nonlin_bpm_mod
  use sim_bpm_mod
  use tune_tracker_mod
  use tt_tools_mod
  use sim_ac_and_cbar_mod
  use sim_eta_mod
  use bsim_cesr_interface
  use ptc_layout_mod
  use mode3_mod

  implicit none
  
  integer, parameter :: n_ma_max = 100


!==========================================================================
! Parameters for the input file

  integer :: n_corrections_max = 4
  character*200 lattice_file, output_file, comment
  character(200) :: predefined_error_file = ''
  real(rp) :: alignment_multiplier=1.
  real(rp) :: sigma_cutoff=3.
  integer :: seed, seed_used, n_iterations, n_lm_iterations, n_bad_seeds
  integer, allocatable :: seed_list(:)
  character(8) :: seed_write
  logical :: write_orbits=.false.
  real(rp) :: key_value1=0, key_value2=0
  type(ma_struct) :: ma(n_ma_max)
  type(correct_struct) :: correct(10)
  logical :: verbose = .true. ! for verbose writeout during running
  logical :: write_all_lats = .false.
  character(60) :: qtune_regex = ''

  ! Optional ele_noise parameters
!  type(ele_noise_struct) :: ele_noise(n_ele_noise_max) ! n_ele_noise_max defined in ele_noise_mod

  ! job monitoring parameters:
  real(rp) :: max_frac_failed = 0.2 ! fraction of failed jobs
  integer :: min_iter = 20 ! minimum number of iterations to run before determining if a job may have failed
  integer :: bad_corrections = 0 ! initializing total number of failed jobs

  ! "broken" elements:
  logical :: use_veto_list = .false.
  character(200), allocatable :: veto_eles(:)

!==========================================================================
! For simulation of data:
  integer :: n_bpms = 0
  real(rp) :: current = 1000000
  type(measurement_struct), allocatable :: meas(:)
  character(40) :: lat_type = ""
  integer :: obs_points(3) = 0
  character(40) :: obs_point_names(3) = ''

! For TBT tracking:
  integer :: n_turns_orbit = 1024, n_damping_orbit = 0 
  real(rp) :: init_vec(6) = 0., init_vec_orig(6) = 0.
  real(rp), allocatable :: phi(:,:,:), cbar_meas(:,:,:), Amp(:,:,:), on_off_save(:)
  

! Related to phase/coupling simulation
  integer :: n_turns_phase = 32768, n_damping_phase = 1.e5
  real(rp), allocatable :: delta_phi_ax(:),delta_phi_by(:)
  real(rp) :: dphi_ax_avg = 0, dphi_by_avg = 0
  real(rp) :: eta_delta = 1.e-3

! parameters for misalignment; to be read in and loaded into bpm_error_sigmas and ma_params structures
  character(40) :: bpm_mask = "^DET\_[0-9]{2}[ewEW]$" ! default to CESR-type BPM naming
  character(200) :: use_bpm_file = "", bpm_filename ! for masking some BPMs, based on s-location
  type(det_error_struct) :: bpm_error_sigmas
  type(det_struct), allocatable :: bpm(:), bpm_baseline(:)
  type(det_data_struct), allocatable :: bpm_save(:,:)
  real(rp) :: bpm_noise = 0.
  real(rp) :: frac_bpm_failed = 0.

! parameters for tune tracker
  TYPE(tt_param_struct) :: tt_params(max_tt) !max_tt=3 set on tune_tracker module
  INTEGER :: n_TTs = 2     ! number of tune trackers
  INTEGER id(max_tt)
  real(rp) :: target_tunes(1:3) = 0. ! desired tunes for tracking studies, in fractional tune
  real(rp), allocatable :: dk1(:)

! parameters for ac_eta calculation
    character(6) :: eta_method = 'dc' ! options are 'ac_old', 'ac_new', or 'dc'
    logical :: use_new_method = .false. ! default = .false.
    type (cesr_all_data_struct) ac_eta 

!==========================================================================
! other things to be defined:

  type data_struct
     real(rp) emit_x, emit_y, emit_z, rms_x, rms_y, rms_eta_x, rms_eta_y, rms_cbar22, &
          rms_cbar12, rms_cbar11, rms_phi_x, rms_phi_y, rms_param
  end type data_struct

  type(data_struct), allocatable :: datablock(:,:), meas_datablock(:,:), meas_actual_datablock(:,:)

  real(rp) rms_param, start_time, now_time

  type(lat_struct) :: design_ring           ! straight from the lattice file
  type(lat_struct) :: ma_ring, ma_ring_old, ma_ring_init  ! misaligned
  type(ele_struct), pointer :: ele

  type(coord_struct), allocatable :: co(:), co_old(:), orb(:), eta(:)
  type(normal_modes_struct) modes, modes_old
  integer rad_cache, i_ele, i_correction, ix_suffix, time(8), i_iter, i, ix, jx
  integer ix_db, ix_lat, dat_ix, loc, status, corr_status, n_seed
  integer :: tune_write = 1

  character(200) init_file, base_name, out_file, out_meas_file, out_meas_actual_file, opt_file
  logical :: everything_ok = .true.
  real(rp) harvest
  integer :: lun_list(20), lun_data(3) ! for file handling

  ! intelligent bookkeeping toggle:
  logical :: auto_bookkeeper = .false.
  logical err, err1, err2, err3, err4

  ! toggle between ring_ma 'classic' and new sim_meas_mod tracking version:
  character(20) :: meas_type = 'basic' ! alternative is 'tracking'

  ! for PTC:
  type(normal_modes_struct) norm_mode
  real(rp) :: sigma_mat(6,6)
  type(coord_struct) closed_orb
  logical :: use_ptc = .false.

 

!=========================================================================

  namelist /parameters/ lattice_file, predefined_error_file, seed, verbose, &
       n_iterations, output_file, comment, n_lm_iterations, correct, n_corrections_max, &
       write_orbits, key_value1, key_value2, ma, alignment_multiplier, sigma_cutoff, &
       tt_params, n_TTs, init_vec, target_tunes, bpm_mask, bpm_error_sigmas, bpm_noise, &
       n_turns_orbit, n_turns_phase, n_damping_orbit, n_damping_phase, eta_method, &
       max_frac_failed, min_iter, auto_bookkeeper, meas_type, &
       obs_point_names, use_ptc, eta_delta, use_veto_list, write_all_lats, &
       use_bpm_file, frac_bpm_failed, qtune_regex

  namelist /veto_ele_in/ veto_eles


!==========================================================================
! Get init file from command line

  if (command_argument_count() .ne. 1) then
     !write(*,*) "usage: ring_ma <init_file>"
     !stop
     init_file = 'ring_ma2.in'
  else
     call get_command_argument(1, init_file)
  endif
  open(1, file=init_file)
  read(1, nml=parameters)
  close(1)

  ! intelligent bookkeeping toggle:
  bmad_com%auto_bookkeeper = auto_bookkeeper

  write(*,*) " " 
  if (trim(meas_type) .eq. 'tracking') then
     write(*,*) "***Using sim_meas_mod to simulate measurements via tracking"
  elseif(trim(meas_type) .eq. 'basic') then
     write(*,*) "***Using Bmad-calculated values as base for measurement simulations"
  else
     write(*,*) "meas_type ", trim(meas_type), "not known. Stopping here."
     stop
  endif

!  if (write_orbits .eqv. .true.) then
!     write(*,*) "WARNING: write_orbits is not supported at this time."
!     write(*,*) "Set write_orbits = .false. or remove from input file."
!     write(*,*) "Stopping here..."
!     stop
!  endif

  if (use_ptc) then
     write(*,*) "***Using PTC for emittance calculations"
  else
     write(*,*) "***Using radiation_integrals for emittance calculations"
  endif


!==========================================================================
! Transfer the parameters from the input file to the appropriate structures

  dr_misalign_params%alignment_multiplier = alignment_multiplier
  dr_misalign_params%sigma_cutoff         = sigma_cutoff
  correct_ring_params%sigma_cutoff        = sigma_cutoff
  correct_ring_params%n_lm_iterations     = n_lm_iterations

  dr_misalign_params%accumulate_errors    = .false.
  dr_misalign_params%tie_dup_ele          = .true.
  correct_ring_params%eta_delta_e_e       = eta_delta
  correct_ring_params%write_elements      = .true.
  correct_ring_params%skip_dup_ele        = .true.   

  init_vec_orig(:) = init_vec(:)


!==========================================================================
! Read in design ring from lattice and do some intialization


  call bmad_parser(lattice_file, design_ring)

  ! read in list of eles to veto:
  allocate(veto_eles(design_ring%n_ele_max))
  veto_eles(:) = ''
  if (use_veto_list .eqv. .true.) then
	  open(1, file=init_file)
	  read(1, nml=veto_ele_in)
	  close(1)

	  do ix = 1, len(veto_eles) ! loop over input list of eles
	     if (len(trim(veto_eles(ix))) .gt. 0) then ! if user provided a name
	        ele_loop: do jx = 1, design_ring%n_ele_max ! loop over all elements to find ele with that name
	           if (match_reg(trim(design_ring%ele(jx)%name), trim(veto_eles(ix)))) then
	              exit ele_loop
	           elseif (jx .eq. design_ring%n_ele_max) then ! no match found
	              write(*,*) "veto_eles(", ix, ") = ", trim(veto_eles(ix)), " has no match! Stopping here..."
	              stop
	           endif ! if match name w/ ele in lat              
	        enddo ele_loop
	     endif ! if user provided a name
	  enddo ! loop over input list of eles
  endif


  ! set up wigglers:
  do i_ele = 1, design_ring%n_ele_max
     ele => design_ring%ele(i_ele)
     if (ele%key == wiggler$) &
          ele%taylor_map_includes_offsets = .false.
  end do

  ! count the number of correction iterations
  n_corrections_max = 0
  ix_loop: do ix = 1, size(correct)
     do jx = 1, n_detcor_groups
        if (len_trim(correct(ix)%det(jx)%param_name) .gt. 0) then
           n_corrections_max = n_corrections_max + 1
           cycle ix_loop
        endif
     enddo
  enddo ix_loop

  ! initialize corrector calibrations:
  call set_custom_attribute_name('CALIB', err, 2)
  do ix = 1, n_corrections_max
     call init_corrector_calib(design_ring, correct(ix))
  enddo

  ! observation point initialization:
  obs_point_loop: do jx = 1,3
     if (trim(obs_point_names(jx)) .eq. '') cycle
     do ix = 1, design_ring%n_ele_track
        
        if (trim(upcase(obs_point_names(jx))) .eq. trim(upcase(design_ring%ele(ix)%name))) then
           obs_points(jx) = ix
           cycle obs_point_loop
        endif
        
     enddo
  enddo obs_point_loop

  ! get information about indexing for eventual write-out:
  call assign_unique_ele_ids(design_ring)

  rad_cache = 0
  
  call set_on_off(rfcavity$, design_ring, off$)
  call closed_orbit_calc(design_ring, co, 4)
  call lat_make_mat6(design_ring, -1, co)
  call twiss_at_start(design_ring)
  call twiss_propagate_all(design_ring)
  call closed_orbit_calc(design_ring, co_old, 4) ! mostly just to get allocation correct
  call radiation_integrals(design_ring, co, modes, rad_cache)
  call release_rad_int_cache(rad_cache)


  write(*,'(A,2es12.4)') 'Initial normal-mode emittances:', modes%a%emittance, modes%b%emittance !, modes%z%emittance

  if (any(target_tunes .gt. 1.e-12)) then ! only qtune if target_tunes set
     call set_on_off(rfcavity$, design_ring, save_state$, saved_values = on_off_save)
     call set_on_off(rfcavity$, design_ring, off$)
     everything_ok = set_tune3(design_ring, target_tunes, verbose, tune_write, qtune_regex)
     call closed_orbit_calc(design_ring, co, 4)
     call lat_make_mat6(design_ring, -1, co)
     call twiss_at_start(design_ring)
     call twiss_propagate_all(design_ring)
     call set_on_off(rfcavity$, design_ring, restore_state$, saved_values = on_off_save)
  endif

  !--------------------------------------------------------------------
  ! Locate BPMs we wish to study:
  if (correct(1)%bpm_mask .ne. '') then
     write(*,*) ""
     write(*,*) "***********" 
     write(*,*) "***WARNING: correct(1)%bpm_mask is deprecated in favor of bpm_mask."
     write(*,*) "***Will auto-populate bpm_mask = correct(1)%bpm_mask and continue."
     write(*,*) "***********" 
     write(*,*)
     bpm_mask = correct(1)%bpm_mask
  endif
  
  ! find all BPMs within regex match:
  call find_bpms(design_ring, bpm_mask, bpm)
  n_bpms = ubound(bpm,1)
  allocate(bpm_baseline(n_bpms))

  ! toggle BPMs we want off for ALL seeds:
  call toggle_bpms(bpm, design_ring, use_bpm_file)

  if (n_bpms .gt. 0) then
     write(*,'(i0,a)') n_bpms, " BPMs found."
  else ! no BPMs located
     write(*,'(3a)') "No BPMs located with detector mask ", trim(bpm_mask), "! Stopping here..."
     stop
  endif

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



  !------------------------------------------------------
 

  allocate(datablock(n_iterations, 0:n_corrections_max), meas_datablock(n_iterations, 0:n_corrections_max), meas_actual_datablock(n_iterations, 0:n_corrections_max))
  call reallocate_coord(cr_model_co, design_ring%n_ele_track)

  ma_ring_init = design_ring

  ! assume 6x6 tracking. therefore, RF is on
  if ((n_TTs .gt. 0) .and. (trim(meas_type) .eq. 'tracking')) then ! use damping and excitation
     bmad_com%radiation_fluctuations_on = .true. ! tune trackers will dominate well over these effects
     bmad_com%radiation_damping_on = .true.
  else
     bmad_com%radiation_fluctuations_on = .false.
     bmad_com%radiation_damping_on = .false.
  endif

  if (len(trim(predefined_error_file)) .gt. 0) then
     write(*,*) "Reading predefined errors from ", trim(predefined_error_file)
     call bmad_parser2(predefined_error_file, ma_ring_init)
  endif

  ! initialize tune trackers for each correction 
  ! Check if the element for each tune tracker's kicker has kick attributes.

  if (trim(meas_type) .eq. 'tracking') then
     call twiss3_at_start(ma_ring_init,everything_ok) 
     call twiss3_propagate_all(ma_ring_init)
     write(*,*) " "
     call check_tt_kickers(ma_ring_init, tt_params, n_TTs)
     call tt_init(ma_ring_init, co, tt_params, n_TTs, id)
     write(*,*) " " 
  endif

  global_com%exit_on_error = .false.

  !rad_cache = -2

  ! Initialize random-number generator from time or with specific value, for EVERY iteration
  if (seed < 1) then
     seed = 0 ! seed randomizer with CPU time
  end if
  call ran_seed_put(seed)
  call ran_seed_get(seed_used)

  write(*,*) "seed      = ", seed
  write(*,*) "seed_used = ", seed_used

  ! Do a little gymnastics to generate a list of seeds to be used:
  allocate(seed_list(1:n_iterations*5)) ! intentionally allocate more seeds than iterations, in case a seed fails.
  seed_list(1) = seed_used
  do ix = 2, n_iterations*5
     call ran_seed_put(seed_list(ix-1))
     call ran_seed_get(seed_used)
     call ran_gauss(harvest)
     seed_list(ix) = int(harvest * 1.e7 + 1.e8) ! don't want to use abs(); this ensures all positive-definite results
  enddo

!==========================================================================
! Set up file names

  ix_suffix = index(init_file, ".", .true.)
  if (ix_suffix > 0) then
     base_name = init_file(:ix_suffix-1)
  else
     base_name = init_file
  end if
!  if (GRID_mode .eqv. .false.) then
     out_file = trim(base_name) // ".bmad.out"
     out_meas_file = trim(base_name) // ".meas.out"
     out_meas_actual_file = trim(base_name) // ".meas-bmad.out"
!  else ! GRID_mode .eqv. .true.
!     write(seed_write,'(i8.8)') seed_list(1)
!     out_file = trim(base_name) // seed_write // ".bmad.out"
!     out_meas_file = trim(base_name) // seed_write // ".meas.out"
!     out_meas_actual_file = trim(base_name) // seed_write // ".meas-bmad.out"
!  endif

  lun_data(1) = lunget()
  open(lun_data(1),file=out_file, recl=250)
  write(lun_data(1),'("#",A)') trim(comment)

  lun_data(2) = lunget()
  open(lun_data(2),file=out_meas_file, recl=250)
  write(lun_data(2),'("#",A)') trim(comment)

  lun_data(3) = lunget()
  open(lun_data(3),file=out_meas_actual_file, recl=250)
  write(lun_data(3),'("#",A)') trim(comment)


!==========================================================================
! Start iterations
  
  if (write_orbits) then
     opt_file = trim(base_name)//".design"
     write(*,*) "Writing design orbit file: ", trim(opt_file)
     open(2, file=opt_file, recl=250)
     call write_opt(design_ring, co, bpm)
     close(2)
  end if

  n_seed = 0

  iter_loop: do i_iter = 1, n_iterations
     call run_timer('START')

     write(*,*) " "
     write(*,*) " "
     write(*,*) " "
     write(*,*) "============================================"
     write(*,'(a,i0,a,i0)') "STARTING RANDOM SEED ITERATION ", i_iter, "/", n_iterations
     
     !------------------------------------
     ! Generate lattice errors. 
     ! Catch when things are too far off
     n_bad_seeds = 0

     do
        ma_ring = ma_ring_init ! if no predefined errors exist, ma_ring_init = design_ring

        n_seed = n_seed + 1

        ! mask out any BPMs we want to model as "disabled".
        ! Note that all 'is_on' information is stored in 'ring',
        ! therefore it gets reset at the start of every iteration.
        call toggle_bpms(bpm, ma_ring, use_bpm_file, frac_bpm_failed, seed_list(n_seed))
        
        ! reset random number generator for main misalignment
        call ran_seed_put(seed_list(n_seed))
        call ran_seed_get(seed_used)

        ! Misalign lattice; initial measurements
        call dr_misalign(ma_ring, ma)
        co(0)%vec = 0.
        call set_on_off(rfcavity$, ma_ring, save_state$, saved_values = on_off_save)
        call set_on_off(rfcavity$, ma_ring, off$)
        call closed_orbit_calc(ma_ring, co, 4, err_flag = err1)
        call lat_make_mat6(ma_ring, -1, co, err_flag = err2)
        call twiss_at_start(ma_ring, status = status)
        call twiss_propagate_all(ma_ring, err_flag = err3)
        call radiation_integrals(ma_ring, co, modes, rad_cache)
        call release_rad_int_cache(rad_cache)
        call set_on_off(rfcavity$, ma_ring, restore_state$, saved_values = on_off_save)

        if ((modes%a%emittance .lt. 0.) .or. (modes%b%emittance .lt. 0.) .or. (modes%z%emittance .lt. 0.)) err1 = .true. ! set error flag
        if (isnan(modes%a%emittance) .or. isnan(modes%b%emittance) .or. isnan(modes%z%emittance)) err1 = .true. ! set error flag

        ! exit loop if lattice is stable:
        if ((err1 .eqv. .false.) .and. (err2 .eqv. .false.) .and. (err3 .eqv. .false.)  .and. (status == ok$)) exit
        
        ! the following only runs if the lattice is unstable:
        n_bad_seeds = n_bad_seeds + 1

        write(*,'(A,i0,A,i0,a)') "Random seed ", seed_used, " failed. ", n_bad_seeds, " seeds failed so far."
        write(*,*) ""

        if (n_bad_seeds > 5) then
           write(*,'(A)') "Couldn't find any good seeds."
           write(*,'(A)') "Perhaps the misalignment parameters are too severe."
           stop
        end if

     end do

     write(*,'(a,i0)') "Random seed used: ", seed_used

     if (write_all_lats) then
        call write_misalignments(ma_ring, ma, seed_used)
     else
        call write_misalignments(ma_ring, ma, -1)
     endif

     !-----------------
     ! Generate BPM errors and write to file
     lun_list(1) = lunget()
     
     if (write_all_lats) then
        write(bpm_filename,'(a,i9.9,a)') 'bpms.', seed_used, '.ma'
     else
        bpm_filename = 'bpms.ma'
     endif

     open(unit=lun_list(1), file=bpm_filename, status='replace')
     write(lun_list(1), '(a7,12a14)') "!index", 'is_on', "x-offset", "y-offset", "tilt", "g1", "g2", "g3", "g4", "shear_x", "timing"
     write(lun_list(1), '(a7, 11a14)') ""
     
     
     do i=1,n_bpms
        call bpm_errors(bpm(i), bpm_error_sigmas, ma_ring%ele(bpm(i)%ix_quad))
        write(lun_list(1), '(i7,L14,12es14.4)') bpm(i)%ix_db, ma_ring%ele(bpm(i)%ix_lat)%is_on, &
             bpm(i)%x_offset, bpm(i)%y_offset, bpm(i)%tilt, bpm(i)%gain(:), bpm(i)%shear_x, &
             bpm(i)%timing(:)
     enddo
     close(unit=lun_list(1))
     
     !-------------------

     i_correction = 0

     if (any(target_tunes .gt. 1.e-12)) then
        call set_on_off(rfcavity$, ma_ring, save_state$, saved_values = on_off_save)
        call set_on_off(rfcavity$, ma_ring, off$)
        everything_ok = set_tune3(ma_ring, target_tunes, verbose, tune_write, qtune_regex)
        call closed_orbit_calc(ma_ring, co, 4)
        call lat_make_mat6(ma_ring, -1, co)
        call twiss_at_start(ma_ring)
        call twiss_propagate_all(ma_ring)
        call radiation_integrals(ma_ring, co, modes, rad_cache)
        call release_rad_int_cache(rad_cache)
        call set_on_off(rfcavity$, ma_ring, restore_state$, saved_values = on_off_save)
     endif
     !call twiss3_at_start(ma_ring, everything_ok)
     !call twiss3_propagate_all(ma_ring)

     if (write_orbits) then
        write(opt_file, '(A,".opt.",i3.3)') trim(base_name), i_iter
        write(*,*) "Writing optimization file: ", trim(opt_file)
        open(2, file=opt_file, recl=250)
        call write_opt(ma_ring, co, bpm)
     end if     
     
     ! Store initial values
     call ring_to_data(ma_ring, design_ring, co, modes, datablock(i_iter, 0), bpm, 'bmad      ', obs_points)
     call write_data(datablock, lun_data(1), i_iter, 0, seed_used = seed_used)

     ! Apply correction(s)
     corr_loop: do i_correction = 1, n_corrections_max

        ! for outputting statistics on measured values
        ma_ring_old = ma_ring 
        modes_old = modes
        co_old = co

        ! make sure the first corrector is populated, otherwise cycle
        if ((len(trim(correct(i_correction)%cor(1)%attrib_name)) .eq. 0)) cycle corr_loop
        
        call run_timer ('READ', start_time)

        write(*,*) " " 
        if (i_correction .gt. 1) write(*,*) "============================================"
        write(*,'(a,i0,a,i0)') "Seed iteration #", i_iter, "   |   Correction iteration #", i_correction
        write(*,*) " "
               
        write(*,'(A,2es12.4)') 'Normal-mode emittances:', modes%a%emittance, modes%b%emittance !, modes%z%emittance


        if (.not. any(init_vec_orig(:) > 1.e-12)) then ! floating point comparison; assume < 1.e-12 is negligible
           init_vec(:) = co(0)%vec(:)
        else
           init_vec(:) = init_vec_orig(:)
        endif

        rms_param = 0.
        cr_model_ring = design_ring

        !!!=============
        !!! correct_ring call:
        call correct_ring(ma_ring, tt_params, id, bpm, correct(i_correction), rms_param, n_turns_phase, n_turns_orbit, n_damping_phase, &
             n_damping_orbit, init_vec, eta_method, verbose, meas_type, meas, seed_used, write_all_lats, corr_status)         
        !!!=============
        !!!=============

        call set_on_off(rfcavity$, ma_ring, save_state$, saved_values = on_off_save)
        call set_on_off(rfcavity$, ma_ring, off$)
        call closed_orbit_calc(ma_ring, co, 4, err_flag = err1)
        call lat_make_mat6(ma_ring, -1, co, err_flag = err2)
        call twiss_at_start(ma_ring, status = status)
        call twiss_propagate_all(ma_ring, err_flag = err4)
        call radiation_integrals(ma_ring, co, modes, rad_cache)
        call release_rad_int_cache(rad_cache)
        call set_on_off(rfcavity$, ma_ring, restore_state$, saved_values = on_off_save)

        if (err1 .or. err2 .or. status /= ok$ .or. err4 .or. corr_status /= 0) then 
           bad_corrections = bad_corrections + 1
           if ( ( (bad_corrections*1.0 / i_iter) .gt. max_frac_failed) .and. (i_iter .gt. min_iter) ) then 
              ! too many failed iterations, cancel the job
              write(*,'(a,i0,a,i0,a)') "***WARNING: ", bad_corrections, " jobs failed out of ", i_iter, &
                   " seeds tried. Perhaps misalignments are too severe. Canceling job here..."
              stop
           else 
              write(*,'(a,i0,a,i0,a)') "***Correction failed for seed iteration ", i_iter, ", correction ", i_correction, "***"
              write(*,'(a,i0,a,i0,a)') "Failure rate for this job: ", bad_corrections, "/", i_iter, ". Cycling iter_loop..."
              cycle iter_loop
           endif
        endif ! everything_ok

        ! q_tune to user-specified tunes EVERY ITERATION - important for tune trackers!
        if (any(target_tunes .gt. 1.e-12)) then
           call set_on_off(rfcavity$, ma_ring, save_state$, saved_values = on_off_save)
           call set_on_off(rfcavity$, ma_ring, off$)
           everything_ok = set_tune3(ma_ring, target_tunes, verbose, tune_write, qtune_regex)
           call closed_orbit_calc(ma_ring, co, 4, err_flag = err1)
           call lat_make_mat6(ma_ring, -1, co, err_flag = err2)
           call twiss_at_start(ma_ring, status = status)
           call twiss_propagate_all(ma_ring, err_flag = err4)
           call radiation_integrals(ma_ring, co, modes, rad_cache)
           call release_rad_int_cache(rad_cache)
           call set_on_off(rfcavity$, ma_ring, restore_state$, saved_values = on_off_save)
        endif

        if (err1 .or. err2 .or. status /= ok$ .or. err4) then 
           bad_corrections = bad_corrections + 1
           if ( ( (bad_corrections*1.0 / i_iter) .gt. max_frac_failed) .and. (i_iter .gt. min_iter) ) then 
              ! too many failed iterations, cancel the job
              write(*,'(a,i0,a,i0,a)') "***WARNING: ", bad_corrections, " jobs failed out of ", i_iter, &
                   " seeds tried. Perhaps misalignments are too severe. Canceling job here..."
              stop
           else 
              write(*,'(a,i0,a,i0,a)') "***Correction failed for seed iteration ", i_iter, ", correction ", i_correction, "***"
              write(*,'(a,i0,a,i0,a)') "Failure rate for this job: ", bad_corrections, "/", i_iter, ". Cycling iter_loop..."
              cycle iter_loop
           endif
        endif ! everything_ok

! redundant:
!!        rad_cache = -2
!        call set_on_off(rfcavity$, ma_ring, save_state$, saved_values = on_off_save)
!        call set_on_off(rfcavity$, ma_ring, on$)
!        call closed_orbit_calc(ma_ring, co, 6)
!        call lat_make_mat6(ma_ring, -1, co)
!        call twiss_at_start(ma_ring)
!        call twiss_propagate_all(ma_ring)
!        call radiation_integrals(ma_ring, co, modes, rad_cache)
!        call release_rad_int_cache(rad_cache)
!        call set_on_off(rfcavity$, ma_ring, restore_state$, saved_values = on_off_save)


        if (write_orbits) then
           !write(opt_file, '(A,".opt.",i3.3)') trim(base_name), i_iter
           !write(*,*) "Writing optimization file: ", trim(opt_file)
           !open(2, file=opt_file, recl=250)
           call write_opt(ma_ring, co, bpm)     
           !close(2)
        endif

        ! for exact ring parameters:
        call ring_to_data(ma_ring, design_ring, co, modes, datablock(i_iter, i_correction), bpm, 'bmad      ', obs_points)
        call write_data(datablock, lun_data(1), i_iter, i_correction)

        ! for measured data: need to write previous iteration's results now, after measurements have been done
        ! at the start of the current iteration
        call ring_to_data(ma_ring_old, design_ring, co_old, modes_old, meas_datablock(i_iter, i_correction-1), bpm, &
             'meas      ', obs_points, meas)
        call ring_to_data(ma_ring_old, design_ring, co_old, modes_old, meas_actual_datablock(i_iter, i_correction-1), bpm, &
             'meas-bmad ',obs_points,meas)
        call write_data(meas_datablock, lun_data(2), i_iter, i_correction-1)
        call write_data(meas_actual_datablock, lun_data(3), i_iter, i_correction-1)

        call run_timer ('READ', now_time)
        print '(a, f10.0)', 'Correction loop time (sec):', now_time - start_time

     !!!===========================
     !!! END OF CORR_LOOP
     !!!===========================

     end do corr_loop ! correction iterations
     
     
     !===============================================================================
     ! simulate orbit, phase/coupling, ac_eta measurements AFTER loading ALL corrections.

     write(*,*) " "
     write(*,*) "*** Corrections complete. Making final measurements... ***"
     write(*,'(A,2es12.4)') 'Normal-mode emittances:', modes%a%emittance, modes%b%emittance !, modes%z%emittance
     write(*,*) " "

     if (.not. allocated(meas)) allocate(meas(1:size(bpm)))

    if (trim(meas_type) .eq. 'tracking') then
       call set_on_off(rfcavity$, ma_ring, save_state$, saved_values = on_off_save)
       call set_on_off(rfcavity$, ma_ring, on$)
       call twiss3_at_start(ma_ring, err) 
       call twiss3_propagate_all(ma_ring)
       call set_on_off(rfcavity$, ma_ring, restore_state$, saved_values = on_off_save)
       ! init_vec has already been checked in main ring_ma2 routine
       
       ! reset each tune tracker now, after computing closed orbit
       ! for this specific correction on this misaligned lattice.
       
       tt_params(1)%offset = co(tt_params(1)%bpm_loc)%vec(1)
       tt_params(2)%offset = co(tt_params(2)%bpm_loc)%vec(3)
       !    tt_params(3)%offset = co(tt_params(3)%bpm_loc)%vec(5) ! don't change longitudinal; seems to cause problems
       do ix=1, n_TTs 
          call reset_dTT(id(ix), tt_params(ix))
       enddo
       
       ! simulate the data:
       write(*,*) "Simulating orbit measurement..."
       call sim_tbt_data(ma_ring, bpm, n_turns_orbit, .true., 0, tt_params, id, n_damping_orbit, &
            init_vec, bpm_save) ! disable tunetrackers for orbit measurement
       write(*,*) " "
       if (trim(eta_method) .ne. 'dc') then
          write(*,*) "Simulating phase/cbar and ac_eta measurements after correction..."
          call sim_ac_and_cbar_data(ma_ring, bpm, n_turns_phase, 3, tt_params, id, n_damping_phase, &
               init_vec, use_new_method, ac_eta, phi, cbar_meas, Amp)
       else ! use DC dispersion
          write(*,*) "Simulating phase/cbar measurements after correction..."
          call sim_cbar_data(ma_ring, bpm, n_turns_phase, .true., 2, tt_params, id, n_damping_phase, &
               init_vec, phi, cbar_meas, Amp)
          write(*,*) ""
          write(*,*) "Simulating DC dispersion after correction..."
          call sim_dc_eta_data(ma_ring, bpm, n_turns_orbit, &
               correct_ring_params%eta_delta_e_e, eta)
       endif
       !call set_on_off(rfcavity$, ma_ring, off$)

       !do ix=1, n_TTs 
       !   call reset_dTT(id(ix), tt_params(ix))
       !enddo
       
    else ! must be meas_type = 'basic'; already filtered out bad options at startup

       call sim_eta_data(ma_ring, bpm, correct_ring_params%eta_delta_e_e, eta)

       if (.not. allocated(cbar_meas)) allocate(cbar_meas(1:n_bpms,2,2))

       ! orbit and cbar:
       bpm(:)%vec(1) = co(bpm(:)%ix_lat)%vec(1)
       bpm(:)%vec(3) = co(bpm(:)%ix_lat)%vec(3)
       do jx = 1, n_bpms

          ele => ma_ring%ele(bpm(jx)%ix_lat)

          call apply_bpm_errors(bpm(jx), current, lat_type)  ! note: also transfers vec(6)-->amp(4) button signals
          
          ! transfer back to co for now:
          co(bpm(jx)%ix_lat)%vec(:) = bpm(jx)%vec

          ! note: this also populates the cbar_meas struct-- don't comment out!
          call rotate_cbar(ele, bpm(jx)%tilt, cbar_meas(jx,:,:))

       enddo

    endif ! meas_type

     
     if (.not. allocated(orb)) allocate(orb(size(co)))
     
     ! load data into relevant data structures:
     do ix = 1, size(bpm)
        ele => ma_ring%ele(bpm(ix)%ix_lat)
        loc = bpm(ix)%ix_lat
        
        if (trim(meas_type) .eq. 'tracking') then
           
           ! new values:
           meas(ix)%orbit_x = sum(bpm_save(:,ix)%vec(1))/size(bpm_save(:,ix)) ! average orbit over several turns
           meas(ix)%orbit_y = sum(bpm_save(:,ix)%vec(3))/size(bpm_save(:,ix))
           if (trim(eta_method) .ne. 'dc') then
              meas(ix)%eta_x   = ac_eta%ac_eta_x(ix)%value
              meas(ix)%eta_y   = ac_eta%ac_eta_y(ix)%value
           else
              meas(ix)%eta_x = eta(loc)%vec(1)
              meas(ix)%eta_y = eta(loc)%vec(3)
           endif
           meas(ix)%phi_a   = phi(ix,1,1)  ! already phase-offset-corrected
           meas(ix)%phi_b   = phi(ix,2,2)  
           meas(ix)%cbar22  = cbar_meas(ix,2,2)
           meas(ix)%cbar12  = cbar_meas(ix,1,2)
           meas(ix)%cbar11  = cbar_meas(ix,1,1)

           ! turn off tune trackers, otherwise emittance will 
           ! be reported wrong
           do jx=1, n_TTs 
              call reset_dTT(id(jx), tt_params(jx))
           enddo
           
        else ! must be 'basic':
           
           meas(ix)%orbit_x = co(loc)%vec(1)
           meas(ix)%orbit_y = co(loc)%vec(3)
           meas(ix)%eta_x   = eta(loc)%vec(1)
           meas(ix)%eta_y   = eta(loc)%vec(3)
           meas(ix)%phi_a = ele%a%phi
           meas(ix)%phi_b = ele%b%phi
           
           meas(ix)%cbar22  = cbar_meas(ix,2,2)
           meas(ix)%cbar12  = cbar_meas(ix,1,2)
           meas(ix)%cbar11  = cbar_meas(ix,1,1)
           
        endif
     enddo
     
     call ring_to_data(ma_ring, design_ring, co, modes, meas_datablock(i_iter, n_corrections_max), bpm, 'meas      ', obs_points, meas)
     call write_data(meas_datablock, lun_data(2), i_iter, n_corrections_max)
     call ring_to_data(ma_ring, design_ring, co, modes, meas_actual_datablock(i_iter, n_corrections_max), bpm, 'meas-bmad ', obs_points, meas)
     call write_data(meas_actual_datablock, lun_data(3), i_iter, n_corrections_max)
     
     if (write_orbits) close(2)
     
     write(*,'(a,i0,a,i0)') "Failure rate for this job: ", bad_corrections, "/", i_iter


     call run_timer ('READ', now_time)
     print '(a, f10.0)', 'Iteration loop time (sec):', now_time

  !!!===================
  !!! END OF ITERATION LOOP
  !!!===================
  end do iter_loop
  
  ! destruct dTTs:
  if (trim(meas_type) .eq. 'tracking') then
     do ix = 1, n_TTs
        call dest_dTT(id(ix))           
     enddo
  endif

  ! Write summary
  call write_data(datablock, lun_data(1))
  call write_data(meas_datablock, lun_data(2))
  call write_data(meas_actual_datablock, lun_data(3))
  close(lun_data(1))
  close(lun_data(2))
  close(lun_data(3))


contains
!==========================================================================
! Routine for writing out the complete element-by-element
! parameters for every seed. Can generate a LOT of output.

  subroutine write_opt(ring, co, bpm)
    implicit none
    type(lat_struct) :: ring
    type(det_struct) :: bpm(:)
    type(coord_struct) :: co(:)
    integer i_ele, i
    real(rp) cbar(2,2)

    write(2,'(2A5,9A13,A18)') '# cor', 'ele', 's', 'x', 'y', 'eta_x', 'eta_y', &
         'phi_a', 'phi_b', 'cbar12', 'length', 'name'
    do i = 1, size(bpm)
       i_ele = bpm(i)%ix_lat ! only print out BPMs, as data only EXISTS for BPMs.
       call c_to_cbar(ring%ele(i_ele), cbar)
       write(2,'(2i5,9es13.4,a18)') i_correction, i_ele, ring%ele(i_ele)%s, &
            co(i_ele)%vec(1), co(i_ele)%vec(3), &
            ring%ele(i_ele)%x%eta, ring%ele(i_ele)%y%eta, &
            ring%ele(i_ele)%a%phi - design_ring%ele(i_ele)%a%phi, &
            ring%ele(i_ele)%b%phi - design_ring%ele(i_ele)%b%phi, &
            cbar(1,2), value_of_attribute(ring%ele(i_ele),'L',err1), trim(ring%ele(i_ele)%name)
    end do
    write(2,*)
    write(2,*)
  end subroutine write_opt

!==========================================================================
! Routine to transfer relevant parameters from a ring to the DATABLOCK

  subroutine ring_to_data(ring, design_ring, co, modes, data, bpm, data_type, obs_points, meas)

    implicit none

    type(lat_struct) :: ring, design_ring
    type(coord_struct), intent(in) :: co(:)
    type(det_struct), intent(in) :: bpm(:)
    type(data_struct), intent(out) :: data
    type(measurement_struct), optional :: meas(:)
    real(rp) l, ltot, cbar(2,2)
    integer, intent(in) :: obs_points(:)
    integer i_ele, ix, i, j, k
    character(10), intent(in) :: data_type

    integer :: rad_cache = 0

    ! for PTC:
    real(rp) emit_norm(3), emit_projected(3)
    type (normal_modes_struct), intent(inout) :: modes
    real(rp) sigma_mat(6,6), sigma_beta_mat(6,6)
    type (coord_struct) closed_orb

    ! ------------------------------------------
    ! ------------------------------------------
    ! Emittances: compute using PTC:
    
    !-----------------------------------------
    !data%emit_x = modes%a%emittance
    !data%emit_y = modes%b%emittance
    !data%emit_z = 0.
   
    !------------------------------------------
!    if (use_ptc) then
!      if (trim(data_type) .eq. 'bmad') then ! only need to recalculate this once
!        call set_on_off(rfcavity$, ring, save_state$)
!        call set_on_off(rfcavity$, ring, on$)
!        call lat_to_ptc_layout(ring)
!        
!        emit_norm(:) = 0.
!        emit_projected(:) = 0.
!   
!        ! Only allow one observation point for now, due to computation time for PTC 
!        !do ix = 1,3
!    
!           ! compute emittance at desired source point:    
!           call ptc_emit_calc(ring%ele(obs_points(1)), modes, sigma_mat, closed_orb)
!           !call ptc_emit_calc(ring%ele(obs_points(ix)), modes, sigma_mat, closed_orb)
!        endif
!     else ! don't use PTC:
!        rad_cache = 0
!        call radiation_integrals(ring, co, modes, rad_cache)
!        call release_rad_int_cache(rad_cache)
!     endif

     ! populate emit_norm for both PTC and radiation_integrals:
     !if (ix .eq. 1) then 
     emit_norm(1) = modes%a%emittance
     !elseif (ix .eq. 3) then 
     emit_norm(2) = modes%b%emittance
     !elseif (ix .eq. 3) then 
     emit_norm(3) = modes%z%emittance
     !endif
     
!     if (use_ptc) then    
!        if (trim(data_type) .eq. 'bmad') then 
!           call kill_ptc_layouts(ring)
!           call set_on_off(rfcavity$, ring, restore_state$)
!        endif
!  
!      ! next block derived from bsim_cesr/code/make_sigma_mat.f90:
!      DO j=1,6
!         DO k=j,6        
!            sigma_beta_mat(j,k) = sigma_mat(j,k) - sigma_mat(j,6)*sigma_mat(k,6)/sigma_mat(6,6)
!            IF( j .ne. k ) THEN
!               sigma_beta_mat(k,j) = sigma_beta_mat(j,k)
!            ENDIF
!         ENDDO
!      ENDDO
!     
!      emit_projected(1) = sqrt(sigma_beta_mat(1,1)*sigma_beta_mat(2,2) - sigma_beta_mat(1,2)**2)
!      emit_projected(2) = sqrt(sigma_beta_mat(3,3)*sigma_beta_mat(4,4) - sigma_beta_mat(3,4)**2)
!      emit_projected(3) = emit_norm(3) ! for now
!
!    else ! don't use PTC:
      emit_projected(:) = emit_norm(:) ! temporary
!    endif
   
!    !enddo ! for multiple observation points
    
    ! load into necessary structures:
    if (trim(data_type) .eq. 'bmad') then
       data%emit_x = emit_norm(1)
       data%emit_y = emit_norm(2)
       !data%emit_z = emit_norm(3)
       data%emit_z = 0. ! temporary
    elseif (trim(data_type) .eq. 'meas') then
       data%emit_x = emit_projected(1)
       data%emit_y = emit_projected(2)
       !data%emit_z = emit_projected(3)
       data%emit_z = 0. ! temporary
    elseif (trim(data_type) .eq. 'meas-bmad') then
       data%emit_x = emit_projected(1) - emit_norm(1)
       data%emit_y = emit_projected(2) - emit_norm(2)
       !data%emit_z = emit_projected(3) - emit_norm(3)
       data%emit_z = 0. ! temporary
    endif


    ! Initialize RMS values
    data%rms_x = 0.
    data%rms_y = 0.
    data%rms_eta_x = 0.
    data%rms_eta_y = 0.
    data%rms_cbar22 = 0.
    data%rms_cbar12 = 0.
    data%rms_cbar11 = 0.
    data%rms_phi_x = 0.
    data%rms_phi_y = 0.

    ltot = 0.


    if (trim(data_type) .eq. 'meas-bmad') then
       
       do i = 1, size(bpm)
          i_ele = bpm(i)%ix_lat
          ele => ring%ele(i_ele)
          call c_to_cbar(ele, cbar)
          
          !! Take a simple average over BPMs
          data%rms_x      = data%rms_x      + (meas(i)%orbit_x - co(i_ele)%vec(1))**2
          data%rms_y      = data%rms_y      + (meas(i)%orbit_y - co(i_ele)%vec(3))**2
          data%rms_eta_x  = data%rms_eta_x  + (meas(i)%eta_x  - ring%ele(i_ele)%x%eta)**2 
          data%rms_eta_y  = data%rms_eta_y  + (meas(i)%eta_y  - ring%ele(i_ele)%y%eta)**2 
          data%rms_cbar22 = data%rms_cbar22 + (meas(i)%cbar22 - cbar(2,2))**2
          data%rms_cbar12 = data%rms_cbar12 + (meas(i)%cbar12 - cbar(1,2))**2
          data%rms_cbar11 = data%rms_cbar11 + (meas(i)%cbar11 - cbar(1,1))**2
          data%rms_phi_x  = data%rms_phi_x  + (meas(i)%phi_a - ring%ele(i_ele)%a%phi)**2
          data%rms_phi_y  = data%rms_phi_y  + (meas(i)%phi_b - ring%ele(i_ele)%b%phi)**2
       end do

       data%rms_x      = sqrt(data%rms_x / size(bpm))
       data%rms_y      = sqrt(data%rms_y / size(bpm))
       data%rms_eta_x  = sqrt(data%rms_eta_x / size(bpm))
       data%rms_eta_y  = sqrt(data%rms_eta_y / size(bpm))
       data%rms_cbar22 = sqrt(data%rms_cbar22 / size(bpm))
       data%rms_cbar12 = sqrt(data%rms_cbar12 / size(bpm))
       data%rms_cbar11 = sqrt(data%rms_cbar11 / size(bpm))
       data%rms_phi_x  = sqrt(data%rms_phi_x / size(bpm))
       data%rms_phi_y  = sqrt(data%rms_phi_y / size(bpm))
       data%rms_param  = 0. ! not relevant for this output
       
    elseif (trim(data_type) .eq. 'meas') then

       ! foo - 
       !do i = 1, size(bpm)
       !   i_ele = bpm(i)%ix_lat
       !   write(17,'(i,3es18.8)') i, meas(i)%eta_x, ring%ele(i_ele)%x%eta, design_ring%ele(i_ele)%x%eta
       !enddo
       do i = 1, size(bpm)
          i_ele = bpm(i)%ix_lat
          ele => ring%ele(i_ele)
          call c_to_cbar(ele, cbar)

          !! Take a simple average over BPMs
          data%rms_x      = data%rms_x      + (meas(i)%orbit_x)**2
          data%rms_y      = data%rms_y      + (meas(i)%orbit_y)**2
          data%rms_eta_x  = data%rms_eta_x  + (meas(i)%eta_x - design_ring%ele(i_ele)%x%eta)**2
          data%rms_eta_y  = data%rms_eta_y  + (meas(i)%eta_y - design_ring%ele(i_ele)%y%eta)**2
          data%rms_cbar22 = data%rms_cbar22 + (meas(i)%cbar22)**2
          data%rms_cbar12 = data%rms_cbar12 + (meas(i)%cbar12)**2
          data%rms_cbar11 = data%rms_cbar11 + (meas(i)%cbar11)**2
          data%rms_phi_x  = data%rms_phi_x  + (mod(meas(i)%phi_a - design_ring%ele(i_ele)%a%phi, twopi))**2
          data%rms_phi_y  = data%rms_phi_y  + (mod(meas(i)%phi_b - design_ring%ele(i_ele)%b%phi, twopi))**2
       end do

       data%rms_x      = sqrt(data%rms_x / size(bpm))
       data%rms_y      = sqrt(data%rms_y / size(bpm))
       data%rms_eta_x  = sqrt(data%rms_eta_x / size(bpm))
       data%rms_eta_y  = sqrt(data%rms_eta_y / size(bpm))
       data%rms_cbar22 = sqrt(data%rms_cbar22 / size(bpm))
       data%rms_cbar12 = sqrt(data%rms_cbar12 / size(bpm))
       data%rms_cbar11 = sqrt(data%rms_cbar11 / size(bpm))
       data%rms_phi_x  = sqrt(data%rms_phi_x / size(bpm))
       data%rms_phi_y  = sqrt(data%rms_phi_y / size(bpm))
       data%rms_param  = 0. ! not relevant for this output 


    elseif (trim(data_type) .eq. 'bmad') then ! must be data_type .eq. 'bmad'
       
       do i = 1, size(bpm)
          i_ele = bpm(i)%ix_lat
          ele => ring%ele(i_ele)
          call c_to_cbar(ele, cbar)
          
          !! Take a simple average over BPMs
          data%rms_x      = data%rms_x      + co(i_ele)%vec(1)**2
          data%rms_y      = data%rms_y      + co(i_ele)%vec(3)**2
          data%rms_eta_x  = data%rms_eta_x  + (ring%ele(i_ele)%x%eta - design_ring%ele(i_ele)%x%eta)**2 
          data%rms_eta_y  = data%rms_eta_y  + (ring%ele(i_ele)%y%eta - design_ring%ele(i_ele)%y%eta)**2 
          data%rms_cbar22 = data%rms_cbar22 + cbar(2,2)**2
          data%rms_cbar12 = data%rms_cbar12 + cbar(1,2)**2
          data%rms_cbar11 = data%rms_cbar11 + cbar(1,1)**2
          data%rms_phi_x  = data%rms_phi_x  + (mod(ring%ele(i_ele)%a%phi - design_ring%ele(i_ele)%a%phi, twopi))**2
          data%rms_phi_y  = data%rms_phi_y  + (mod(ring%ele(i_ele)%b%phi - design_ring%ele(i_ele)%b%phi, twopi))**2
       end do
       data%rms_x      = sqrt(data%rms_x / size(bpm))
       data%rms_y      = sqrt(data%rms_y / size(bpm))
       data%rms_eta_x  = sqrt(data%rms_eta_x / size(bpm))
       data%rms_eta_y  = sqrt(data%rms_eta_y / size(bpm))
       data%rms_cbar22 = sqrt(data%rms_cbar22 / size(bpm))
       data%rms_cbar12 = sqrt(data%rms_cbar12 / size(bpm))
       data%rms_cbar11 = sqrt(data%rms_cbar11 / size(bpm))
       data%rms_phi_x  = sqrt(data%rms_phi_x / size(bpm))
       data%rms_phi_y  = sqrt(data%rms_phi_y / size(bpm))
       data%rms_param  = rms_param

    else ! data_type not recognized
       write(*,*) "ring_to_data - data type not supported."
       stop
       
    endif

  end subroutine ring_to_data





!==========================================================================
! Routine to write out the results of an individual seed, OR to write
! a summary table

  subroutine write_data(data, lun_num, iter, cor, seed_used)
    implicit none
    type(data_struct), allocatable, intent(in) :: data(:,:)
    integer, intent(in), optional :: iter, cor, seed_used
    integer, intent(in) :: lun_num
    integer i_cor

    if (present(iter)) then
       if (iter==1 .and. cor==0) then
          write(lun_num,'("#",2a6,12a14)') "iter", "cor", "emit_a", "emit_b", "emit_z", "rms_x", "rms_y", &
               "rms_eta_x", "rms_eta_y", "rms_cbar22", "rms_cbar12", "rms_cbar11", "rms_phi_x", "rms_phi_y" !, &
               !"param_rms" , "key_val1", "key_val2"
       endif

       if ((cor == 0) .and. present(seed_used)) write(lun_num,'(a,i0)') "# Random seed: ", seed_used

       write(lun_num,'(" ",2i6,12es14.5)') iter, cor, &
            data(iter, cor)%emit_x, &
            data(iter, cor)%emit_y, &
            data(iter, cor)%emit_z, &
            data(iter, cor)%rms_x, &
            data(iter, cor)%rms_y, &
            data(iter, cor)%rms_eta_x, &
            data(iter, cor)%rms_eta_y, &
            data(iter, cor)%rms_cbar22, &
            data(iter, cor)%rms_cbar12, &
            data(iter, cor)%rms_cbar11, &
            data(iter, cor)%rms_phi_x, &
            data(iter, cor)%rms_phi_y
            !data(iter, cor)%rms_param, &
            !key_value1, key_value2
    else

       write(lun_num,*)
       write(lun_num,'(a)') "# Summary"
       write(lun_num, '("#",a6,a10,6a14)') "cor", "param", "key_val1", "key_val2", "mean", "sigma", "50pct", "95pct"
       do i_cor = 0, n_corrections_max
          if (i_cor > 0) then
             if (all(len_trim(correct(i_cor)%cor(:)%attrib_name) .eq. 0)) cycle
          end if
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "emit_a",  key_value1, key_value2, data_line(data(:,i_cor)%emit_x)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "emit_b",  key_value1, key_value2, data_line(data(:,i_cor)%emit_y)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "emit_z",  key_value1, key_value2, data_line(data(:,i_cor)%emit_z)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "orbit_x", key_value1, key_value2, data_line(data(:,i_cor)%rms_x)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "orbit_y", key_value1, key_value2, data_line(data(:,i_cor)%rms_y)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "eta_x",   key_value1, key_value2, data_line(data(:,i_cor)%rms_eta_x)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "eta_y",   key_value1, key_value2, data_line(data(:,i_cor)%rms_eta_y)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "cbar22",  key_value1, key_value2, data_line(data(:,i_cor)%rms_cbar22)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "cbar12",  key_value1, key_value2, data_line(data(:,i_cor)%rms_cbar12)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "cbar11",  key_value1, key_value2, data_line(data(:,i_cor)%rms_cbar11)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "phi_x",   key_value1, key_value2, data_line(data(:,i_cor)%rms_phi_x)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "phi_y",   key_value1, key_value2, data_line(data(:,i_cor)%rms_phi_y)
          write(lun_num, '(" ",i6,a10,6es14.5)') i_cor, "param",   key_value1, key_value2, data_line(data(:,i_cor)%rms_param)
       end do
    end if
  end subroutine write_data
  
!==========================================================================
! Routine to compute statistical summary information for an array

  function data_line(array)
    implicit none

    real(rp) :: data_line(4)
    real(rp), intent(in) :: array(:)
    integer indx(size(array)), k50, k95

    ! Compute mean and standard deviation
    if (size(array) > 1) then
       call avevar(array, data_line(1), data_line(2))
    else ! only one data point; not enough to do statistics on
       data_line(1) = array(1)
       data_line(2) = 0
    endif

    data_line(2) = sqrt(data_line(2))

    ! Compute quantiles if we have enough seeds
    if (size(array) > 10) then
       call indexer(array, indx)
       k50 = ceiling(.50 * size(array))
       k95 = ceiling(.95 * size(array))
       data_line(3) = array(indx(k50))
       data_line(4) = array(indx(k95))
    else
       data_line(3:4) = -1
    end if
  end function data_line

end program ring_ma2