Now Tao chrom data evaluated using local dpz. (#1220)

bmad/code/chrom_calc.f90

@@ -1,5 +1,5 @@
 !+
-! Subroutine chrom_calc (lat, delta_e, chrom_x, chrom_y, err_flag, &
+! Subroutine chrom_calc (lat, delta_e, chrom_x, chrom_y, err_flag,
 !                        pz, low_E_lat, high_E_lat, low_E_orb, high_E_orb, ix_branch, orb0)
 !
 ! Subroutine to calculate the chromaticities by computing the tune change when the energy is changed.

tao/code/tao_evaluate_a_datum.f90

@@ -61,7 +61,7 @@ recursive subroutine tao_evaluate_a_datum (datum, u, tao_lat, datum_value, valid
 real(rp) datum_value, mat6(6,6), vec0(6), angle, px, py, vec2(2)
 real(rp) eta_vec(4), v_mat(4,4), v_inv_mat(4,4), a_vec(4), mc2, charge
 real(rp) beta_gamma, one_pz, xi_sum, xi_diff, w0_mat(3,3), w_mat(3,3), vec3(3), value, s_len, n0(3)
-real(rp) dz, dx, cos_theta, sin_theta, zz_pt, xx_pt, zz0_pt, xx0_pt, dE, s_offset
+real(rp) dz, dx, cos_theta, sin_theta, zz_pt, xx_pt, zz0_pt, xx0_pt, dpz, s_offset
 real(rp) zz_center, xx_center, xx_wall, phase, amp, dalpha, dbeta, aa, bb, g2
 real(rp) xx_a, xx_b, dxx1, dzz1, drad, ang_a, ang_b, ang_c, dphi, amp_a, amp_b
 real(rp), allocatable :: value_vec(:)
@@ -719,8 +719,6 @@ recursive subroutine tao_evaluate_a_datum (datum, u, tao_lat, datum_value, valid
 
   !----
 
-  dE = 2 * s%global%delta_e_chrom  ! Actually this is the change in pz
-
   select case (data_type)
 
   case ('chrom.dtune.a', 'chrom.a')
@@ -736,63 +734,71 @@ recursive subroutine tao_evaluate_a_datum (datum, u, tao_lat, datum_value, valid
   case ('chrom.dbeta.a')
     if (data_source == 'lat') then
       do i = ix_start, ix_ele
-        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%a%beta - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%a%beta) / (tao_lat%lat%ele(i)%a%beta * dE)
+        dpz = tao_branch%high_E_orb(i)%vec(6) - tao_branch%low_E_orb(i)%vec(6)
+        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%a%beta - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%a%beta) / (tao_lat%lat%ele(i)%a%beta * dpz)
       end do
       call tao_load_this_datum (value_vec, ele_ref, ele_start, ele, datum_value, valid_value, datum, branch, why_invalid)
     endif
 
   case ('chrom.dbeta.b')
     if (data_source == 'lat') then
       do i = ix_start, ix_ele
-        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%b%beta - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%b%beta) / (tao_lat%lat%ele(i)%b%beta * dE)
+        dpz = tao_branch%high_E_orb(i)%vec(6) - tao_branch%low_E_orb(i)%vec(6)
+        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%b%beta - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%b%beta) / (tao_lat%lat%ele(i)%b%beta * dpz)
       end do
       call tao_load_this_datum (value_vec, ele_ref, ele_start, ele, datum_value, valid_value, datum, branch, why_invalid)
     endif
 
   case ('chrom.dphi.a')
     if (data_source == 'lat') then
       do i = ix_start, ix_ele
-        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%a%phi - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%a%phi)/ dE
+        dpz = tao_branch%high_E_orb(i)%vec(6) - tao_branch%low_E_orb(i)%vec(6)
+        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%a%phi - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%a%phi)/ dpz
       end do
       call tao_load_this_datum (value_vec, ele_ref, ele_start, ele, datum_value, valid_value, datum, branch, why_invalid)
     endif
 
   case ('chrom.dphi.b')
     if (data_source == 'lat') then
       do i = ix_start, ix_ele
