program bmadz

  use bmad
  use bunchcross_mod
  use zquad_lens_mod
  use pretz_mod
  use twiss_max_mod
  use constraints_mod
  use nonlin_mod
  use bmadz_utils
  use opt_com_mod
  use bmadz_interface
  use opt_mod
  use geodesic_lm_mod

  implicit none
  
  type (indep_var_struct)    indep_var, sext_var
  type (constraint_struct)   con
  type (lat_struct)         lat_0, lat_1
  type (lat_struct):: save, lat_in, lat_out, lat_temp
  type (coord_struct)        coords_0
  type (coord_struct), allocatable :: co(:)
  type (coord_array_struct) co_arr(0:0)
  type (db_struct)          cesr_i
  type (pc_struct)           pc
  type (twiss_max_struct)    twiss_max
  type (pretz_struct)        pretz
  type (global_struct)       global
  type (zquad_struct)        quad
  type (nonlin_ele_struct)    nonlin
  type (moment_struct)         moments

  procedure(track1_custom_def) :: track1_custom
  procedure(track1_preprocess_def) :: track1_preprocess

  character*60 con_file, line
  character*20 plot_type
  character*200 lattice

  real(rp) ftol, fret
  real(rp) current

  integer lun, i
  integer ix, istat, lib$get_foreign, nret
  integer iter
  integer j,k
  integer terms
  integer particle, i_train, j_car

  logical error_type
  logical use_taylor/.false./
  logical rec_taylor



  integer nargs, iargc
  character arg*80 

! parse command line arguments

  track1_custom_ptr => track1_custom
  track1_preprocess_ptr => track1_preprocess

  call get_cl_input(con_file,con)

  if (con%plot) then
    print*, "Plotting enabled."
  else
    print*, "Plotting disabled."
  end if

! get lat structure and make hybrid lat with only the elements needed.
  call bmad_parser( con%lat_file, lat_0 )

  if(con%define_special_variables /= "") &
       call bmad_parser2( con%define_special_variables, lat_0 )


  if (con%current(1) /= 0.) then
    particle = positron$
    i_train = 1
    rec_taylor = .true.
    lat_in = lat_0
    call lrbbi_setup( &
         lat_in, lat_temp, particle, i_train, con%j_car, &
         con%n_trains, con%n_cars, con%current(1), rec_taylor)
    print *,' lat_temp%n_ele_track = ',lat_temp%n_ele_track
    allocate( co(0:lat_0%n_ele_max) )
    
    particle = electron$
    do i=1,lat_temp%n_ele_track
      if (lat_temp%ele(i)%key == beambeam$ .and. &
           lat_temp%ele(i)%type == 'POSITRON') &
           lat_temp%ele(i)%value(charge$) = 0
    end do
    call lrbbi_setup( &
         lat_temp, lat_out, particle, i_train, con%j_car, &
         con%n_trains, con%n_cars, con%current(1), rec_taylor )
    lat_0 = lat_out
    ! default will be to turn on LRBBI for positrons and off for electrons

    lat_0%param%particle = positron$
    do i=1,lat_0%n_ele_track
      if(lat_0%ele(i)%key == beambeam$)then
        if (lat_0%ele(i)%type == 'ELECTRON') lat_0%ele(i)%value(charge$) = 0
        if (lat_0%ele(i)%type == 'POSITRON') lat_0%ele(i)%value(charge$) = -1
      endif
    end do


  endif



  print *, ' Energy= ', lat_0%ele(0)%value(E_TOT$)

  allocate( co_arr(0)%orbit(0:lat_0%n_ele_max) )

  co_arr(0)%orbit(:)%vec(1) = 0
  co_arr(0)%orbit(:)%vec(2) = 0
  co_arr(0)%orbit(:)%vec(3) = 0
  co_arr(0)%orbit(:)%vec(4) = 0
  co_arr(0)%orbit(:)%vec(5) = 0
  co_arr(0)%orbit(:)%vec(6) = 0

