subroutine init_lattice (u, read_lat_file, make_groups)

use cesrv_init_groups               
use taylor2_mod
use bookkeeper_mod
use z_tune_mod
implicit none

type (universe_struct), target :: u


integer i, j, ix, ixc, ie, d_counts(5), status
integer make_groups, ix_last, ix_ele
integer, external :: vnumbr

real(rp) beta_a, beta_b, s_last
real(rp), save, target :: dummy = 0
real(rp) dE_temp

character(200) bmad_file, dir, base_name

logical read_lat_file, energy_var_temp, err

!------------------------------------------------------------------
! default var settings

u%var(:)%exists = .false.
u%var(:)%good_var = .false.
u%var(:)%good_user = .true.
u%var(:)%dvar_dcu = 0.0
u%var(:)%do_limit_calc = .true.
u%var(:)%makeup_method = make_mat6$
u%var(:)%high_target_lim =  1e20
u%var(:)%low_target_lim = -1e20
u%var(:)%ix_ele = 0
u%var(:)%dvar_dcu = 0

do i = lbound(u%var, 1), ubound(u%var, 1)
  u%var(i)%model => dummy    ! just to initially point somewhere
enddo

! init lattice

print *, 'Initializing Lattice and Data'

if (read_lat_file) then
  call lattice_to_bmad_file_name (logic%lattice, bmad_file)
  call set_ptc (taylor_order = 3) ! reset to standard so bmad_parser will no complain
  call bmad_parser(bmad_file, u%ring, err_flag = err)
  if (err) stop
  print *, 'Lattice file: ', trim(bmad_file)
endif

! Initialize orb

if (allocated(u%orb)) then
  deallocate(u%orb)
endif
call reallocate_coord (u%orb, u%ring%n_ele_max)  

! turn off the rf and calculate the design twiss parameters.

call set_on_off (rfcavity$, u%ring, off$)

if (u%ring%lattice /= " ") then
  logic%lattice = u%ring%lattice
  call upcase_string(logic%lattice)
endif
logic%ring_initialized = .true.
call twiss_at_start (u%ring)
call twiss_propagate_all (u%ring)

! check that the input file name is consistant with the lattice name

bmad_file = u%ring%input_file_name
call str_upcase (bmad_file, bmad_file)
ix = splitfilename (bmad_file, dir, base_name)
ix = index(base_name, 'BMAD_')
if (ix /= 0) base_name = base_name(ix+5:)
ix = index(base_name, '.LAT')
if (ix /= 0) base_name = base_name(:ix-1)

print *
print *, 'Lattice is: ', logic%lattice
if (base_name /= logic%lattice) then
  print *, 'Bmad file name: ', trim(base_name)
  print *
  print *, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NOTE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
  print *, 'THE LATTICE NAME DOES NOT MATCH THE NAME OF THE INPUT LATTICE FILE!'
  print *, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
  print *
endif

print *, 'Species is: ', particle_name(u%ring%param%particle)

! setup data base groups...
! with reestablish_groups$ see if groups already exist. If so use them.

call bmad_to_db (u%ring, u%db, 'NOW')

call setup_groups (u, make_groups)

! Save design values

call init_lat (u%design, u%ring%n_ele_max)
u%design = u%ring

!--------------------------------------------------------------------
! variable inits

if (associated(u%dm_dv)) deallocate(u%dm_dv)

u%hsep_kick%v(:)%good_user = .false.
u%bpm_tilt%v(:)%good_user = .false.
u%x_kick_quad%v(:)%good_user = .false.
u%y_kick_quad%v(:)%good_user = .false.
u%hsteer_kick%v(111:)%good_user = .false.
u%vsteer_kick%v(101:)%good_user = .false.

u%hsep_kick%v(:)%do_limit_calc = .false.

u%hsteer_kick%v(:)%makeup_method = normal_makeup$
u%vsteer_kick%v(:)%makeup_method = normal_makeup$
u%hsep_kick%v(:)%makeup_method = normal_makeup$
u%energy_var%v(1)%makeup_method = no_makeup$
u%bpm_tilt%v(:)%makeup_method = no_makeup$

