Better set_tune_via_group_knobs algorithm.

bmad/code/set_tune_via_group_knobs.f90

@@ -31,7 +31,7 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err)
 
 real(rp) phi_set(2), dphi_a, dphi_b, dQ_max
 real(rp) phi_a, phi_b, d_a1, d_a2, d_b1, d_b2, det
-real(rp) d1, d2, del0
+real(rp) d1, d2, del0, factor
 real(rp) :: phi_array(2)
 real(rp), allocatable :: deriv1(:), deriv2(:), kinit(:)
 
@@ -48,6 +48,7 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err)
 dQ_max = 0.001
 del0 = 0.001
 ok = .false.
+factor = 1
 rf_on = rf_is_on(branch)
 
 allocate(kinit(branch%n_ele_track), deriv1(branch%n_ele_track), deriv2(branch%n_ele_track))
@@ -58,7 +59,7 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err)
 
 ! Q tune
 
-do i = 1, 10
+do i = 1, 20
   call lattice_bookkeeper(branch%lat)
 
   if (rf_on) then
@@ -72,7 +73,19 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err)
   call lat_make_mat6 (branch%lat, -1, orb, branch%ix_branch)
 
   call twiss_at_start(branch%lat, status, branch%ix_branch, print_err)
-  if (status /= ok$) return
+  if (status /= ok$) then
+    if (i == 1) return
+    group1%control%var(1)%value = group1%control%var(1)%value - factor * d1
+    group2%control%var(1)%value = group2%control%var(1)%value - factor * d2
+    factor = 0.5 * factor
+    group1%control%var(1)%value = group1%control%var(1)%value + factor * d1
+    group2%control%var(1)%value = group2%control%var(1)%value + factor * d2
+    call set_flags_for_changed_attribute (group1)
+    call set_flags_for_changed_attribute (group2)
+    cycle
+  else
+    factor = 1.0
+  endif
 
   call twiss_propagate_all (branch%lat, branch%ix_branch, err)
   if (err) return
@@ -113,8 +126,8 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err)
 
   ! Put in the changes
 
-  group1%control%var(1)%value = group1%control%var(1)%value + d1
-  group2%control%var(1)%value = group2%control%var(1)%value + d2
+  group1%control%var(1)%value = group1%control%var(1)%value + factor * d1
+  group2%control%var(1)%value = group2%control%var(1)%value + factor * d2
 
   call set_flags_for_changed_attribute (group1)
   call set_flags_for_changed_attribute (group2)

bmad/multiparticle/beam_file_io.f90

@@ -10,22 +10,26 @@ module beam_file_io
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
 !+
-! Subroutine write_beam_file (file_name, beam, new_file, file_format, lat)
+! Subroutine write_beam_file (file_name, beam, new_file, file_format, lat, cols)
 !
 ! Routine to write a beam file.
 !
 ! A '.h5' suffix will be appended to the created file if hdf5$ format is used and file_name does not
 ! already have a '.h5' or '.hdf5' suffix.
 !
+! The ASCII format is ASCII::4 except if cols is present. In this case ASCII::5 is used.
+!
 ! Input:
 !   file_name     -- character(*): Name of file.
 !   beam          -- beam_struct: Beam to write
 !   new_file      -- logical, optional: New file or append? Default = True.
 !   file_format   -- logical, optional: ascii$, or hdf5$ (default).
 !   lat           -- lat_struct, optional: If present, lattice info will be writen to hdf5 files.
+!   cols          -- character(*): Table columns. Possibilities are coord_struct component names.
+!                      column names should be separated by spaces.
 !-
 
-subroutine write_beam_file (file_name, beam, new_file, file_format, lat)
+subroutine write_beam_file (file_name, beam, new_file, file_format, lat, cols)
 
 type (beam_struct), target :: beam
 type (bunch_struct), pointer :: bunch
@@ -36,6 +40,7 @@ subroutine write_beam_file (file_name, beam, new_file, file_format, lat)
 integer, optional :: file_format
 
 character(*) file_name
+character(*), optional :: cols
 character(200) full_name
 character(*), parameter :: r_name = 'write_beam_file'
 
@@ -58,6 +63,11 @@ subroutine write_beam_file (file_name, beam, new_file, file_format, lat)
 
 !
 
+if (present(cols)) then
+  call write_ascii4_beam_file(file_name, beam, new_file, cols, lat)
+  return
+endif
+
 if (logic_option(.true., new_file)) then
   open (iu, file = full_name)
   write (iu, '(a)') '!ASCII::3'