!

  call set_on_off( rfcavity$, lat_0, off$, co_arr(0)%orbit )
  call calc_z_tune( lat_0%branch(0) )

  if (con%linearize_custom) then
    do i = 1,lat_0%n_ele_track
      if (lat_0%ele(i)%key == custom$) lat_0%ele(i)%type = 'LINEAR'
    end do
  endif

  if (con%ele_names(1) == ' ') then
     con%ele_names(1) = 'V_SEP*'
     con%ele_names(2) = 'H_SEP*'
     con%ele_names(3) = 'IP_L3*'
     con%ele_names(4) = 'Q*'
     con%ele_names(5) = 'B*'
     con%ele_names(6) = 'WIG*'
     con%ele_names(7) = 'SOL*'
     con%ele_names(8) = 'IP_L0*'
     con%ele_names(9) = 'SK*'
     con%ele_names(10) = 'SE*'
     con%ele_names(11) = 'RF*'
     con%ele_names(12) = 'DQ*'
  endif

  call name_to_list(lat_0, con%ele_names(1:12))

  call read_indep_var(con_file, lat_0, 'USE_LIST', indep_var, .true.)
  do i = 1, indep_var%n_var
    ix = indep_var%v(i)%rindex
    lat_0%ele(ix)%select = .true.
  enddo

  forall (i=1:lat_0%n_ele_track) co_arr(0)%orbit(i)%vec = 0
  call make_hybrid_lat(lat_0, lat_1, use_taylor, co_arr)
  
  if (.not. con%hybridize) lat_1 = lat_0

  call bmad_to_cesr (lat_1, cesr_i)

  !! call separator_symmetry(cesr_i, lat_1)   ! Separators removed for chess-u upgrade.
  if (con%nonlinearity) then
    call set_nonlin( nonlin, con%disp_file, lat_1%n_ele_max )
    nonlin%n_det_calc = con%n_det_calc
    nonlin%energy_offset(1) = con%energy_offset(1)
    nonlin%energy_offset(2) = con%energy_offset(2)
    nonlin%energy_offset(3) = con%energy_offset(3)
  else
    call allocate_nonlin_ele( nonlin, lat_1%n_ele_max )
  endif

! do some bookkeeping

  if (con%pre_init_set_file /= "") &
       call bmad_parser2( con%pre_init_set_file, lat_1 )

  call init_ele_init( con, lat_1, co_arr(0)%orbit )
  if (con%initial_lat_file == "") then
    lat_0 = lat_1
  else
    call bmad_parser( con%initial_lat_file, lat_0 )
    if (con%pre_init_set_file /= "") &
         call bmad_parser2( con%pre_init_set_file, lat_0 )
    if (con%current(1) /= 0.) then
      current = 0.001
      lat_in = lat_0
      call lrbbi_setup( &
           lat_in, lat_out, particle, i_train, con%j_car, &
           con%n_trains, con%n_cars, current, rec_taylor )
      lat_0 = lat_out
    endif

    co_arr(0)%orbit(0)%species = lat_0%param%particle

    call twiss_at_start( lat_0 )
    call closed_orbit_calc( lat_0, co_arr(0)%orbit, 4, +1 )
    call track_all( lat_0, co_arr(0)%orbit )
    call lat_make_mat6( lat_0, -1, co_arr(0)%orbit )
    call twiss_at_start( lat_0 )
    call twiss_propagate_all( lat_0 )
  endif

  bmad_com%aperture_limit_on = .false.  !dlr 1/13/05
  call identify_quads( lat_1, quad )
  
  call bunchcross( lat_1, con, pc )
  if (con%pretzel) call parasitic_crossing_points( lat_1, pc )

  call calcf( lat_1, con, pc, quad, co_arr(0)%orbit, global, twiss_max, pretz, nonlin )
  call calc_z_tune( lat_1%branch(0) )

  call map_constraint_to_lat_element( con, lat_1, co_arr(0)%orbit, lat_0 )
  call read_indep_var( con_file, lat_1, 'USE_LIST', indep_var, .false. )
  call read_indep_var( con_file, lat_1, 'USE_LIST_SEXT', sext_var, .true. )

  call element_list( lat_1 )

  if (con%plot .and. con%n_loops /= 0) &
       call plotdo_bmadz( '/xs', lat_1, co_arr(0)%orbit, pc )

  if (con%post_init_set_file /= "") &
       call bmad_parser2( con%post_init_set_file, lat_1 )

  call allocate_moment( moments, lat_1%n_ele_max )

  call output_direct( -1, .false. )