-        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%b%phi - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%b%phi)/ dE
+        dpz = tao_branch%high_E_orb(i)%vec(6) - tao_branch%low_E_orb(i)%vec(6)
+        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%b%phi - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%b%phi)/ dpz
       end do
       call tao_load_this_datum (value_vec, ele_ref, ele_start, ele, datum_value, valid_value, datum, branch, why_invalid)
     endif
 
   case ('chrom.deta.x')
     if (data_source == 'lat') then
       do i = ix_start, ix_ele
-        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%x%eta - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%x%eta)/ dE
+        dpz = tao_branch%high_E_orb(i)%vec(6) - tao_branch%low_E_orb(i)%vec(6)
+        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%x%eta - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%x%eta)/ dpz
       end do
       call tao_load_this_datum (value_vec, ele_ref, ele_start, ele, datum_value, valid_value, datum, branch, why_invalid)
     endif
 
   case ('chrom.deta.y')
     if (data_source == 'lat') then
       do i = ix_start, ix_ele
-        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%y%eta - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%y%eta)/ dE
+        dpz = tao_branch%high_E_orb(i)%vec(6) - tao_branch%low_E_orb(i)%vec(6)
+        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%y%eta - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%y%eta)/ dpz
       end do
       call tao_load_this_datum (value_vec, ele_ref, ele_start, ele, datum_value, valid_value, datum, branch, why_invalid)
     endif
 
   case ('chrom.detap.x')
     if (data_source == 'lat') then
       do i = ix_start, ix_ele
-        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%x%etap - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%x%etap)/ dE
+        dpz = tao_branch%high_E_orb(i)%vec(6) - tao_branch%low_E_orb(i)%vec(6)
+        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%x%etap - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%x%etap)/ dpz
       end do
       call tao_load_this_datum (value_vec, ele_ref, ele_start, ele, datum_value, valid_value, datum, branch, why_invalid)
     endif
 
   case ('chrom.detap.y')
     if (data_source == 'lat') then
       do i = ix_start, ix_ele
-        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%y%etap - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%y%etap)/ dE
+        dpz = tao_branch%high_E_orb(i)%vec(6) - tao_branch%low_E_orb(i)%vec(6)
+        value_vec(i) = (tao_lat%high_E_lat%branch(ix_branch)%ele(i)%y%etap - tao_lat%low_E_lat%branch(ix_branch)%ele(i)%y%etap)/ dpz
       end do
       call tao_load_this_datum (value_vec, ele_ref, ele_start, ele, datum_value, valid_value, datum, branch, why_invalid)
     endif
@@ -809,8 +815,9 @@ recursive subroutine tao_evaluate_a_datum (datum, u, tao_lat, datum_value, valid
           z1 => tao_lat%low_E_lat%branch(ix_branch)%ele(i)%b
           z0 => branch%ele(i)%b
         endif
-        dalpha = (z2%alpha - z1%alpha) / dE
-        dbeta  = (z2%beta - z1%beta) / dE
+        dpz = tao_branch%high_E_orb(i)%vec(6) - tao_branch%low_E_orb(i)%vec(6)
+        dalpha = (z2%alpha - z1%alpha) / dpz
+        dbeta  = (z2%beta - z1%beta) / dpz
         aa = dalpha - z0%alpha * dbeta / z0%beta
         bb = dbeta / z0%beta
         value_vec(i) = sqrt(aa**2 + bb**2)

tao/code/tao_init.f90

@@ -322,7 +322,7 @@ subroutine tao_init (err_flag)
 
       if (branch%param%geometry == closed$ .and. tao_branch%track_state == moving_forward$) then
         call chrom_calc (tao_lat%lat, s%global%delta_e_chrom, tao_branch%a%chrom, tao_branch%b%chrom, err, &
-                   tao_branch%orbit(0)%vec(6), low_E_lat=tao_lat%low_E_lat, high_E_lat=tao_lat%high_E_lat, ix_branch = ib)
+              tao_branch%orbit(0)%vec(6), tao_lat%low_E_lat, tao_lat%high_E_lat, tao_branch%low_E_orb, tao_branch%high_E_orb, ib)
         call emit_6d(branch%ele(0), .false., tao_branch%modes_6d, sigma, tao_branch%orbit)
         call emit_6d(branch%ele(0), .true., tao_branch%modes_6d, sigma, tao_branch%orbit)
         tao_branch%modes_6d%momentum_compaction = momentum_compaction(branch)