! Set variable limits, etc.
! Limit horizontal hard bend minimums to 1000 cu to make sure PS regulation
! is OK.

u%db%csr_hbnd_cur%cu_low_lim = 1000

call init_var_info (u%hsteer_kick, u%db%csr_horz_cur, u%ring)
call init_var_info (u%hsteer_kick, u%db%csr_hbnd_cur, u%ring)
call init_var_info (u%hsteer_kick, u%db%und_cntg_cur, u%ring)
call init_var_info (u%hsteer_kick, u%db%und_cntgtrim, u%ring)

call init_var_info (u%vsteer_kick, u%db%csr_vert_cur, u%ring)
call init_var_info (u%vsteer_kick, u%db%und_vert_cur, u%ring)

call init_var_info (u%hsep_kick, u%db%csr_hsp_volt, u%ring)

call init_var_info (u%quad_k1,       u%db%quad, u%ring)
call init_var_info (u%x_kick_quad, u%db%quad, u%ring, hkick$, .false.)
call init_var_info (u%y_kick_quad, u%db%quad, u%ring, vkick$, .false.)
u%x_kick_quad%v%name = 'X_KICK_' // u%x_kick_quad%v%name
u%y_kick_quad%v%name = 'Y_KICK_' // u%y_kick_quad%v%name

call init_var_info (u%skew_quad_k1, u%db%csr_sqewquad, u%ring)

call init_var_info (u%sex_k2, u%db%csr_sext_cur, u%ring)

call init_var_info (u%skew_sex_k2, u%db%csr_sqewsext, u%ring)

call init_var_info (u%oct_k3, u%db%csr_octu_cur, u%ring)

do i = lbound(u%hsep_kick%v, 1), ubound(u%hsep_kick%v, 1)
  u%hsep_kick%v(i)%model = 0        ! turn off
  call var_bookkeeper (u%hsep_kick%v(i), u%ring, u%orb)
enddo

! Zero v_seps if they exist

do i = 1, size(u%db%csr_vsp_volt)
  ix = u%db%csr_vsp_volt(i)%ix_lat
  if (ix < 1) cycle
  u%ring%ele(ix)%value(vkick$) = 0
enddo

! init init_orb for an open ring

u%init_orb%v(:)%name = (/ "X Position", "X Velocity", &
                          "Y Position", "Y Velocity", &
                          "Z Position", "Z'  (dE/E)" /)
u%init_orb%v(:)%good_var = .true.
u%init_orb%v(:)%good_user = .false.

! Energy error variable

u%energy_var%v(1)%name = 'Energy var'
u%energy_var%v(1)%design = 0.0
u%energy_var%v(1)%model => u%orb(0)%vec(6)  ! point to orbit vector
u%energy_var%v(1)%exists = .true.
u%energy_var%v(1)%good_var = .false.
u%energy_var%v(1)%good_user = .false. ! don't use initially

!

u%skew_quad_k1%v(:)%good_user = .false.
u%hsep_kick%v(:)%good_user = .false.

! calc cu_design. 
! Turn off vars with dvar_dcu = 0

do i = 1, size(u%var)
  if (.not. u%var(i)%exists) cycle
  if (u%var(i)%dvar_dcu == 0) then
    if (.not. u%var(i)%good_var) cycle
    if (u%var(i)%name(1:8) == 'X_KICK_Q') cycle
    if (u%var(i)%name(1:8) == 'Y_KICK_Q') cycle
    print *, 'Note: dVar_dCU = 0 so Good_var veto of: ' // u%var(i)%name
    u%var(i)%good_var = .false.
  else
    u%var(i)%cu_design = (u%var(i)%design - u%var(i)%base_cu0) / u%var(i)%dvar_dcu
  endif
enddo

u%var(:)%saved = u%var(:)%design
call saved_and_lim_calc (u%var)
call set_var_useit (u)

!------------------------------------------------------------------
! data inits

u%data(:)%exists = .false.
u%data(:)%good_user = .true.

