New write field command.

bmad/code/em_field_calc.f90

@@ -1,5 +1,5 @@
 !+
-! Subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, field, calc_dfield, err_flag, &
+! Subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, field, calc_dfield, err_flag,
 !               calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles, print_err)
 !
 ! Routine to calculate the E and B fields at a particular place in an element.
@@ -653,9 +653,9 @@ recursive subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, f
       do i = 0, ix_pole_max
         if (a_pole(i) == 0 .and. b_pole(i) == 0) cycle
         if (do_df_calc) then
-          call ab_multipole_kick(a_pole(i), b_pole(i), i, local_orb%species, 0, local_orb, kx, ky, dkm)
+          call ab_multipole_kick(a_pole(i), b_pole(i), i, ele%ref_species, 0, local_orb, kx, ky, dkm)
         else
-          call ab_multipole_kick(a_pole(i), b_pole(i), i, local_orb%species, 0, local_orb, kx, ky)
+          call ab_multipole_kick(a_pole(i), b_pole(i), i, ele%ref_species, 0, local_orb, kx, ky)
         endif
         field%B(1) = field%B(1) + f_p0c * ky / ele%value(l$)
         field%B(2) = field%B(2) - f_p0c * kx / ele%value(l$)
@@ -1378,7 +1378,7 @@ recursive subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, f
 
         if (logic_option(.false., calc_potential)) then
           if (r /= 0) then
-            abs_tol = abs(1e-10_rp * r * orbit%p0c * (1 + orbit%vec(6)) / (c_light * charge_of(orbit%species)))
+            abs_tol = abs(1e-10_rp * r * orbit%p0c * (1 + orbit%vec(6)) / (c_light * charge_of(ele%ref_species)))
             inte = super_qromb(rb_field, 0.0_rp, r, 1e-12_rp, abs_tol, 2, err) / r
             field%A(1:2) = field%A(1:2) + inte * [-y, x] / r
           endif

bmad/code/offset_particle.f90