tao/code/tao_lattice_calc.f90

@@ -151,8 +151,8 @@ subroutine tao_lattice_calc (calc_ok, print_err)
       tao_lat%chrom_calc_ok = .false.
       if (s%com%force_chrom_calc .or. u%calc%chrom_for_data .or. u%calc%chrom_for_plotting) then
         call chrom_calc (tao_lat%lat, s%global%delta_e_chrom, tao_branch%a%chrom, tao_branch%b%chrom, err, &
-                tao_branch%orbit(0)%vec(6), low_E_lat=tao_lat%low_E_lat, high_E_lat=tao_lat%high_E_lat, &
-                ix_branch = ib, orb0 = tao_branch%orbit(0))
+                tao_branch%orbit(0)%vec(6), tao_lat%low_E_lat, tao_lat%high_E_lat, &
+                tao_branch%low_E_orb, tao_branch%high_E_orb, ib, orb0 = tao_branch%orbit(0))
         tao_lat%chrom_calc_ok = (.not. err)
       endif
 

tao/code/tao_show_this.f90

@@ -5826,7 +5826,7 @@ subroutine tao_show_this (what, result_id, lines, nl)
         'Note2: Instability with respect to radiation excitations does not affect such things as the closed orbit calculation.')
     endif
     call chrom_calc (lat, s%global%delta_e_chrom, tao_branch%a%chrom, tao_branch%b%chrom, &
-                          pz = tao_branch%orbit(0)%vec(6), ix_branch = ix_branch)
+                                              pz = tao_branch%orbit(0)%vec(6), ix_branch = ix_branch)
 
     mode_d => design_tao_branch%modes_6d
     mode_m => tao_branch%modes_6d

tao/code/tao_struct.f90

@@ -936,6 +936,7 @@ module tao_struct
   type (normal_modes_struct) modes_6d                     ! 6D radiation matrices.
   type (ptc_normal_form_struct) ptc_normal_form
   type (bmad_normal_form_struct) bmad_normal_form
+  type (coord_struct), allocatable :: high_E_orb(:), low_E_orb(:)
   real(rp) :: cache_x_min = 0, cache_x_max = 0
   real(rp) :: comb_ds_save = -1                           ! Master parameter for %bunch_params_comb(:)%ds_save
   real(rp) :: comb_max_ds_save = -1                       ! Master parameter for %bunch_params_comb(:)%max_ds_save

tao/code/tao_write_cmd.f90

@@ -429,66 +429,8 @@ subroutine tao_write_cmd (what)
 
 case ('curve')
 