u%q2x%a_in%d(:)%good_user = .true.
u%q2x%a_out%d(:)%good_user = .true.
u%q2y%a_in%d(:)%good_user = .true.
u%q2y%a_out%d(:)%good_user = .true.
u%qx_plus_qy%a_in%d(:)%good_user = .true.
u%qx_plus_qy%a_out%d(:)%good_user = .true.
u%qx_minus_qy%a_in%d(:)%good_user = .true.
u%qx_minus_qy%a_out%d(:)%good_user = .true.

u%beta%x%d(:)%good_user = .false.
u%beta%y%d(:)%good_user = .false.

u%e_xray%x%d(:)%good_user = .false.
u%e_xray%y%d(:)%good_user = .false.

u%p_xray%x%d(:)%good_user = .false.
u%p_xray%y%d(:)%good_user = .false.

u%cmat_a%m12%d(:)%good_user = .false.
u%cmat_a%m22%d(:)%good_user = .false.
u%cmat_b%m11%d(:)%good_user = .false.
u%cmat_b%m12%d(:)%good_user = .false.

u%mode_eta%x%d(:)%good_user = .true.
u%mode_eta%y%d(:)%good_user = .true.

u%energy_data%d1%d(1)%good_user = .true.
u%energy_data%d1%d(1)%exists = .true.
u%energy_data%d1%d(1)%good_dat = .true.
u%energy_data%d1%d(1)%good_ref = .true.

u%chrom%x%d(1)%good_user = .false.
u%chrom%x%d(1)%exists = .true.

u%chrom%y%d(1)%good_user = .false.
u%chrom%y%d(1)%exists = .true.

u%q2x%a_in%d(:)%weight = 1.0e-4;  u%q2x%a_out%d(:)%weight = 1.0e-4
u%q2y%a_in%d(:)%weight = 1.0e-4;  u%q2y%a_out%d(:)%weight = 1.0e-4
u%qx_plus_qy%a_in%d(:)%weight = 1.0e-4;  u%qx_plus_qy%a_out%d(:)%weight = 1.0e-4
u%qx_minus_qy%a_in%d(:)%weight = 1.0e-4; u%qx_minus_qy%a_out%d(:)%weight = 1.0e-4
u%beta%x%d(:)%weight = 1.0e0;   u%beta%y%d(:)%weight = 1.0e0
u%phase%x%d(:)%weight = 1.0e2;  u%phase%y%d(:)%weight = 1.0e2
u%tune%x%d(:)%weight = 1.0e3;   u%tune%y%d(:)%weight = 1.0e3
u%orbit%x%d(:)%weight = 1.0e6;  u%orbit%y%d(:)%weight = 1.0e6
u%e_xray%x%d(:)%weight = 1.0e6;  u%e_xray%y%d(:)%weight = 1.0e6
u%p_xray%x%d(:)%weight = 1.0e6;  u%p_xray%y%d(:)%weight = 1.0e6
u%eta%x%d(:)%weight = 0;        u%eta%y%d(:)%weight = 2.0e2
u%ac_eta%x%d(:)%weight = 0;     u%ac_eta%y%d(:)%weight = 2.0e2
u%ac_eta_c12%a%d(:)%weight = 2.0e2;     u%ac_eta_c12%b%d(:)%weight= 2.0e2
u%ac_eta_yx%yxcos%d(:)%weight = 2.0e2;   u%ac_eta_yx%yxsin%d(:)%weight  = 2.0e2
u%mode_eta%x%d(:)%weight = 0;   u%mode_eta%y%d(:)%weight = 2.0e2
u%energy_data%d1%d(1)%weight = 1e10
u%cbar%m12%d(:)%weight = 1.0e3
u%cbar%m11%d(:)%weight = 1.0e3
u%cbar%m22%d(:)%weight = 1.0e3

u%cmat_a%m12%d(:)%weight = 1.0e3
u%cmat_a%m22%d(:)%weight = 1.0e3
u%cmat_b%m11%d(:)%weight = 1.0e3
u%cmat_b%m12%d(:)%weight = 1.0e3

