More lcavity devel. (bmad-sim#1600)

* More lcavity devel.

Author: DavidSagan

bmad/code/em_field_calc.f90

@@ -322,7 +322,7 @@ recursive subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, f
       if (present(rf_time)) then
         time = rf_time
       else
-        time = particle_rf_time(orbit, ele, .false., s_body)
+        time = particle_rf_time(orbit, ele, .false., s_body, rf_freq = ele%value(rf_frequency$))
       endif
       phase = twopi * (ele%value(phi0$) + ele%value(phi0_multipass$) + ele%value(phi0_autoscale$) - &
                       (time - rf_ref_time_offset(ele) - s_body/c_light) * ele%value(rf_frequency$))
@@ -347,7 +347,7 @@ recursive subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, f
       if (present(rf_time)) then
         time = rf_time
       else
-        time = particle_rf_time(orbit, ele, .true., s_body)
+        time = particle_rf_time(orbit, ele, .true., s_body, rf_freq = ele%value(rf_frequency$))
       endif
       phase = (ele%value(phi0$) + ele%value(phi0_multipass$) + ele%value(phi0_err$) + ele%value(phi0_autoscale$))
       field%e(3) = e_accel_field (ele, gradient$) * cos(twopi * (time * ele%value(rf_frequency$) + phase)) / ref_charge
@@ -425,7 +425,7 @@ recursive subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, f
     if (present(rf_time)) then
       time = rf_time
     else
-      time = particle_rf_time(orbit, ele, .true., s_body)
+      time = particle_rf_time(orbit, ele, .true., s_body, rf_freq = ele%value(rf_frequency$))
     endif
 
     if (nint(ele%value(cavity_type$)) == traveling_wave$) then
@@ -782,13 +782,6 @@ recursive subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, f
 ! FieldMap
 
 case(fieldmap$)
-
-  if (present(rf_time)) then
-    time = rf_time
-  else
-    time = particle_rf_time(orbit, ele, .false., s_body)
-  endif
-
   if (.not. associated(ele%cylindrical_map) .and. .not. associated(ele%cartesian_map) .and. &
       .not. associated(ele%gen_grad_map) .and. .not. associated(ele%grid_field)) then
     call out_io (s_fatal$, r_name, 'No associated fieldmap (cartesican_map, grid_field, etc) FOR: ' // ele%name) 
@@ -1108,6 +1101,12 @@ recursive subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, f
         freq0 = ele%value(rf_frequency$)
         freq = ele%value(rf_frequency$) * cl_map%harmonic
 
+        if (present(rf_time)) then
+          time = rf_time
+        else
+          time = particle_rf_time(orbit, ele, .false., s_body, rf_freq = freq0)
+        endif
+
         if (freq0 == 0) then
           call out_io (s_fatal$, r_name, 'Element frequency is zero but cylindrical_map harmonic is not in: ' // ele%name)
           if (global_com%exit_on_error) call err_exit
@@ -1314,6 +1313,12 @@ recursive subroutine em_field_calc (ele, param, s_pos, orbit, local_ref_frame, f
           return  
         endif
 
+        if (present(rf_time)) then
+          time = rf_time
+        else
+          time = particle_rf_time(orbit, ele, .false., s_body, rf_freq = freq0)
+        endif
+
         t_ref = (ele%value(phi0$) + ele%value(phi0_multipass$) + ele%value(phi0_err$) + &
                                                   phi0_autoscale + g_field%phi0_fieldmap) / freq0
         if (ele%key == rfcavity$) t_ref = 0.25/freq0 - t_ref

bmad/code/particle_rf_time.f90

@@ -1,5 +1,5 @@
 !+
-! Function particle_rf_time (orbit, ele, reference_active_edge, s_rel, time_coords, rf_freq, rf_clock_harmonic, abs_time) result (time)
+! Function particle_rf_time (orbit, ele, reference_active_edge, s_rel, time_coords, rf_freq, abs_time) result (time)
 !
 ! Routine to return the reference time used to calculate the phase of
 ! time-dependent EM fields.
@@ -19,17 +19,18 @@
 !   s_rel             -- real(rp), optional: Longitudinal position relative to the upstream edge of the element.
 !                         Needed for relative time tracking when the particle is inside the element. Default is 0.
 !   time_coords       -- logical, optional: Default False. If True then orbit is using time based phase space coordinates.
-!   rf_freq           -- real(rp), optional: If present and non-zero, the returned time shifted by an integer multiple
-!                          of 1/rf_freq to be in the range [-1/2*rf_freq, 1/2*rf_freq].
 !   rf_clock_harmonic -- integer, optional: Used with the rf clock in cases where an element has multiple frequencies.
+!   rf_freq           -- real(rp), optional: If present, the returned time shifted by an integer multiple
+!                          of 1/rf_freq to be in the range [-1/2*rf_freq, 1/2*rf_freq]. This is useful to
+!                          avoid round-off errors.
 !   abs_time          -- real(rp), optional: If False (default) use setting of bmad_com%absolute_time_tracking.
 !                           If True, use absolute time instead of relative time.
 !
 ! Ouput:
 !   time      -- Real(rp): Current time.
 !-
 
-function particle_rf_time (orbit, ele, reference_active_edge, s_rel, time_coords, rf_freq, rf_clock_harmonic, abs_time) result (time)
+function particle_rf_time (orbit, ele, reference_active_edge, s_rel, time_coords, rf_freq, abs_time) result (time)
 
 use bmad_routine_interface, dummy_except => particle_rf_time
 use attribute_mod, only: has_attribute
@@ -42,8 +43,8 @@ function particle_rf_time (orbit, ele, reference_active_edge, s_rel, time_coords
 type (ele_pointer_struct), allocatable :: chain(:)
 
 real(rp), optional :: s_rel, rf_freq
-real(rp) time, s_hard_offset, beta0, freq
-integer, optional :: rf_clock_harmonic
+real(qp) time
+real(rp) s_hard_offset, beta0, freq
 integer n, ix_pass, n_links, harmonic
 logical, optional :: reference_active_edge, time_coords, abs_time
 
@@ -113,10 +114,9 @@ function particle_rf_time (orbit, ele, reference_active_edge, s_rel, time_coords
   time = time - s_hard_offset / (c_light * beta0)
 endif
 
-!
-
-freq = real_option(0.0_rp, rf_freq)
-if (freq /= 0) time = modulo2(time, 0.5_rp/freq)
+if (present(rf_freq)) then
+  if (rf_freq /= 0) time = modulo2(time, 0.5_qp/rf_freq)
+endif
 
 end function particle_rf_time
 

bmad/code/track1.f90

@@ -226,14 +226,17 @@ recursive subroutine track1 (start_orb, ele, param, end_orb, track, err_flag, ig
 case (bmad_standard$)
   if (end_orb%species == photon$) then
     call track1_bmad_photon (end_orb, ele, param, err)
+    if (err) return
+
   else
     call track1_bmad (end_orb, ele, param, err, track, mat6 = ele%mat6, make_matrix = make_map1)
+    if (err) return
 
     select case (ele%key)
     case (beambeam$, crab_cavity$, patch$);  do_spin_tracking = .false.
+    case (lcavity$);                         do_spin_tracking = (do_spin_tracking .and. nint(ele%value(n_rf_steps$)) < 1)
     end select
   endif
-  if (err) return
 
 case (runge_kutta$, fixed_step_runge_kutta$) 
   call track1_runge_kutta (end_orb, ele, param, err, track, mat6 = ele%mat6, make_matrix = make_map1)

bmad/doc/charged-tracking.tex

@@ -1172,7 +1172,8 @@ \subsection{drift-kick model}
 kicks are applied at the ends of the element and in between sections with the kicks at the ends
 being half that of the kicks applied in the interior. The voltage at a kick point is:
 \begin{equation}
-  V_{kick} = \frac{r_q \, \kappa \, V_{tot}}{N_s} \, \cos(2 \pi \, (\phi_t + \phi_\REF))
+  V_txt{kick} = \frac{r_q \, \kappa \, V_{tot}}{N_s} \, \cos(2 \pi \, (\phi_t + \phi_\REF))
+  \label{vkick}
 \end{equation}
 See \Eq{dergl} for a description of the parameters in the equation. Here $\kappa$ is 0.5 for the end
 kick points and is unity for interior kick points.
@@ -1193,6 +1194,9 @@ \subsection{drift-kick model}
 in each section will affect the plotting of phase space momentum coordinates since these
 are normalized to $P_0$ (\sref{s:phase.space}).
 
+For spin tracking through the energy kick, the BMT equation is used with the integrated kick field
+given by Eq.~{vkick} modeled as an electric field in the $z$-direction.
+
 %---------------------------------------------------------------------------------
 \subsection{Edge kick}
 
@@ -1205,7 +1209,7 @@ \subsection{Edge kick}
   H_f = \mp \frac{q}{2 \, P_0 \, c} G_t \, \cos(\phi_t + \phi_\REF) \, (x^2 + y^2)
 \end{equation}
 where the minus sign is for the entrance end and the plus sign is for the exit end, $q$ is the
-particle charge, $P_0$ is the reference momentum.  Since $\phi_t = 2 \pi f_\txt{rf} \, t$ is time
+particle charge, $P_0$ is the reference momentum. Since $\phi_t = 2 \pi f_\txt{rf} \, t$ is time
 dependent, it is easiest to use this Hamiltonian in energy-time phase space coordinates
 (\sref{energy.phase.space}) to give the fringe kicks:
 \begin{align}
@@ -1223,7 +1227,8 @@ \subsection{Edge kick}
 oscillation also delays the particle in time compared to the case where there is no oscillation
 conponent.
 
-For spin tracking, the BMT equation is used with the integrated edge field (voltage) $\bfV_f$:
+For spin tracking through the edge field, the BMT equation is used with the integrated edge field
+(voltage) $\bfV_f$ parallel to the radius vector $\bfr$:
 \begin{equation}
   \bfV_f = -G_t \, \bfr
 \end{equation}
@@ -1267,10 +1272,12 @@ \subsection{Standing wave focusing: the pondermotive force}
 forward propagating particle is rapidly varying. The formulas above are for the kick averaged over
 all phases and therefore the phase dependence has been eliminated.
 
-For spin tracking, the focusing kick is modeled as a radial integrated field (voltage) $\bfV_f$
-\begin{equation}
-  \bfV_f = -G_t \, \bfr
-\end{equation}
+For spin tracking, to zeroth order, the averaged field is zero for the backwards wave so no spin
+kick is applied for the backwards wave.
+\begin{align}
+  \bfE(\bfr)   &= \CRNO
+  \bfB(\theta) &= 
+\end{align}
 
 %-----------------------------------------------------------------------------------------
 %-----------------------------------------------------------------------------------------

bmad/doc/methods.tex

@@ -528,14 +528,14 @@ \section{Spin Tracking Methods}
 \begin{sideways}\vn{Magnus}\end{sideways} &
 \begin{sideways}\vn{Sprint}\end{sideways} &
 \begin{sideways}\vn{Symp_Lie_PTC}\end{sideways} &
-\begin{sideways}\vn{Tracking}\end{sideways}
+\begin{sideways}\vn{Tracking}\end{sideways} &
 \begin{sideways}\vn{Transverse_Kick}\end{sideways} &
 \\ \midrule
 %                                      O   C   M  Spt SLP Trk TK
   \vn{ab_multipole and multipole}    & X & X &   &   &   & D &   \\  
   \vn{ac_kicker}                     & X & X &   &   &   & D &   \\  
   \vn{beambeam}                      & X & X &   &   &   & D &   \\  
-  \vn{bends: rbend and sbend}        & X & X & X  & X & X & D &   \\ 
+  \vn{bends: rbend and sbend}        & X & X & X & X & X & D &   \\ 
   \vn{converter}                     & X & X &   &   &   & D &   \\  
   \vn{crab_cavity}                   & X & X &   &   &   & D &   \\  
   \vn{custom}                        & X & D &   &   &   & X &   \\  
@@ -547,10 +547,10 @@ \section{Spin Tracking Methods}
   \vn{fiducial}                      & X & X &   &   & X & D &   \\  
   \vn{floor_shift}                   & X & X &   &   & X & D &   \\  
   \vn{hkicker}                       & X & X & X & X & X & D &   \\  
-  \vn{instrument, monitor and pipe}  & X & X & X  & X & X & D &   \\  
+  \vn{instrument, monitor and pipe}  & X & X & X & X & X & D &   \\  
   \vn{kicker}                        & X & X & X & X & X & D &   \\  
-  \vn{lcavity and rfcavity}          & X & X & X &   & X & D &   \\  
-  \vn{marker}                        & X & X & X &   & X & D &   \\  
+  \vn{lcavity and rfcavity}          & X & X &   &   & X & D &   \\  
+  \vn{marker}                        & X & X &   &   &   & D &   \\  
   \vn{match}                         & D &   &   &   &   &   & X \\  
   \vn{octupole}                      & X & X &   & X & X & D &   \\ 
   \vn{patch}                         & X & X &   &   & X & D &   \\ 

bmad/low_level/ac_kicker_amp.f90

@@ -37,7 +37,8 @@ function ac_kicker_amp(ele, orbit, true_time) result (ac_amp)
 !
 
 ref_ele => ele
-if (ref_ele%slave_status == super_slave$ .or. ele%slave_status == slice_slave$) ref_ele => pointer_to_super_lord (ref_ele)
+if (ref_ele%slave_status == super_slave$ .or. ele%slave_status == slice_slave$) &
+                                                                      ref_ele => pointer_to_super_lord (ref_ele)
 
 ac_amp = 1
 if (ele%key /= ac_kicker$) return
@@ -55,12 +56,12 @@ function ac_kicker_amp(ele, orbit, true_time) result (ac_amp)
 if (allocated(ac%frequency)) then
   ac_amp = 0
   do i = 1, size(ac%frequency)
-    t = real_option(particle_rf_time(orbit, ele, rf_clock_harmonic = ac%frequency(i)%rf_clock_harmonic), true_time)
+    t = real_option(real(particle_rf_time(orbit, ele, rf_freq = ac%frequency(i)%f), rp), true_time)
     ac_amp = ac_amp + ac%frequency(i)%amp * cos(twopi*(ac%frequency(i)%f * t + ac%frequency(i)%phi))
   enddo
 
 else
-  t = real_option(particle_rf_time(orbit, ele), true_time)
+  t = real_option(real(particle_rf_time(orbit, ele), rp), true_time)
   ac_amp = knot_interpolate(ac%amp_vs_time%time, ac%amp_vs_time%amp, t, nint(ele%value(interpolation$)), err_flag)
   if (err_flag) then
     call out_io (s_fatal$, r_name, 'INTERPOLATION PROBLEM FOR AC_KICKER: ' // ele%name)

bmad/low_level/track_a_crab_cavity.f90

@@ -81,7 +81,7 @@ subroutine track_a_crab_cavity (orbit, ele, param, mat6, make_matrix)
   s_here = (i - 0.5_rp) * dl
 
   phase0 = twopi * (ele%value(phi0$) + ele%value(phi0_multipass$) + ele%value(phi0_autoscale$) - &
-          (particle_rf_time (orbit, ele, .false.) - rf_ref_time_offset(ele, s_here)) * &
+          (particle_rf_time (orbit, ele, .false., rf_freq = ele%value(rf_frequency$)) - rf_ref_time_offset(ele, s_here)) * &
           ele%value(rf_frequency$))
   if (ele%orientation == -1) phase0 = phase0 + twopi * ele%value(rf_frequency$) * dt_length
   phase = phase0