Better spin tracking method bookkeeping and added spin_taylor with fi…

…nite differences. (#1320)

* Better spin tracking method bookkeeping and added spin_taylor with finite differences.

bmad/code/init_taylor_series.f90

@@ -0,0 +1,60 @@
+!+
+! Subroutine init_taylor_series (bmad_taylor, n_term, save_old)
+!
+! Subroutine to initialize or extend a Bmad Taylor series (6 of these series make
+! a Taylor map). Note: This routine does not zero the terms.
+!
+! Input:
+!   bmad_taylor -- taylor_struct: Old structure.
+!   n_term      -- integer: Number of terms to allocate. 
+!                    n_term < 0 => bmad_taylor%term pointer will be disassociated.
+!   save_old    -- logical, optional: If True then save any old terms and ref orbit when
+!                    bmad_taylor is resized. If False zero the ref orbit. Default is False.
+!
+! Output:
+!   bmad_taylor -- Taylor_struct: Initalized structure.
+!-
+
+subroutine init_taylor_series (bmad_taylor, n_term, save_old)
+
+use bmad_routine_interface, dummy => init_taylor_series
+
+implicit none
+
+type (taylor_struct) bmad_taylor
+type (taylor_term_struct), pointer :: term(:)
+integer n_term
+integer n
+logical, optional :: save_old
+
+!
+
+if (.not. logic_option (.false., save_old)) bmad_taylor%ref = 0
+
+if (n_term < 0) then
+  if (associated(bmad_taylor%term)) deallocate(bmad_taylor%term)
+  return
+endif
+
+if (.not. associated (bmad_taylor%term)) then
+  allocate (bmad_taylor%term(n_term))
+  return
+endif
+
+if (size(bmad_taylor%term) == n_term) return
+
+!
+
+if (logic_option (.false., save_old) .and. n_term > 0 .and. size(bmad_taylor%term) > 0) then
+  n = min (n_term, size(bmad_taylor%term))
+  term => bmad_taylor%term
+  allocate (bmad_taylor%term(n_term))
+  bmad_taylor%term(1:n) = term(1:n)
+  deallocate (term)
+
+else
+  deallocate (bmad_taylor%term)
+  allocate (bmad_taylor%term(n_term))
+endif
+
+end subroutine init_taylor_series

bmad/code/make_hybrid_lat.f90

@@ -268,7 +268,7 @@ subroutine make_hybrid_lat (lat_in, lat_out, use_taylor, orb0_arr)
             c0%vec = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
           endif
           if (.not. associated(ele_out%taylor(1)%term)) then ! construct taylor
-            call ele_to_taylor (ele_out, b_in%param, c0)
+            call ele_to_taylor (ele_out, c0)
           endif
         endif
 

bmad/code/make_mat6.f90

@@ -132,7 +132,7 @@ recursive subroutine make_mat6 (ele, param, start_orb, end_orb, err_flag)
   call make_mat6_custom_ptr (ele, param, a_start_orb, a_end_orb, err)
 
 case (taylor$)
-  call make_mat6_taylor (ele, param, a_start_orb, a_end_orb, err)
+  call make_mat6_taylor (ele, a_start_orb, a_end_orb, err)
 
 case (bmad_standard$)
   if (a_start_orb%species == photon$) then
@@ -142,7 +142,7 @@ recursive subroutine make_mat6 (ele, param, start_orb, end_orb, err_flag)
   endif
 
 case (symp_lie_ptc$)
-  call make_mat6_symp_lie_ptc (ele, param, a_start_orb, a_end_orb)
+  call make_mat6_symp_lie_ptc (ele, a_start_orb, a_end_orb)
 
 case (symp_lie_bmad$)
   a_end_orb = a_start_orb
@@ -194,25 +194,29 @@ recursive subroutine make_mat6 (ele, param, start_orb, end_orb, err_flag)
   endif
 endif
 
-! Spin
 
-if (bmad_com%spin_tracking_on) then
-  if (ele%spin_tracking_method == sprint$ .and. .not. associated(ele%spin_taylor(0)%term)) then
-    call sprint_spin_taylor_map(ele, ele%taylor%ref)
-  endif
+! Finish up mat6 calc
 
-  if (associated(ele%spin_taylor(0)%term)) then
-    ele%spin_q = spin_taylor_to_linear (ele%spin_taylor, .true., a_start_orb%vec-ele%spin_taylor_ref_orb_in, ele%is_on)
-  endif
-endif
+ele%map_ref_orb_out = a_end_orb
+if (ele%bookkeeping_state%mat6 == stale$) ele%bookkeeping_state%mat6 = ok$
 
-! Finish up
+! Spin
 
-ele%map_ref_orb_out = a_end_orb
+!if (bmad_com%spin_tracking_on .and. a_start_orb%species /= photon$) then
+!  if (.not. associated(ele%spin_taylor(0)%term)) then
+!    call ele_to_spin_taylor(ele, param, a_start_orb)
+!  endif
+!
+!  if (associated(ele%spin_taylor(0)%term)) then
+!    ele%spin_q = spin_taylor_to_linear (ele%spin_taylor, .true., a_start_orb%vec-ele%spin_taylor_ref_orb_in, ele%is_on)
+!  else
+!    call out_io (s_error$, r_name, 'CANNOT CONSTRUCT SPIN TAYLOR MAP FOR ELEMENT: ' // ele_full_name(ele))
+!  endif
+!endif
 
-bmad_com%radiation_fluctuations_on = rad_fluct_save
+! And end
 
-if (ele%bookkeeping_state%mat6 == stale$) ele%bookkeeping_state%mat6 = ok$
+bmad_com%radiation_fluctuations_on = rad_fluct_save
 if (present(err_flag)) err_flag = .false.
 
 end subroutine

bmad/code/track1.f90

@@ -243,7 +243,7 @@ recursive subroutine track1 (start_orb, ele, param, end_orb, track, err_flag, ig
   call track1_linear (end_orb, ele, param)
 
 case (taylor$) 
-  call track1_taylor (end_orb, ele, param)
+  call track1_taylor (end_orb, ele)
 
 case (symp_lie_bmad$) 
   call symp_lie_bmad (ele, param, end_orb, track, mat6 = ele%mat6, make_matrix = make_map1)

bmad/code/valid_tracking_method.f90

@@ -60,7 +60,10 @@ function valid_tracking_method (ele, species, tracking_method) result (is_valid)
 case (ac_kicker$)
   select case (method)
   case (bmad_standard$, runge_kutta$, time_runge_kutta$, linear$, custom$, symp_lie_ptc$, taylor$)
-    is_valid = ((method /= symp_lie_ptc$ .and. method /= taylor$) .or. allocated(ele%ac_kick%frequency))
+    is_valid = (method /= symp_lie_ptc$ .and. method /= taylor$)
+    if (.not. is_valid .and. allocated(ele%ac_kick%frequency)) then
+      if (size(ele%ac_kick%frequency) == 1) is_valid = .true.
+    endif
   end select
 
 case (beambeam$)

bmad/low_level/make_mat6_symp_lie_ptc.f90

@@ -1,11 +1,10 @@
 !+
-! Subroutine make_mat6_symp_lie_ptc (ele, param, start_orb, end_orb)
+! Subroutine make_mat6_symp_lie_ptc (ele, start_orb, end_orb)
 !
 ! Subroutine to make the 6x6 transfer matrix for an element. 
 !
 ! Input:
 !   ele        -- Ele_struct: Element with transfer matrix
-!   param      -- lat_param_struct: Parameters are needed for some elements.
 !   start_orb  -- Coord_struct: Coordinates at the beginning of element. 
 !
 ! Output:
@@ -15,15 +14,14 @@
 !   end_orb    -- Coord_struct: Coordinates at end of element.
 !-
 
-subroutine make_mat6_symp_lie_ptc (ele, param, start_orb, end_orb)
+subroutine make_mat6_symp_lie_ptc (ele, start_orb, end_orb)
 
 use ptc_interface_mod, except_dummy => make_mat6_symp_lie_ptc
 
 implicit none
 
 type (ele_struct), target :: ele
 type (coord_struct) :: start_orb, end_orb, s_orb
-type (lat_param_struct)  param
 type (taylor_struct) bmad_taylor(6), spin_taylor(0:3)
 
 real(rp) dtime_ref, quat(0:3)
@@ -32,7 +30,7 @@ subroutine make_mat6_symp_lie_ptc (ele, param, start_orb, end_orb)
 
 s_orb = start_orb
 
-call ele_to_taylor(ele, param, start_orb, .true., orbital_taylor = bmad_taylor, spin_taylor = spin_taylor)
+call ele_to_taylor(ele, start_orb, .true., orbital_taylor = bmad_taylor, spin_taylor = spin_taylor)
 call taylor_to_mat6 (bmad_taylor, start_orb%vec, ele%vec0, ele%mat6)
 
 end_orb = start_orb

bmad/low_level/make_mat6_taylor.f90

@@ -1,11 +1,10 @@
 !+
-! Subroutine make_mat6_taylor (ele, param, start_orb, end_orb, err_flag)
+! Subroutine make_mat6_taylor (ele, start_orb, end_orb, err_flag)
 !
 ! Subroutine to make the 6x6 transfer matrix for an element. 
 !
 ! Input:
 !   ele      -- Ele_struct: Element to track through.
-!   param     -- lat_param_struct: Parameters are needed for some elements.
 !   start_orb -- coord_struct: Starting coords.
 !
 ! Output:
@@ -16,7 +15,7 @@
 !   err       -- Logical, optional: Set True if there is an error. False otherwise.
 !-
 
-subroutine make_mat6_taylor (ele, param, start_orb, end_orb, err_flag)
+subroutine make_mat6_taylor (ele, start_orb, end_orb, err_flag)
 
 use bmad_interface, dummy => make_mat6_taylor
 use ptc_interface_mod, only: ele_to_taylor
@@ -25,20 +24,19 @@ subroutine make_mat6_taylor (ele, param, start_orb, end_orb, err_flag)
 
 type (ele_struct), target :: ele
 type (coord_struct) :: start_orb, end_orb
-type (lat_param_struct)  param
 
 logical, optional :: err_flag
 
 !
 
 if (present(err_flag)) err_flag = .false.
 
-if (.not. associated(ele%taylor(1)%term)) call ele_to_taylor(ele, param, start_orb)
+if (.not. associated(ele%taylor(1)%term)) call ele_to_taylor(ele, start_orb)
 
 call mat_make_unit (ele%mat6)
 
 end_orb = start_orb
-call track1_taylor (end_orb, ele, param, mat6 = ele%mat6, make_matrix = .true.)
+call track1_taylor (end_orb, ele, mat6 = ele%mat6, make_matrix = .true.)
 
 ele%vec0 = end_orb%vec - matmul(ele%mat6, start_orb%vec)
 

bmad/low_level/make_mat6_tracking.f90

@@ -1,5 +1,5 @@
 !+
-! Subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
+! Subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag, spin_only)
 !
 ! Subroutine to make the 6x6 transfer matrix for an element using the
 ! Present tracking method.
@@ -12,16 +12,18 @@
 !   ele       -- Ele_struct: Element with transfer matrix
 !   param     -- lat_param_struct: Parameters are needed for some elements.
 !   start_orb -- Coord_struct: Coordinates at the beginning of element. 
+!   spin_only -- logical, optional: Default False. If True, only calculate ele%spin_taylor. 
 !
 ! Output:
 !   ele       -- Ele_struct: Element with transfer matrix.
-!     %vec0      -- 0th order map component
-!     %mat6      -- 6x6 transfer matrix.
+!     %vec0         -- 0th order map component
+!     %mat6         -- 6x6 transfer matrix.
+!     %spin_taylor  -- Spin taylor map (if spin tracking is on).
 !   end_orb   -- Coord_struct: Coordinates at the end of element.
 !   err_flag  -- logical: Set True if there is an error. False otherwise.
 !-
 
-subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
+subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag, spin_only)
 
 use bmad_interface, except_dummy => make_mat6_tracking
 
@@ -32,8 +34,11 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
 type (lat_param_struct)  param
 
 real(rp) del_orb(6), dorb6, abs_p, mat6(6,6)
+real(rp) q1(0:3), q2(0:3), q_map(0:3,0:6)
 integer i
-logical err_flag
+logical err_flag, do_spin
+logical, optional :: spin_only
+
 character(*), parameter :: r_name = 'make_mat6_tracking'
 
 ! This computation is singular when start%vec(6) = -1 (zero starting velocity).
@@ -43,6 +48,7 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
 del_orb = bmad_com%d_orb
 abs_p = max(abs(start_orb%vec(2)) + abs(del_orb(2)), abs(start_orb%vec(4)) + abs(del_orb(4)), abs(del_orb(6)))
 bmad_private%random_on = .false.
+do_spin = (logic_option(.false., spin_only) .or. bmad_com%spin_tracking_on)
 
 ! The factor of 1.01 is used to avoid roundoff problems.
 ! Note: init_coord is avoided since init_coord will make z and t consistent with the element's t_ref.
@@ -51,7 +57,7 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
 dorb6 = max(0.0_rp, 1.01 * (abs_p - (1 + start_orb%vec(6))))   ! Shift in start%vec(6) to apply.
 start_orb0 = start_orb
 
-call track1 (start_orb0, ele, param, end_orb)
+call track_this (start_orb0, ele, param, end_orb, q_map(:,0), .false.)
 if (end_orb%state /= alive$) then
   call out_io (s_error$, r_name, 'PARTICLE LOST IN TRACKING CENTRAL PARTICLE. MATRIX NOT CALCULATED FOR ELEMENT: ' // ele%name)
   bmad_private%random_on = .true.
@@ -64,8 +70,7 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
   start = start_orb0
   start%vec(6) = start%vec(6) + dorb6
   start%vec(i) = start%vec(i) + del_orb(i)
-  call adjust_this
-  call track1 (start, ele, param, end2)
+  call track_this (start, ele, param, end2, q1, .true.)
   if (end2%state /= alive$) then
     call out_io (s_error$, r_name, 'PARTICLE LOST IN TRACKING (+). MATRIX NOT CALCULATED FOR ELEMENT: ' // ele%name)
     bmad_private%random_on = .true.
@@ -75,46 +80,78 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
   start = start_orb0
   start%vec(6) = start%vec(6) + dorb6
   start%vec(i) = start%vec(i) - del_orb(i)
-  call adjust_this
-  call track1 (start, ele, param, end1)
+  call track_this (start, ele, param, end1, q2, .true.)
   if (end1%state /= alive$) then
     call out_io (s_error$, r_name, 'PARTICLE LOST IN TRACKING (-). MATRIX NOT CALCULATED FOR ELEMENT: ' // ele%name)
     bmad_private%random_on = .true.
     return
   endif
 
   mat6(1:6, i) = (end2%vec - end1%vec) / (2 * del_orb(i))
+  if (do_spin) q_map(:, i) = (q2 - q1) / (2 * del_orb(i))
 enddo
 
 ! vestart_orb calc
 
-ele%vec0 = end_orb%vec - matmul(mat6, start_orb%vec)
-ele%mat6 = mat6
+if (.not. logic_option(.false., spin_only)) then
+  ele%vec0 = end_orb%vec - matmul(mat6, start_orb%vec)
+  ele%mat6 = mat6
+endif
+
+if (do_spin) then
+  ele%spin_q = q_map
+  call linear_to_spin_taylor(q_map, ele%spin_taylor)
+endif
 
 bmad_private%random_on = .true.
 err_flag = .false.
 
 !------------------------------------------------------
 contains
 
-subroutine adjust_this
+subroutine track_this(starting, ele, param, ending, q, adjust)
 
-if (start_orb%species == photon$) then
- call init_coord(start, start, ele, start_end$, start_orb%species)
- return
-endif
+type (coord_struct) starting, ending
+type (ele_struct) ele
+type (lat_param_struct) param
+
+real(rp) q(0:3), m(3,3)
+integer i
+logical adjust
 
 !
 
-call convert_pc_to (start%p0c * (1 + start%vec(6)), start%species, beta = start%beta)
+if (adjust) then
+  if (start_orb%species == photon$) then
+   call init_coord(starting, starting, ele, start_end$, start_orb%species)
+   return
+  endif
+
+  !
+
+  call convert_pc_to (starting%p0c * (1 + starting%vec(6)), starting%species, beta = starting%beta)
+
+  if (start_orb%beta == 0) then
+    starting%t = start_orb%t - starting%vec(5) / (c_light * starting%beta)
+  else
+    starting%t = start_orb%t - starting%vec(5) / (c_light * starting%beta) + start_orb%vec(5) / (c_light * start_orb%beta)
+  endif
+endif
+
+if (do_spin) then
+  do i = 1, 3
+    starting%spin = 0
+    starting%spin(i) = 1
+    call track1(starting, ele, param, ending)
+    m(:,i) = ending%spin
+  enddo
+  q = w_mat_to_quat(m)
 
-if (start_orb%beta == 0) then
-  start%t = start_orb%t - start%vec(5) / (c_light * start%beta)
 else
-  start%t = start_orb%t - start%vec(5) / (c_light * start%beta) + start_orb%vec(5) / (c_light * start_orb%beta)
+  call track1(starting, ele, param, ending)
 endif
 
-end subroutine adjust_this
+end subroutine track_this
 
 end subroutine