u%chrom%x%d(:)%weight = 1.0e1;   u%chrom%y%d(:)%weight = 1.0e1

! Orbit, Phase, Eta and Cbar data

do i = 0, 120

  ix = u%db%detector(i)%ix_lat

  u%phase%x%d(i)%ix_ele = ix
  u%phase%y%d(i)%ix_ele = ix

  u%q2x%a_in%d(i)%ix_ele = ix
  u%q2x%a_out%d(i)%ix_ele = ix

  u%q2y%a_in%d(i)%ix_ele = ix
  u%q2y%a_out%d(i)%ix_ele = ix

  u%qx_plus_qy%a_in%d(i)%ix_ele = ix
  u%qx_plus_qy%a_out%d(i)%ix_ele = ix

  u%qx_minus_qy%a_in%d(i)%ix_ele = ix
  u%qx_minus_qy%a_out%d(i)%ix_ele = ix

  u%beta%x%d(i)%ix_ele = ix
  u%beta%y%d(i)%ix_ele = ix

  u%orbit%x%d(i)%ix_ele = ix
  u%orbit%y%d(i)%ix_ele = ix

  u%eta%x%d(i)%ix_ele = ix
  u%eta%y%d(i)%ix_ele = ix

  u%ac_eta%x%d(i)%ix_ele = ix
  u%ac_eta%y%d(i)%ix_ele = ix

  u%ac_eta_c12%a%d(i)%ix_ele = ix
  u%ac_eta_c12%b%d(i)%ix_ele = ix

  u%ac_eta_yx%yxcos%d(i)%ix_ele = ix
  u%ac_eta_yx%yxsin%d(i)%ix_ele = ix

  u%mode_eta%x%d(i)%ix_ele = ix
  u%mode_eta%y%d(i)%ix_ele = ix

  u%cbar%m11%d(i)%ix_ele = ix
  u%cbar%m12%d(i)%ix_ele = ix
  u%cbar%m21%d(i)%ix_ele = ix
  u%cbar%m22%d(i)%ix_ele = ix

  u%cmat_a%m12%d(i)%ix_ele = ix
  u%cmat_a%m22%d(i)%ix_ele = ix
  u%cmat_b%m11%d(i)%ix_ele = ix
  u%cmat_b%m12%d(i)%ix_ele = ix

  u%energy_data%d1%d(i)%ix_ele = ix

  if (ix == 0) cycle  ! Does not exist

  u%phase%x%d(i)%s     = u%ring%ele(ix)%s
  u%phase%y%d(i)%s     = u%ring%ele(ix)%s
  u%q2x%a_in%d(i)%s     = u%ring%ele(ix)%s
  u%q2x%a_out%d(i)%s   = u%ring%ele(ix)%s
  u%q2y%a_in%d(i)%s     = u%ring%ele(ix)%s
  u%q2y%a_out%d(i)%s   = u%ring%ele(ix)%s
  u%qx_plus_qy%a_in%d(i)%s     = u%ring%ele(ix)%s
  u%qx_plus_qy%a_out%d(i)%s   = u%ring%ele(ix)%s
  u%qx_minus_qy%a_in%d(i)%s    = u%ring%ele(ix)%s
  u%qx_minus_qy%a_out%d(i)%s  = u%ring%ele(ix)%s
  u%beta%x%d(i)%s      = u%ring%ele(ix)%s
  u%beta%y%d(i)%s      = u%ring%ele(ix)%s
  u%eta%x%d(i)%s       = u%ring%ele(ix)%s
  u%eta%y%d(i)%s       = u%ring%ele(ix)%s
  u%ac_eta%x%d(i)%s    = u%ring%ele(ix)%s
  u%ac_eta%y%d(i)%s    = u%ring%ele(ix)%s

  u%ac_eta_c12%a%d(i)%s    = u%ring%ele(ix)%s
  u%ac_eta_c12%b%d(i)%s    = u%ring%ele(ix)%s

  u%ac_eta_yx%yxcos%d(i)%s    = u%ring%ele(ix)%s
  u%ac_eta_yx%yxsin%d(i)%s    = u%ring%ele(ix)%s

  u%mode_eta%x%d(i)%s  = u%ring%ele(ix)%s
  u%mode_eta%y%d(i)%s  = u%ring%ele(ix)%s
  u%orbit%x%d(i)%s     = u%ring%ele(ix)%s
  u%orbit%y%d(i)%s     = u%ring%ele(ix)%s
  u%cbar%m11%d(i)%s = u%ring%ele(ix)%s
  u%cbar%m12%d(i)%s = u%ring%ele(ix)%s
  u%cbar%m21%d(i)%s = u%ring%ele(ix)%s
  u%cbar%m22%d(i)%s = u%ring%ele(ix)%s
  u%cmat_a%m12%d(i)%s = u%ring%ele(ix)%s
  u%cmat_a%m22%d(i)%s = u%ring%ele(ix)%s
  u%cmat_b%m11%d(i)%s = u%ring%ele(ix)%s
  u%cmat_b%m12%d(i)%s = u%ring%ele(ix)%s
  u%energy_data%d1%d(i)%s = u%ring%ele(ix)%s

  if (i > 99) cycle   ! Veto all dets above 99.

  u%phase%x%d(i)%exists     = .true.
  u%phase%y%d(i)%exists     = .true.
  u%q2x%a_in%d(i)%exists     = .true.
  u%q2x%a_out%d(i)%exists   = .true.
  u%q2y%a_in%d(i)%exists     = .true.
  u%q2y%a_out%d(i)%exists   = .true.
  u%qx_plus_qy%a_in%d(i)%exists     = .true.
  u%qx_plus_qy%a_out%d(i)%exists   = .true.
  u%qx_minus_qy%a_in%d(i)%exists    = .true.
  u%qx_minus_qy%a_out%d(i)%exists  = .true.
  u%beta%x%d(i)%exists      = .true.
  u%beta%y%d(i)%exists      = .true.
  u%eta%x%d(i)%exists       = .true.
  u%eta%y%d(i)%exists       = .true.
  u%ac_eta%x%d(i)%exists    = .true.
  u%ac_eta%y%d(i)%exists    = .true.
  u%ac_eta_c12%a%d(i)%exists    = .true.
  u%ac_eta_c12%b%d(i)%exists    = .true.
  u%ac_eta_yx%yxcos%d(i)%exists    = .true.
  u%ac_eta_yx%yxsin%d(i)%exists    = .true.
  u%mode_eta%x%d(i)%exists  = .true.
  u%mode_eta%y%d(i)%exists  = .true.
  u%orbit%x%d(i)%exists     = .true.
  u%orbit%y%d(i)%exists     = .true.
  u%cbar%m11%d(i)%exists    = .true.
  u%cbar%m12%d(i)%exists    = .true.
  u%cbar%m21%d(i)%exists    = .true.
  u%cbar%m22%d(i)%exists    = .true.
  u%cmat_a%m12%d(i)%exists  = .true.
  u%cmat_a%m22%d(i)%exists  = .true.
  u%cmat_b%m11%d(i)%exists  = .true.
  u%cmat_b%m12%d(i)%exists  = .true.
  u%energy_data%d1%d(i)%exists = .true.

