Merge pull request #552 from bmad-sim/devel/step23

More bookkeeping correction.

bmad/spin/spin_concat_linear_maps.f90

@@ -1,5 +1,5 @@
 !+
-! Subroutine spin_concat_linear_maps (map1, branch, n1, n2, q_ele, orbit, excite_zero)
+! Subroutine spin_concat_linear_maps (err_flag, map1, branch, n1, n2, q_ele, orbit, excite_zero)
 !
 ! Routine to concatenate element spin/orbit maps in the range branch%ele(n1+1:n2)
 ! This routine will wrap around the ends of the lattice so n2 may be less than n1.

tao/code/tao_spin_polarization_calc.f90

@@ -29,10 +29,10 @@ subroutine tao_spin_polarization_calc (branch, tao_branch, excite_zero, ignore_k
 
 type (branch_struct), target :: branch
 type (tao_lattice_branch_struct), target :: tao_branch
+type (tao_spin_polarization_struct), pointer :: tspin
 type (coord_struct), pointer :: orbit(:)
-type (spin_orbit_map1_struct) :: q_1turn
 type (spin_orbit_map1_struct), pointer :: q1
-type (spin_orbit_map1_struct), target :: q_ele(0:branch%n_ele_track)
+
 type (ele_struct), pointer :: ele
 type (ele_pointer_struct), allocatable :: eles(:)
 
@@ -71,9 +71,10 @@ subroutine tao_spin_polarization_calc (branch, tao_branch, excite_zero, ignore_k
 !
 
 orbit => tao_branch%orbit
+tspin => tao_branch%spin
 if (.not. allocated(tao_branch%spin_ele)) allocate (tao_branch%spin_ele(0:branch%n_ele_track))
 tao_branch%spin_ele(:)%valid = .false.
-tao_branch%spin%valid = .false.
+tspin%valid = .false.
 if (branch%param%geometry == closed$ .and. orbit(branch%n_ele_track)%state /= alive$) return
 
 !
@@ -90,13 +91,14 @@ subroutine tao_spin_polarization_calc (branch, tao_branch, excite_zero, ignore_k
 endif
 
 if (.not. tao_branch%spin_map_valid) then
- call spin_concat_linear_maps (err, q_1turn, branch, 0, branch%n_ele_track, q_ele, orbit, excite_zero)
+ if (.not. allocated(tspin%q_ele)) allocate(tspin%q_ele(0:branch%n_ele_track))
+ call spin_concat_linear_maps (err, tspin%q_1turn, branch, 0, branch%n_ele_track, tspin%q_ele, orbit, excite_zero)
  if (err) return
  tao_branch%spin_map_valid = .true.
 endif
 
 if (branch%param%geometry == closed$) then 
- tao_branch%spin%tune = 2.0_rp * atan2(norm2(q_1turn%spin_q(1:3,0)), q_1turn%spin_q(0,0))
+ tspin%tune = 2.0_rp * atan2(norm2(tspin%q_1turn%spin_q(1:3,0)), tspin%q_1turn%spin_q(0,0))
  n0 = 0
  s_vec = [0.0_rp, 0.0_rp, 1.0_rp]
  dn_dpz = 0
@@ -128,21 +130,21 @@ subroutine tao_spin_polarization_calc (branch, tao_branch, excite_zero, ignore_k
  ele => branch%ele(ie)
 
  if (branch%param%geometry == closed$) then
- q_1turn = q_ele(ie) * q_1turn * map1_inverse(q_ele(ie))
- dn_dpz = spin_dn_dpz_from_qmap(q_1turn%orb_mat, q_1turn%spin_q, partial, partial2, err)
+ tspin%q_1turn = tspin%q_ele(ie) * tspin%q_1turn * map1_inverse(tspin%q_ele(ie))
+ dn_dpz = spin_dn_dpz_from_qmap(tspin%q_1turn%orb_mat, tspin%q_1turn%spin_q, partial, partial2, err)
  if (err) exit
  tao_branch%spin_ele(ie)%valid = .true.
  tao_branch%spin_ele(ie)%dn_dpz%vec = dn_dpz
  tao_branch%spin_ele(ie)%dn_dpz%partial = partial
  tao_branch%spin_ele(ie)%dn_dpz%partial2 = partial2
- n0 = q_1turn%spin_q(1:3, 0)
+ n0 = tspin%q_1turn%spin_q(1:3, 0)
  n0 = n0 / norm2(n0)
  if (ie == 0) cycle
 
  else
  if (ie == 0) cycle
- dn_dpz = quat_sym(ele, q_ele(ie)%spin_q, n0) + quat_rotate(q_ele(ie)%spin_q(:,0), tao_branch%spin_ele(ie-1)%dn_dpz%vec)
- n0 = quat_rotate(q_ele(ie)%spin_q(:,0), n0)
+ dn_dpz = quat_sym(ele, tspin%q_ele(ie)%spin_q, n0) + quat_rotate(tspin%q_ele(ie)%spin_q(:,0), tao_branch%spin_ele(ie-1)%dn_dpz%vec)
+ n0 = quat_rotate(tspin%q_ele(ie)%spin_q(:,0), n0)
  partial = 0
  partial2 = 0
  tao_branch%spin_ele(ie)%valid = .true.
@@ -194,31 +196,31 @@ subroutine tao_spin_polarization_calc (branch, tao_branch, excite_zero, ignore_k
 ! So test that integral_1ns is non-zero.
 
 if (integral_1ns == 0 .or. err) then
- tao_branch%spin%pol_limit_st = 0
- tao_branch%spin%pol_limit_dk = 0
- tao_branch%spin%pol_limit_dk_partial = 0
- tao_branch%spin%pol_limit_dk_partial2 = 0
- tao_branch%spin%pol_rate_bks = 0
- tao_branch%spin%depol_rate = 0
- tao_branch%spin%depol_rate_partial = 0
- tao_branch%spin%depol_rate_partial2 = 0
+ tspin%pol_limit_st = 0
+ tspin%pol_limit_dk = 0
+ tspin%pol_limit_dk_partial = 0
+ tspin%pol_limit_dk_partial2 = 0
+ tspin%pol_rate_bks = 0
+ tspin%depol_rate = 0
+ tspin%depol_rate_partial = 0
+ tspin%depol_rate_partial2 = 0
 else
  cm_ratio = charge_to_mass_of(branch%param%particle)
  call convert_pc_to ((1 + orbit(0)%vec(6)) * orbit(0)%p0c, branch%param%particle, gamma = gamma)
  f = f_rate * gamma**5 * cm_ratio**2 / branch%param%total_length
- tao_branch%spin%pol_limit_st = abs(f_limit * integral_bn / integral_1ns)
- tao_branch%spin%pol_limit_dk = abs(f_limit * (integral_bn - integral_bdn) / (integral_1ns + integral_dn2))
- tao_branch%spin%pol_limit_dk_partial = abs(f_limit * (integral_bn - integral_bdn_partial) / (integral_1ns + integral_dn2_partial))
- tao_branch%spin%pol_limit_dk_partial2 = abs(f_limit * (integral_bn - integral_bdn_partial2) / (integral_1ns + integral_dn2_partial2))
- tao_branch%spin%pol_rate_bks = f * integral_1ns
- tao_branch%spin%depol_rate = f * integral_dn2
- tao_branch%spin%depol_rate_partial = f * integral_dn2_partial
- tao_branch%spin%depol_rate_partial2 = f * integral_dn2_partial2
- tao_branch%spin%integral_bn = integral_bn
- tao_branch%spin%integral_bdn = integral_bdn
- tao_branch%spin%integral_1ns = integral_1ns
- tao_branch%spin%integral_dn2 = integral_dn2
- tao_branch%spin%valid = .true.
+ tspin%pol_limit_st = abs(f_limit * integral_bn / integral_1ns)
+ tspin%pol_limit_dk = abs(f_limit * (integral_bn - integral_bdn) / (integral_1ns + integral_dn2))
+ tspin%pol_limit_dk_partial = abs(f_limit * (integral_bn - integral_bdn_partial) / (integral_1ns + integral_dn2_partial))
+ tspin%pol_limit_dk_partial2 = abs(f_limit * (integral_bn - integral_bdn_partial2) / (integral_1ns + integral_dn2_partial2))
+ tspin%pol_rate_bks = f * integral_1ns
+ tspin%depol_rate = f * integral_dn2
+ tspin%depol_rate_partial = f * integral_dn2_partial
+ tspin%depol_rate_partial2 = f * integral_dn2_partial2
+ tspin%integral_bn = integral_bn
+ tspin%integral_bdn = integral_bdn
+ tspin%integral_1ns = integral_1ns
+ tspin%integral_dn2 = integral_dn2
+ tspin%valid = .true.
  if (present(err_flag)) err_flag = .false.
 endif
 

tao/code/tao_struct.f90

@@ -892,6 +892,8 @@ module tao_struct
  real(rp) :: integral_1ns = real_garbage$ ! Integral of g^3 (1 - 2(n * s_hat)/9)
  real(rp) :: integral_dn2 = real_garbage$ ! Integral of g^3 * 11 (dn/ddelta)^2 / 9
  logical :: valid = .false.
+ type (spin_orbit_map1_struct) :: q_1turn ! Save results from spin_concat_linear_maps in tao_spin_polarization.
+ type (spin_orbit_map1_struct), allocatable :: q_ele(:) ! Save results from spin_concat_linear_maps in tao_spin_polarization.
 end type
 
 ! For caching lattice calculations associated with plotting.