@@ -87,6 +97,64 @@ subroutine write_beam_file (file_name, beam, new_file, file_format, lat)
 
 end subroutine write_beam_file
 
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+!+
+! Subroutine write_ascii4_beam_file (file_name, beam, new_file, cols, lat)
+!
+! Routine to write a beam file in ASCII::4 format.
+!
+! A '.h5' suffix will be appended to the created file if hdf5$ format is used and file_name does not
+! already have a '.h5' or '.hdf5' suffix.
+!
+! The ASCII format is ASCII::4 except if cols is present. In this case ASCII::5 is used.
+!
+! Input:
+!   file_name     -- character(*): Name of file.
+!   beam          -- beam_struct: Beam to write
+!   new_file      -- logical, optional: New file or append? Default = True.
+!   cols          -- character(*): Table columns. Possibilities are coord_struct component names.
+!                      column names should be separated by spaces.
+!   lat           -- lat_struct, optional: If present, lattice info will be writen to hdf5 files.
+!-
+
+subroutine write_ascii4_beam_file (file_name, beam, new_file, cols, lat)
+
+type (beam_struct), target :: beam
+type (bunch_struct), pointer :: bunch
+type (coord_struct), pointer :: p
+type (lat_struct), optional :: lat
+
+integer j, iu, ib, ip, ix_ele, n, n0
+
+character(*) file_name
+character(*), optional :: cols
+character(200) full_name
+character(*), parameter :: r_name = 'write_ascii4_beam_file'
+
+logical, optional :: new_file
+logical error, append
+
+!
+
+iu = lunget()
+call fullfilename (file_name, full_name)
+
+if (logic_option(.true., new_file)) then
+  open (iu, file = full_name)
+  write (iu, '(a)') '!ASCII::3'
+else
+  open (iu, file = full_name, access = 'append')
+endif
+
+!
+
+
+close (iu)
+
+end subroutine write_ascii4_beam_file
+
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
@@ -636,10 +704,15 @@ subroutine read_beam_ascii4 (iu, file_name, beam, beam_init, err_flag, ele, prin
       call read_switch(line(:ix), p0%species, str, err)
       if (err) return
 
-    case ('s_position');  p0%s           = read_param(line)
-    case ('time');        p0%t           = read_param(line)
+    case ('r');           p0%r           = read_param(line)
+    case ('s');           p0%s           = read_param(line)
+    case ('t');           p0%t           = read_param(line)
     case ('p0c');         p0%p0c         = read_param(line)
+    case ('dt_ref');      p0%dt_ref      = read_param(line)
     case ('charge');      p0%charge      = read_param(line)
+    case ('spin');        call read_params(line, p0%spin)
+    case ('field');       call read_params(line, p0%field)
+    case ('phase');       call read_params(line, p0%phase)
     case ('time_dir');    p0%time_dir    = nint(read_param(line))
     case ('direction');   p0%direction   = nint(read_param(line))
     case ('ix_ele');      p0%ix_ele      = nint(read_param(line))
@@ -839,6 +912,24 @@ end function read_param
 !---------------------------------------------------------------------------------------------------
 ! contains
 
+subroutine read_params(line, param)
+character(*) line
+real(rp) param(:)
+integer ix, ios
+
+!
+
+ix = index(line, '=')
+read(line(ix+1:), *, iostat = ios) param
+if (ios /= 0 .or. ix == 0) then
+  call out_io (s_error$, r_name, 'ERROR READING BEAM FILE PARAMETER!')
+endif
+
+end subroutine read_params
+
+!---------------------------------------------------------------------------------------------------
+! contains
+
 function read_string(line) result (str)
 character(*) line
 character(200) str

bsim/modules/ts_mod.f90

@@ -14,11 +14,11 @@ module ts_mod
   character(100) :: lat_file = '', dat_out_file = '', quad_mask = ''
   character(40) :: group_knobs(2) = ["", ""]
   character(40) :: test = ""
-  real(rp) :: Q_a0 = 0, Q_a1 = 0, dQ_a = 0
-  real(rp) :: Q_b0 = 0, Q_b1 = 0, dQ_b = 0
-  real(rp) :: Q_z0 = 0, Q_z1 = 0, dQ_z = 0
+  real(rp) :: Q_a0 = real_garbage$, Q_a1 = real_garbage$, dQ_a = real_garbage$
+  real(rp) :: Q_b0 = real_garbage$, Q_b1 = real_garbage$, dQ_b = real_garbage$
+  real(rp) :: Q_z0 = real_garbage$, Q_z1 = real_garbage$, dQ_z = real_garbage$
   real(rp) :: pz0 = 0, pz1 = 0, dpz = 0
-  real(rp) :: a_emit = 0, b_emit = 0, sig_pz
+  real(rp) :: a_emit = 0, b_emit = 0, sig_pz = 0
   real(rp) :: a0_amp = 0, b0_amp = 0, pz0_amp = 0
   real(rp) :: timer_print_dtime = 120
   integer :: n_turn = 0
@@ -80,13 +80,31 @@ subroutine ts_init_params (ts, ts_com)
 ts_com%master_input_file = 'tune_scan.init'
 if (n_arg == 1) call get_command_argument(1, ts_com%master_input_file)
 
+bmad_com%auto_bookkeeper = .false.
+
 open (unit= 1, file = ts_com%master_input_file, status = 'old', action = 'read')
 read(1, nml = params)
 close (1)
 
-ts%Q_z0 = abs(ts%Q_z0)
-ts%Q_z1 = abs(ts%Q_z1)
-ts%dQ_z = abs(ts%dQ_z)
+if (ts%Q_a0 == real_garbage$ .or. ts%Q_a1 == real_garbage$ .or. ts%dQ_a == real_garbage$) then
+  print '(a)', 'Error: One of ts%Q_a0, ts%Q_a1, ts%dQ_a is not set! Stopping here.'
+  stop
+endif
+
+if (ts%Q_b0 == real_garbage$ .or. ts%Q_b1 == real_garbage$ .or. ts%dQ_b == real_garbage$) then
+  print '(a)', 'Error: One of ts%Q_b0, ts%Q_b1, ts%dQ_b is not set! Stopping here.'
+  stop
+endif
+
+if (ts%rf_on) then
+  if (ts%Q_z0 == real_garbage$ .or. ts%Q_z1 == real_garbage$ .or. ts%dQ_z == real_garbage$) then
+    print '(a)', 'Error: One of ts%Q_z0, ts%Q_z1, ts%dQ_z is not set! Stopping here.'
+    stop
+  endif
+  ts%Q_z0 = abs(ts%Q_z0)
+  ts%Q_z1 = abs(ts%Q_z1)
+  ts%dQ_z = abs(ts%dQ_z)
+endif
 
 if (ts%dat_out_file == '') call file_suffixer(ts_com%master_input_file, ts%dat_out_file, 'dat', .true.)
 
@@ -201,7 +219,7 @@ recursive subroutine ts_track_particle(ts, ts_com, jtune, ts_dat)
 type (ele_struct), pointer :: ele
 
 real(rp) Jvec0(1:6), Jvec(1:6), init_vec(6), amp(3), r(3), v_mat(4,4)
-integer jtune(3), nt, track_state, status
+integer jtune(3), nt, track_state, status, iqm
 logical ok, error
 
 !
@@ -213,15 +231,19 @@ recursive subroutine ts_track_particle(ts, ts_com, jtune, ts_dat)
 ts_dat%tune(2) = ts_com%int_Qb + ts%Q_b0 + jtune(2)*ts%dQ_b
 if (ts%rf_on) then
   ts_dat%tune(3) = -(ts%Q_z0 + jtune(3)*ts%dQ_z)
+  iqm = 3
 else
   ts_dat%tune(3) = ts%pz0 + jtune(3)*ts%dpz
+  iqm = 2
 endif
 
 ring = ts_com%ring              ! Use copy in case tune setting fails, which may garble the lattice
 closed_orb = ts_com%closed_orb  ! Use copy in case tune setting fails, which may garble the closed orbit
 ele => ring%ele(0)
 
 ok = set_tune_3d (ring%branch(0), ts_dat%tune, ts%quad_mask, ts%use_phase_trombone, ts%rf_on, ts%group_knobs)  ! Tunes in radians.
+if (.not. ok) print '(a, 3f9.4)', 'Cannot set tunes (due to resonance?) for: ', ts_dat%tune(1:iqm)
+
 if (.not. ok) return    ! Tunes could not be set, probably on a resonance.
 
 if (ts%rf_on) then

tao/version/tao_version_mod.f90

@@ -6,5 +6,5 @@
 !-
 
 module tao_version_mod
-character(*), parameter :: tao_version_date = "2023/10/16 23:05:40"
+character(*), parameter :: tao_version_date = "2023/10/24 11:25:06"
 end module