enddo

print *, 'Note: All detectors with index greater than 99 are vetoed.'

! xray

do i = lbound(u%e_xray%y%d, 1), ubound(u%e_xray%y%d, 1)

  ix = u%db%e_chess_monitor_source(i)%ix_lat

  u%e_xray%x%d(i)%ix_ele = ix
  u%e_xray%y%d(i)%ix_ele = ix

  if (ix == 0) cycle  ! Does not exist

  u%e_xray%x%d(i)%s = u%ring%ele(ix)%s
  u%e_xray%y%d(i)%s = u%ring%ele(ix)%s

  u%e_xray%x%d(i)%ix_db = u%db%e_chess_monitor_source(i)%ix_db
  u%e_xray%y%d(i)%ix_db = u%db%e_chess_monitor_source(i)%ix_db

  u%e_xray%x%d(i)%exists = .true.
  u%e_xray%y%d(i)%exists = .true.

enddo

!

do i = lbound(u%p_xray%y%d, 1), ubound(u%p_xray%y%d, 1)

  ix = u%db%p_chess_monitor_source(i)%ix_lat

  u%p_xray%x%d(i)%ix_ele = ix
  u%p_xray%y%d(i)%ix_ele = ix

  if (ix == 0) cycle  ! Does not exist

  u%p_xray%x%d(i)%s = u%ring%ele(ix)%s
  u%p_xray%y%d(i)%s = u%ring%ele(ix)%s

  u%p_xray%x%d(i)%ix_db = u%db%p_chess_monitor_source(i)%ix_db
  u%p_xray%y%d(i)%ix_db = u%db%p_chess_monitor_source(i)%ix_db

  u%p_xray%x%d(i)%exists = .false.
  u%p_xray%y%d(i)%exists = .true.

