diff --git a/bmad/spin/spin_concat_linear_maps.f90 b/bmad/spin/spin_concat_linear_maps.f90 index 79e93c70a7..be4abb430b 100644 --- a/bmad/spin/spin_concat_linear_maps.f90 +++ b/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. diff --git a/tao/code/tao_spin_polarization_calc.f90 b/tao/code/tao_spin_polarization_calc.f90 index 14bb465a97..b92dbfca71 100644 --- a/tao/code/tao_spin_polarization_calc.f90 +++ b/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 diff --git a/tao/code/tao_struct.f90 b/tao/code/tao_struct.f90 index 7cb496ef4e..6b54d77ea9 100644 --- a/tao/code/tao_struct.f90 +++ b/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.