@@ -154,10 +154,10 @@ subroutine offset_particle (ele, set, orbit, set_tilt, set_hvkicks, drift_to_edg
 
 rel_p = (1 + orbit%vec(6))
 
-rel_tracking_charge = rel_tracking_charge_to_mass (orbit, ele%ref_species)
-charge_dir = rel_tracking_charge * sign_z_vel * orbit%time_dir 
-
 if (set_hv) then
+  rel_tracking_charge = rel_tracking_charge_to_mass (orbit, ele%ref_species)
+  charge_dir = rel_tracking_charge * sign_z_vel * orbit%time_dir 
+
   select case (ele%key)
   case (elseparator$, kicker$, hkicker$, vkicker$)
     set_hv1 = .false.

tao/code/tao_locate_all_elements.f90

@@ -8,7 +8,7 @@
 ! Input:
 !   ele_list     -- Character(*): String with element names using element list format.
 !   ignore_blank -- Logical, optional: If present and true then do nothing if
-!     ele_list is blank. otherwise treated as an error.
+!     ele_list is blank. otherwise a blank is treated as an error.
 !
 ! Output:
 !   eles  -- ele_pointer_struct(:), allocatable: Array of elements in the model lat. 

tao/code/tao_write_cmd.f90

@@ -33,7 +33,9 @@ subroutine tao_write_cmd (what)
 type (branch_struct), pointer :: branch
 type (ele_pointer_struct), allocatable :: eles(:)
 type (ele_struct), pointer :: ele
+type (lat_param_struct), pointer :: param
 type (coord_struct), pointer :: p
+type (coord_struct)  :: orbit
 type (tao_d2_data_struct), pointer :: d2
 type (tao_d1_data_struct), pointer :: d1
 type (tao_data_struct), pointer :: dat
@@ -43,6 +45,7 @@ subroutine tao_write_cmd (what)
 type (bunch_track_struct), pointer :: bunch_params_comb(:)
 type (bunch_track_struct), pointer :: comb1
 type (bunch_params_struct), pointer :: bpt
+type (em_field_struct) field
 
 real(rp) scale, mat6(6,6)
 real(rp), allocatable :: values(:)
@@ -57,7 +60,7 @@ subroutine tao_write_cmd (what)
 character(200) :: word(40)
 character(*), parameter :: r_name = 'tao_write_cmd'
 
-real(rp) dr(3), r_max(3), r_min(3)
+real(rp) dr(3), r_max(3), r_min(3), rr(3)
 
 integer i, j, k, m, n, ie, ic, id, ix, iu, nd, ii, i_uni, ib, ip, ios, loc, iy, iz
 integer n_max(3), n_min(3), i_chan, ix_beam, ix_word, ix_w2, file_format
@@ -519,6 +522,7 @@ subroutine tao_write_cmd (what)
   n_min = int_garbage$
   n_max = int_garbage$
   file_name = ''
+  ix_word = 0
 
   do
     ix_word = ix_word + 1
@@ -557,7 +561,7 @@ subroutine tao_write_cmd (what)
   endif
 
   if (r_max(1) /= real_garbage$ .and. r_min(1) == real_garbage$) r_min = [-r_max(1), -r_max(2), 0.0_rp]
-  if (n_max(1) /= real_garbage$ .and. n_min(1) == real_garbage$) n_min = [-n_max(1), -n_max(2), 0]
+  if (n_max(1) /= int_garbage$ .and. n_min(1) == int_garbage$) n_min = [-n_max(1), -n_max(2), 0]
 
   if (r_max(1) == real_garbage$) then
     r_max = n_max * dr
@@ -569,19 +573,43 @@ subroutine tao_write_cmd (what)
     dr = (r_max - r_min) / (n_max - n_min)
   endif
 
-  call tao_locate_elements(ele_name, -2, eles, err)
+  call tao_locate_elements (ele_name, s%global%default_universe, eles, err)
   if (err) return
 
   if (size(eles) > 1) then
-    call out_io(s_error$, r_name, 'ELEMENT NAME MATCHES TO MULTIPLE ELEMENTS: ' // ele_name)
-    return
+    call out_io(s_warn$, r_name, 'Element name matches to multiple elements: ' // ele_name, 'Will use first one.')
   elseif (size(eles) == 0) then
     call out_io(s_error$, r_name, 'ELEMENT NAME DOES NOT MATCH TO ANY ELEMENTS: ' // ele_name)
     return
   endif
 
-  do 
+  ele => eles(1)%ele
+
+  !
 
+  open (iu, file = file_name)
+  write (iu, '(9a)') '# ele = ', quote(ele%name) 
+  write (iu, '(9a)') '# nx = [', int_str(n_min(1)), ', ', int_str(n_max(1)), ']' 
+  write (iu, '(9a)') '# ny = [', int_str(n_min(2)), ', ', int_str(n_max(2)), ']' 
+  write (iu, '(9a)') '# nz = [', int_str(n_min(3)), ', ', int_str(n_max(3)), ']' 
+  write (iu, '(9a)') '# rx = [', real_str(r_min(1)), ', ', real_str(r_max(1)), ']' 
+  write (iu, '(9a)') '# ry = [', real_str(r_min(2)), ', ', real_str(r_max(2)), ']' 
+  write (iu, '(9a)') '# rz = [', real_str(r_min(3)), ', ', real_str(r_max(3)), ']' 
+  write (iu, '(9a)') '# dr = [', real_str(dr(1)), ', ', real_str(dr(2)), ', ', real_str(dr(3)), ']'
+  write (iu, '(9a)') '##' 
+  write (iu, '(9a)') '##         x           y           z                  Bx                  By                  Bz' 
+
+  branch => pointer_to_branch(ele)
+
+  do ix = n_min(1), n_max(1)
+  do iy = n_min(2), n_max(2)
+  do iz = n_min(3), n_max(3)
+    rr = [ix, iy, iz] * dr
+    orbit%vec(1:3:2) = rr(1:2)
+    call em_field_calc(ele, branch%param, rr(3), orbit, .false., field)
+    write (iu, '(3f12.6, 3es20.11)') rr, field%b
+  enddo
+  enddo
   enddo
 
 !---------------------------------------------------

tao/doc/command-list.tex

@@ -3846,6 +3846,7 @@ \section{write}\index{commands!write}
   write derivative_matrix ...     ! \sref{s:write.deriv.matrix}
   write digested ...              ! \sref{s:write.digested}
   write elegant ...               ! \sref{s:write.elegant}
+  write field ...                 ! \sref{s:write.field
   write gif ...                   ! \sref{s:write.gif}
   write hard                      ! \sref{s:write.hard}
   write mad8 ...                  ! \sref{s:write.mad8}
@@ -4077,6 +4078,48 @@ \subsection{write elegant}
 Write a lattice file in \vn{Elegant} format. The default file name is \vn{lat_\#.lte} where \vn{\#}
 is replaced by the universe number.
 
+%% write field --------------------------------------------------------------
+
+\subsection{write field}
+\label{s:write.field}
+
+The \vn{write field} command creates a file with a table of magnetic field values for a given lattice
+element. Syntax:
+\begin{example}
+    write field \{-nmax <nx_max>, <ny_max> <nz_max>\} \{-nmin <nx_min>, <ny_min> <nz_min>\}
+                \{-rmax <rx_max>, <ry_max> <rz_max>\} \{-rmin <rx_min>, <ry_min> <rz_min>\} 
+        \{-dr <dr_x> <dr_y <dr_z>\} -ele <eleID> \{<file_name>\}  
+\end{example}
+
+The default file name is \vn{field.dat}.
+
+The fiducial point for the grid is the coordinates at the entrance end of the the lattice element
+given by the \vn{-ele} switch. Laboratory (not body) coordinates are used.
+
+The field grid indexes runs from \vn{nx_min} to \vn{nx_max} for $x$,
+etc. The distance between points is set with the \vn{-dr} switch and the extent of the grid is set
+by \vn{-rmin} and \vn{rmax} switches. 
+
+If \vn{-nmin} is present, so must \vn{-nmax} be present. If \vn{-nmin} is not present but \vn{-nmax}
+is, the default values for \vn{-nmin} are [-nx_max, -ny_max, 0].
+
+If \vn{-rmin} is present, so must \vn{-rmax} be present. If \vn{-rmin} is not present but \vn{-rmax}
+is, the default values for \vn{-rmin} are [-rx_max, -ry_max, 0].
+
+\vn{-dr}, \vn{-nmax}, and \vn{-rmax} values are interdependent. For example:
+\begin{example}
+  rx_max = nx_max * dr_x
+\end{example}
+Given this, exactly two of the three needs to be present.
+
+Examples:
+\begin{example}
+  write field -ele Q1 -dr 0.01 0.02 0.05 -nmax 30 20, 200
+  write field -ele Q1 -dr 0.01 0.02 0.05 -rmax 0.3 0.4 10.0   ! Same as above
+\end{example}
+The above examples both specify the same grid which will have an index range in $x$ of $[-30, 30]$,
+ana range in $y$ of $[-20, 20]$, 
+
 %% write gif --------------------------------------------------------------
 
 \subsection{write gif}