enddo

! spline beta

call local_data_transfer (u%beta%x, u%spline_beta%x)
call local_data_transfer (u%beta%y, u%spline_beta%y)

! tunes

u%tune%x%d(1)%exists = .true.; u%tune%x%d(1)%good_dat = .true.
u%tune%y%d(1)%exists = .true.; u%tune%y%d(1)%good_dat = .true.
u%tune%z%d(1)%exists = .true.; u%tune%z%d(1)%good_dat = .true.
u%tune%x%d(1)%ix_ele = u%ring%n_ele_track
u%tune%y%d(1)%ix_ele = u%ring%n_ele_track
u%tune%z%d(1)%ix_ele = u%ring%n_ele_track
u%tune%x%d(1)%s = u%ring%param%total_length
u%tune%y%d(1)%s = u%ring%param%total_length
u%tune%z%d(1)%s = u%ring%param%total_length
u%tune%z%d(1)%good_user = .false.

! plot index

u%orbit%x%d%ix_plot_index = u%orbit%x%d%ix_index
do i = 101, ubound(u%orbit%x%d, 1)
  if (.not. u%orbit%x%d(i)%exists) cycle
  s_last = 0
  ix_last = 0
  do j = 0, 99
    if (.not. u%orbit%x%d(j)%exists) cycle
    if (u%orbit%x%d(j)%s > u%orbit%x%d(i)%s) then
      u%orbit%x%d(i)%ix_plot_index = ix_last + (u%orbit%x%d(j)%ix_index - ix_last) * &
            (u%orbit%x%d(i)%s - s_last) / (u%orbit%x%d(j)%s - s_last) 
      exit
    endif
    s_last  = u%orbit%x%d(j)%s
    ix_last = u%orbit%x%d(j)%ix_index
  enddo
enddo