! optimize lat

if (con%optimizer == 'GEOLEVMAR') then
      call geodesic_lm_opt( lat_1, indep_var, cesr_i, con, pc, quad,  &
           co_arr(0)%orbit, global, twiss_max, pretz, nonlin, moments, sext_var )
else

  do i=1,con%n_loops

    if (con%optimizer == 'MINIDEL') call bmad_opt (lat_1, indep_var,  &
         cesr_i, con,  pc, quad, co_arr(0)%orbit, global, twiss_max, pretz,  &
         nonlin, moments, sext_var)

    if (con%optimizer == 'MINOP') call bmad_minop( lat_1, indep_var,  &
         cesr_i, con,  pc, quad, co_arr(0)%orbit, global, twiss_max, pretz,  &
         nonlin, moments, sext_var)

    if (con%optimizer == 'LINEAR') call linear_opt( lat_1, indep_var,  &
         cesr_i, con,  pc, quad, co_arr(0)%orbit, global, twiss_max, pretz,  &
         nonlin, moments, sext_var )

    if (con%optimizer == 'FRPRMN') then
      ftol=1.e-7
      call frprmn_opt( lat_1, indep_var, cesr_i, con, pc, quad, co_arr(0)%orbit, global, &
           twiss_max, pretz, nonlin, moments, sext_var, ftol, iter, fret )
      print *,' ITER = ',iter,'    FRET =',fret
    endif

    if (con%optimizer == 'MARQUARDT') then
      call marquardt_opt( lat_1, indep_var, cesr_i, con, pc, quad,  &
           co_arr(0)%orbit, global, twiss_max, pretz, nonlin, moments, sext_var )
    endif

       if (con%special_output /= ' ') then
      call special_output( lat_0, lat_1, indep_var, con )
    endif

    print *, 'End of Loop:', i
        
  end do

endif  

  if (con%optimizer == 'MANUAL') call manual_opt( lat_1, indep_var,  &
       cesr_i, con,  pc, quad, co_arr(0)%orbit, global, twiss_max, pretz,  &
       nonlin, moments, sext_var )



! end stuff

  call calcf( lat_1, con,  pc, quad, co_arr(0)%orbit, global, twiss_max, &
       pretz, nonlin, indep_var )

  call fom( lat_1, cesr_i, con, pc, quad, co_arr(0)%orbit, global, twiss_max,  &
       pretz ,nonlin, moments, indep_var, sext_var )

!  if(con%nonlinearity)call write_moments(lat_1, co_arr(0)%orbit, nonlin, moments)

  if (con%minimize_moments) then

    call calcf( lat_1, con, pc, quad, co_arr(0)%orbit, global, twiss_max, &
         pretz, nonlin, indep_var )

    con%minimize_moments = .false.
    
    call fom( lat_1, cesr_i, con, pc, quad, co_arr(0)%orbit, global, twiss_max,  &
         pretz ,nonlin, moments, indep_var, sext_var )

    con%minimize_moments = .true.

  endif


  call showme( con, lat_1, 0 )

  call string_trim( con_file, con_file, ix )
  call str_upcase( con_file, con_file )
  if (index(con_file, 'VALIDATOR') /= 0) call validator( con, lat_1, 0 )

  call zrep_out( lat_1, co_arr(0)%orbit, pc, global, twiss_max, &
       pretz, cesr_i, nonlin, moments, con )

  if(con%plot) then
    call plotdo_bmadz('/xs', lat_1, co_arr(0)%orbit, pc, nonlin)
    print '(a, $)',' postscript (P) or gif (G) ? (def = none)'
    read(*,'(a)') plot_type
    call string_trim( plot_type, plot_type, ix )
    call str_upcase( plot_type, plot_type )
    if (plot_type == 'P' .or. plot_type == 'G') &
         call plotdo_bmadz( plot_type, lat_1, co_arr(0)%orbit, pc, nonlin )
  endif

  deallocate( co_arr(0)%orbit )
  call deallocate_nonlin_ele( nonlin )
  call deallocate_moment( moments )
  
end program bmadz