-  call tao_find_plots (err, word(1), 'BOTH', curve = curve, blank_means_all = .true.)
-  if (err .or. size(curve) == 0) then
-    call out_io (s_error$, r_name, 'CANNOT FIND CURVE')
-    return
-  endif
-
-  if (size(curve) > 1) then
-    call out_io (s_error$, r_name, 'MULTIPLE CURVES FIT NAME')
-    return
-  endif
-
-  file_name = 'curve.dat'
-  if (word(2) /= ' ') file_name = word(2)
-  call fullfilename (file_name, file_name)
-
-  c => curve(1)%c
-  ok = .false.
-
-  if (c%g%type == "phase_space") then
-    i_uni = c%ix_universe
-    if (i_uni == 0) i_uni = s%global%default_universe
-    beam => s%u(i_uni)%model_branch(c%ix_branch)%ele(c%ix_ele_ref_track)%beam
-    call file_suffixer (file_name, file_name, 'particle_dat', .true.)
-    open (iu, file = file_name)
-    write (iu, '(a, 6(12x, a))') '  Ix', '  x', 'px', '  y', 'py', '  z', 'pz'
-    do i = 1, size(beam%bunch(1)%particle)
-      write (iu, '(i6, 6es15.7)') i, (beam%bunch(1)%particle(i)%vec(j), j = 1, 6)
-    enddo
-    call out_io (s_info$, r_name, 'Written: ' // file_name)
-    close(iu)
-    ok = .true.
-  endif
-
-  if (allocated(c%x_symb) .and. allocated(c%y_symb)) then
-    call file_suffixer (file_name, file_name, 'symbol_dat', .true.)
-    open (iu, file = file_name)
-    write (iu, '(a, 6(12x, a))') '  Ix', '  x', '  y'
-    do i = 1, size(c%x_symb)
-      write (iu, '(i6, 2es15.7)') i, c%x_symb(i), c%y_symb(i)
-    enddo
-    call out_io (s_info$, r_name, 'Written: ' // file_name)
-    close(iu)
-    ok = .true.
-  endif
-
-  if (allocated(c%x_line) .and. allocated(c%y_line)) then
-    call file_suffixer (file_name, file_name, 'line_dat', .true.)
-    open (iu, file = file_name)
-    write (iu, '(a, 6(12x, a))') '  Ix', '  x', '  y'
-    do i = 1, size(c%x_line)
-      write (iu, '(i6, 2es15.7)') i, c%x_line(i), c%y_line(i)
-    enddo
-    call out_io (s_info$, r_name, 'Written: ' // file_name)
-    close(iu)
-    ok = .true.
-  endif
-
-  if (.not. ok) then
-    call out_io (s_info$, r_name, 'No data found in curve to write')
-  endif
+  call out_io (s_info$, r_name, &
+      '"show curve" command superseded by the more versatile "show -write <file> curve ..." command.')
 
 !---------------------------------------------------
 ! derivative_matrix

tao/doc/command-list.tex

@@ -3835,7 +3835,6 @@ \section{write}\index{commands!write}
   write bmad ...                  ! \sref{s:write.bmad}
   write bunch_comb ...            ! \sref{s:write.bunch.comb}
   write covariance_matrix ...     ! \sref{s:write.covar.matrix}
-  write curve ...                 ! \sref{s:write.curve}
   write derivative_matrix ...     ! \sref{s:write.deriv.matrix}
   write digested ...              ! \sref{s:write.digested}
   write elegant ...               ! \sref{s:write.elegant}
@@ -4032,28 +4031,6 @@ \subsection{write covariance_matrix}
 
 The default file name is \vn{covar.matrix}.
 
-%% write curve  --------------------------------------------------------------
-
-\subsection{write curve}
-\label{s:write.curve}
-
-The \vn{write curve} command writes plot curve data to a file.
-Syntax:
-\begin{example}
-    write curve <curve_name> \{<file_name>\} ! Write the curve data
-\end{example}
-
-\vn{write curve} will produce two or three files:
-\begin{example}
-  <file_name>.symbol_dat    ! Symbol coordinates file
-  <file_name>.line_dat      ! Curve coords.
-  <file_name>.particle_dat  ! Particle data file
-\end{example}
-The particle data file is only produced if particle data is associated with the curve.  The curve
-coordinates are the the set of points that are used to draw the (possibly smooth) curve through the
-symbols.
-
-
 %% write derivative_matrix --------------------------------------------------------------
 
 \subsection{write derivative_matrix}

tao/doc/cover-page.tex

@@ -2,7 +2,7 @@
 
 \begin{flushright}
 \large
-Revision: September 28, 2024 \\
+Revision: October 8, 2024 \\
 \end{flushright}
 
 \vfill

tao/doc/data.tex

@@ -1060,6 +1060,7 @@ \section{Tao Data Types}\index{data!data Types}
 
 Chromaticities will be calculated even if the geometry of the lattice branch has an open geometry.
 In this case, dbeta/dpz and dalpha/dpz at the beginning of the branch will be assumed to be zero.
+And $p_z$ is the local $p_z$ at the evaluation point (as opposed to the $p_z$ at the start of the lattice).
 
   %----------------------
   \item[chrom_ptc.a.$N$, chrom_ptc.b.$N$, $N = 0, 1, 2, \ldots$] \Newline \hlabel{chrom.ptc}

tao/version/tao_version_mod.f90

@@ -6,5 +6,5 @@
 !-
 
 module tao_version_mod
-character(*), parameter :: tao_version_date = "2024/09/16 12:29:21"
+character(*), parameter :: tao_version_date = "2024/10/02 00:49:23"
 end module