u%orbit%y%d%ix_plot_index       = u%orbit%x%d%ix_plot_index
u%phase%x%d%ix_plot_index       = u%orbit%x%d%ix_plot_index
u%phase%y%d%ix_plot_index       = u%orbit%x%d%ix_plot_index
u%q2x%a_in%d%ix_plot_index       = u%orbit%x%d%ix_plot_index
u%q2x%a_out%d%ix_plot_index     = u%orbit%x%d%ix_plot_index
u%q2y%a_in%d%ix_plot_index       = u%orbit%x%d%ix_plot_index
u%q2y%a_out%d%ix_plot_index     = u%orbit%x%d%ix_plot_index
u%qx_plus_qy%a_in%d%ix_plot_index       = u%orbit%x%d%ix_plot_index
u%qx_plus_qy%a_out%d%ix_plot_index     = u%orbit%x%d%ix_plot_index
u%qx_minus_qy%a_in%d%ix_plot_index      = u%orbit%x%d%ix_plot_index
u%qx_minus_qy%a_out%d%ix_plot_index    = u%orbit%x%d%ix_plot_index
u%beta%x%d%ix_plot_index        = u%orbit%x%d%ix_plot_index
u%beta%y%d%ix_plot_index        = u%orbit%x%d%ix_plot_index
u%spline_beta%x%d%ix_plot_index = u%orbit%x%d%ix_plot_index
u%spline_beta%y%d%ix_plot_index = u%orbit%x%d%ix_plot_index
u%eta%x%d%ix_plot_index         = u%orbit%x%d%ix_plot_index
u%eta%y%d%ix_plot_index         = u%orbit%x%d%ix_plot_index
u%ac_eta%x%d%ix_plot_index      = u%orbit%x%d%ix_plot_index
u%ac_eta%y%d%ix_plot_index      = u%orbit%x%d%ix_plot_index
u%ac_eta_c12%a%d%ix_plot_index      = u%orbit%x%d%ix_plot_index
u%ac_eta_c12%b%d%ix_plot_index      = u%orbit%x%d%ix_plot_index
u%ac_eta_yx%yxcos%d%ix_plot_index      = u%orbit%x%d%ix_plot_index
u%ac_eta_yx%yxsin%d%ix_plot_index      = u%orbit%x%d%ix_plot_index
u%mode_eta%x%d%ix_plot_index    = u%orbit%x%d%ix_plot_index
u%mode_eta%y%d%ix_plot_index    = u%orbit%x%d%ix_plot_index
u%cbar%m11%d%ix_plot_index     = u%orbit%x%d%ix_plot_index
u%cbar%m12%d%ix_plot_index     = u%orbit%x%d%ix_plot_index
u%cbar%m21%d%ix_plot_index     = u%orbit%x%d%ix_plot_index
u%cbar%m22%d%ix_plot_index     = u%orbit%x%d%ix_plot_index
u%cmat_a%m12%d%ix_plot_index   = u%orbit%x%d%ix_plot_index
u%cmat_a%m22%d%ix_plot_index   = u%orbit%x%d%ix_plot_index
u%cmat_b%m11%d%ix_plot_index   = u%orbit%x%d%ix_plot_index
u%cmat_b%m12%d%ix_plot_index   = u%orbit%x%d%ix_plot_index

! record ele_name

do i = 1, size(u%data)
  ie = u%data(i)%ix_ele
  if (ie == 0) cycle
  u%data(i)%ele_name = u%ring%ele(ie)%name
enddo

! weights and step sizes

u%x_kick_quad%v(:)%weight = 1.0e6;   u%x_kick_quad%v(:)%step = 1.0e-5
u%y_kick_quad%v(:)%weight = 1.0e6;   u%y_kick_quad%v(:)%step = 1.0e-5
u%skew_quad_k1%v(:)%weight = 1.0e4;  u%skew_quad_k1%v(:)%step = 1.0e-4
u%quad_k1%v(:)%weight = 1.0e4;       u%quad_k1%v(:)%step = 1.0e-4
u%custom_var%v(:)%weight = 0.0;      u%custom_var%v(:)%step = 0.0
u%sex_k2%v(:)%weight = 1.0e-1;       u%sex_k2%v(:)%step = 1.0e-3
u%skew_sex_k2%v(:)%weight = 1.0e0;   u%skew_sex_k2%v(:)%step = 1.0e-2
u%hsep_kick%v(:)%weight = 1e6;       u%hsep_kick%v%step = 1.0e-5
u%init_orb%v(:)%weight = 0;          u%init_orb%v(:)%step = 1.0e-5


where (u%hsteer_kick%v%dvar_dcu /= 0) u%hsteer_kick%v%weight = 1.0e-6 / u%hsteer_kick%v%dvar_dcu**2
u%hsteer_kick%v(101:106)%weight = 0     ! hard bends do not have a weight
u%hsteer_kick%v%step = 1.0e-6

where (u%vsteer_kick%v%dvar_dcu /= 0) u%vsteer_kick%v%weight = &
                                   1.0e-6 / u%vsteer_kick%v%dvar_dcu**2
u%vsteer_kick%v%step = 1.0e-6

u%energy_var%v%step = 1e-4

! bpm rotations

do i = lbound(u%bpm_tilt%v, 1), ubound(u%bpm_tilt%v, 1)
  ix_ele = u%orbit%x%d(i)%ix_ele
  u%bpm_tilt%v(i)%model     => u%ring%ele(ix_ele)%value(tilt$)
  u%bpm_tilt%v(i)%name      = 'bpm_tilt'
  u%bpm_tilt%v(i)%ele_name  = u%orbit%x%d(i)%ele_name
  u%bpm_tilt%v(i)%alias     = u%orbit%x%d(i)%alias
  u%bpm_tilt%v(i)%ix_ele    = ix_ele
  u%bpm_tilt%v(i)%attrib_name = 'TILT'
  u%bpm_tilt%v(i)%ix_attrib = tilt$
  u%bpm_tilt%v(i)%s         = u%orbit%x%d(i)%s
  u%bpm_tilt%v(i)%exists    = u%orbit%x%d(i)%exists
  u%bpm_tilt%v(i)%good_var  = u%orbit%x%d(i)%exists
  u%bpm_tilt%v(i)%step      = 1e-5
  u%bpm_tilt%v(i)%weight    = 1e2
enddo

!!!!!!  call taylor2_init (u)

!-------------------------------------------------------
! Calculate data(:)%design and data(:)%design_nonlin.
! first set energy deviation to zero.

energy_var_temp = u%energy_var%v(1)%good_user
dE_temp = u%energy_var%v(1)%model

u%energy_var%v(1)%good_user = .true.
u%energy_data%d1%d(1)%good_user = .false.
u%energy_var%v(1)%model = 0

! turn on pretzel for %design_nonlin calc.

do i = 1, size(u%hsep_kick%v)
  u%hsep_kick%v(i)%model = u%hsep_kick%v(i)%design
enddo

! get data(:)%design_nonlin.

u%data(:)%design_nonlin = u%data(:)%model

! turn off pretzel for %design calc.

do i = 1, size(u%hsep_kick%v)
  u%hsep_kick%v(i)%model = 0
enddo

if (logic%rf_on) call set_on_off (rfcavity$, u%ring, on$)

call closed_orbit_calc (u%ring, u%orb, 4)
call radiation_integrals (u%ring, u%orb, u%global_design)
logic%l_times_alpha = u%global_design%synch_int(1)    ! needed for ring_calc

call ring_calc (u)

!--------------------------------------------------------
! check for consistancy

do i = 1, size(u%var)
  if (u%var(i)%dvar_dcu /= 0 .and. .not. u%var(i)%exists) then
    print *, 'Warning: dVar_dCU not zero but element does not exist! '
    print '(9x, a, 5x, a, i4)', u%var(i)%name, &
                  u%var(i)%db_node_name, u%var(i)%ix_db
  endif
enddo

! get data(:)%design

call ring_calc (u)   ! to make sure ring calc is self-consistant

u%data(:)%design = u%data(:)%model

call set_on_off (rfcavity$, u%ring, on$)
call radiation_integrals (u%ring, u%orb, u%global_design, rad_int_by_ele = u%rad_int_by_ele_design)
if (.not. logic%rf_on) call set_on_off (rfcavity$, u%ring, off$)

call chrom_calc (u%ring, 1.0e-4_rp, u%global_design%a%chrom, u%global_design%b%chrom)

! restore energy control

u%energy_var%v(1)%good_user = energy_var_temp
u%energy_data%d1%d(1)%good_user = .not. energy_var_temp
u%energy_var%v(1)%model = dE_temp

!---------------------------------------------------------------------------
contains

subroutine local_data_transfer (d1_in, d1_out)

type (d1_data_struct) d1_in, d1_out

d1_out%d(:)%ix_index  = d1_in%d(:)%ix_index
d1_out%d(:)%ix_ele    = d1_in%d(:)%ix_ele
d1_out%d(:)%exists    = d1_in%d(:)%exists
d1_out%d(:)%s         = d1_in%d(:)%s
d1_out%d(:)%good_dat  = .false.
d1_out%d(:)%good_user = d1_in%d(:)%good_user

end subroutine

end subroutine