From 90795681f7c8bf1d0d606be6d29eed954b10001b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 12 Jan 2023 14:39:06 -0500 Subject: [PATCH] Fixed issue with rotation velocity calculation causing floating point errors. Made other improvements, including setting collisions to merge if they gain energy. --- src/base/base_module.f90 | 2 +- src/collision/collision_module.f90 | 38 ++++----- src/collision/collision_resolve.f90 | 4 +- src/collision/collision_util.f90 | 15 +++- src/fraggle/fraggle_generate.f90 | 117 +++++++++++++++++----------- 5 files changed, 105 insertions(+), 71 deletions(-) diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index 85fe74169..0c7bd8cf1 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -346,7 +346,7 @@ subroutine base_util_resize_storage(self, nnew) integer(I4B), intent(in) :: nnew !! New size ! Internals class(base_storage_frame), dimension(:), allocatable :: tmp - integer(I4B) :: i, n, nold, nbig + integer(I4B) :: i, nold, nbig nold = self%nframes if (nnew <= nold) return diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 361a27bd7..462043001 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -27,7 +27,12 @@ module collision integer(I4B), parameter :: COLLRESOLVE_REGIME_SUPERCATASTROPHIC = 3 integer(I4B), parameter :: COLLRESOLVE_REGIME_GRAZE_AND_MERGE = 4 integer(I4B), parameter :: COLLRESOLVE_REGIME_HIT_AND_RUN = 5 - character(len=*),dimension(5), parameter :: REGIME_NAMES = ["Merge", "Disruption", "Supercatastrophic", "Graze and Merge", "Hit and Run"] + character(len=NAMELEN),parameter :: REGIME_NAME_MERGE = "Merge" + character(len=NAMELEN),parameter :: REGIME_NAME_DISRUPTION = "Disruption" + character(len=NAMELEN),parameter :: REGIME_NAME_SUPERCATASTROPHIC = "Supercatastrophic" + character(len=NAMELEN),parameter :: REGIME_NAME_GRAZE_AND_MERGE = "Graze and Merge" + character(len=NAMELEN),parameter :: REGIME_NAME_HIT_AND_RUN = "Hit and Run" + character(len=NAMELEN),dimension(5), parameter :: REGIME_NAMES = [REGIME_NAME_MERGE, REGIME_NAME_DISRUPTION, REGIME_NAME_SUPERCATASTROPHIC, REGIME_NAME_GRAZE_AND_MERGE, REGIME_NAME_HIT_AND_RUN] !> Swiftest class for tracking pl-pl close encounters in a step when collisions are possible type, extends(encounter_list) :: collision_list_plpl @@ -52,14 +57,16 @@ module collision !> Class definition for the variables that describe the bodies involved in the collision type, extends(base_object) :: collision_impactors integer(I4B) :: ncoll !! Number of bodies involved in the collision - integer(I4B), dimension(:), allocatable :: id !! Index of bodies involved in the collision + integer(I4B), dimension(:), allocatable :: id !! Index of bodies involved in the collision real(DP), dimension(NDIM,2) :: rb !! Two-body equivalent position vectors of the collider bodies prior to collision in system barycentric coordinates real(DP), dimension(NDIM,2) :: vb !! Two-body equivalent velocity vectors of the collider bodies prior to collision in system barycentric coordinate real(DP), dimension(NDIM,2) :: rc !! Two-body equivalent position vectors of the collider bodies prior to collision in collision center of mass coordinates real(DP), dimension(NDIM,2) :: vc !! Two-body equivalent velocity vectors of the collider bodies prior to collision in collision center of mass coordinates real(DP), dimension(NDIM,2) :: rot !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(2) :: ke_orbit !! Orbital kinetic energy of each individual impactor + real(DP), dimension(2) :: ke_spin !! Spin kinetic energy of each individual impactors real(DP), dimension(NDIM,2) :: Ip !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision real(DP), dimension(2) :: Gmass !! Two-body equivalent G*mass of the collider bodies prior to the collision real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision @@ -113,10 +120,10 @@ module collision real(DP), dimension(NDIM,nbody) :: t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame real(DP), dimension(NDIM,nbody) :: n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame integer(I1B), dimension(nbody) :: origin_body !! Array of indices indicating which impactor body (1 or 2) the fragment originates from - real(DP), dimension(NDIM) :: L_orbit_tot !! Orbital angular momentum vector of all fragments - real(DP), dimension(NDIM) :: L_spin_tot !! Spin angular momentum vector of all fragments - real(DP), dimension(NDIM,nbody) :: L_orbit !! Orbital angular momentum vector of each individual fragment - real(DP), dimension(NDIM,nbody) :: L_spin !! Spin angular momentum vector of each individual fragment + real(DP), dimension(NDIM) :: L_orbit_tot !! Orbital angular momentum vector of all fragments + real(DP), dimension(NDIM) :: L_spin_tot !! Spin angular momentum vector of all fragments + real(DP), dimension(NDIM,nbody) :: L_orbit !! Orbital angular momentum vector of each individual fragment + real(DP), dimension(NDIM,nbody) :: L_spin !! Spin angular momentum vector of each individual fragment real(DP) :: ke_orbit_tot !! Orbital kinetic energy of all fragments real(DP) :: ke_spin_tot !! Spin kinetic energy of all fragments real(DP) :: pe !! Potential energy of all fragments @@ -127,7 +134,6 @@ module collision procedure :: dealloc => collision_util_dealloc_fragments !! Deallocates all allocatable arrays and sets everything else to 0 procedure :: reset => collision_util_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) procedure :: set_coordinate_system => collision_util_set_coordinate_fragments !! Sets the coordinate system of the fragments - final :: collision_final_fragments !! Finalizer deallocates all allocatables end type collision_fragments @@ -189,7 +195,7 @@ module collision integer(I4B) :: stage_dimid !! ID for the stage dimension integer(I4B) :: stage_varid !! ID for the stage variable character(NAMELEN) :: stage_dimname = "stage" !! name of the stage dimension (before/after) - character(len=6), dimension(2) :: stage_coords = ["before", "after"] !! The stage coordinate labels + character(len=6), dimension(2) :: stage_coords = ["before", "after "] !! The stage coordinate labels integer(I4B) :: collision_id_dimid !! ID for the collision event dimension @@ -500,18 +506,6 @@ end subroutine collision_util_set_original_scale_factors contains - subroutine collision_final_fragments(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(collision_fragments(*)), intent(inout) :: self - - if (allocated(self%info)) deallocate(self%info) - return - end subroutine collision_final_fragments - subroutine collision_final_impactors(self) !! author: David A. Minton !! diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 52948c4d2..4b39eb1b7 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -509,7 +509,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) real(DP) :: E_before, E_after, dLmag real(DP), dimension(NDIM) :: L_orbit_before, L_orbit_after, L_spin_before, L_spin_after, L_before, L_after, dL_orbit, dL_spin, dL logical :: lplpl_collision - character(len=STRMAX) :: timestr + character(len=STRMAX) :: timestr, idstr integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision integer(I4B) :: i, loop, ncollisions @@ -562,6 +562,8 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) param%maxid_collision = max(param%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id)) param%maxid_collision = param%maxid_collision + 1 collider%collision_id = param%maxid_collision + write(idstr,*) collider%collision_id + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Collision id " // trim(adjustl(idstr))) ! Get the collision regime call impactors%get_regime(nbody_system, param) diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 797950d88..d33df2cfb 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -233,7 +233,14 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par self%pe(phase_val) = constraint_system%pe / self%Escale self%be(phase_val) = constraint_system%be / self%Escale self%te(phase_val) = constraint_system%te / self%Escale - if (phase_val == 2) then + if (phase_val == 1) then + do concurrent(i = 1:2) + impactors%ke_orbit(i) = 0.5_DP * impactors%mass(i) * dot_product(impactors%vc(:,i), impactors%vc(:,i)) + impactors%ke_spin(i) = 0.5_DP * impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * dot_product(impactors%rot(:,i), impactors%rot(:,i)) + impactors%L_orbit(:,i) = impactors%mass(i) * impactors%rc(:,i) .cross. impactors%vc(:,i) + impactors%L_spin(:,i) = impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * impactors%rot(:,i) + end do + else if (phase_val == 2) then do concurrent(i = 1:nfrag) fragments%ke_orbit(i) = 0.5_DP * fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) fragments%ke_spin(i) = 0.5_DP * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i), fragments%rot(:,i)) @@ -315,6 +322,8 @@ module subroutine collision_util_dealloc_impactors(self) self%rot(:,:) = 0.0_DP self%L_spin(:,:) = 0.0_DP self%L_orbit(:,:) = 0.0_DP + self%ke_spin(:) = 0.0_DP + self%ke_orbit(:) = 0.0_DP self%Ip(:,:) = 0.0_DP self%mass(:) = 0.0_DP self%radius(:) = 0.0_DP @@ -810,6 +819,8 @@ module subroutine collision_util_set_natural_scale_factors(self) impactors%radius(:) = impactors%radius(:) / collider%dscale impactors%L_spin(:,:) = impactors%L_spin(:,:) / collider%Lscale impactors%L_orbit(:,:) = impactors%L_orbit(:,:) / collider%Lscale + impactors%ke_orbit(:) = impactors%ke_orbit(:) / collider%Escale + impactors%ke_spin(:) = impactors%ke_spin(:) / collider%Escale do concurrent(i = 1:2) impactors%rot(:,i) = impactors%L_spin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) @@ -859,6 +870,8 @@ module subroutine collision_util_set_original_scale_factors(self) impactors%vc = impactors%vc * collider%vscale impactors%L_spin = impactors%L_spin * collider%Lscale impactors%L_orbit = impactors%L_orbit * collider%Lscale + impactors%ke_orbit = impactors%ke_orbit * collider%Escale + impactors%ke_spin = impactors%ke_spin * collider%Escale do concurrent(i = 1:2) impactors%rot(:,i) = impactors%L_spin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) end do diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 611a34eb3..438a33462 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -27,6 +27,7 @@ module subroutine fraggle_generate(self, nbody_system, param, t) ! Internals integer(I4B) :: i, ibiggest, nfrag character(len=STRMAX) :: message + logical :: lfailure select type(nbody_system) class is (swiftest_nbody_system) @@ -49,7 +50,12 @@ module subroutine fraggle_generate(self, nbody_system, param, t) call base_util_exit(FAILURE) end select call self%set_mass_dist(param) - call self%disrupt(nbody_system, param, t) + call self%disrupt(nbody_system, param, t, lfailure) + if (lfailure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Treating this as a merger.") + call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge + return + end if associate (fragments => self%fragments) ! Populate the list of new bodies @@ -128,10 +134,12 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call fraggle_generate_pos_vec(self) call fraggle_generate_rot_vec(self, nbody_system, param) call fraggle_generate_vel_vec(self, nbody_system, param, lfailure_local) + call self%get_energy_and_momentum(nbody_system, param, phase="after") dL = self%L_total(:,2)- self%L_total(:,1) dE = self%te(2) - self%te(1) + lfailure_local = (dE > 0.0_DP) call swiftest_io_log_one_message(COLLISION_LOG_OUT, "All quantities in collision system natural units") call swiftest_io_log_one_message(COLLISION_LOG_OUT, "* Conversion factors (collision system units / nbody system units):") @@ -355,7 +363,8 @@ end subroutine fraggle_generate_pos_vec module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! - !! Calculates the spins of a collection of fragments such that they conserve angular momentum + !! Computes an initial "guess" for the rotation states of fragments based on angular momentum and energy constraints. + !! These will be adjusted later when the final fragment velocities are computed in fraggle_generate_vel_vec implicit none ! Arguments class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object @@ -363,34 +372,48 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals real(DP), dimension(NDIM) :: Lbefore, Lafter, L_spin, rotdir - real(DP) :: v_init, v_final, mass_init, mass_final, rotmag + real(DP) :: v_init, v_final, mass_init, mass_final, rotmag, dKE, KE_init, KE_final real(DP), parameter :: random_scale_factor = 0.01_DP !! The relative scale factor to apply to the random component of the rotation vector integer(I4B) :: i associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) - ! Torque the first body based on the change in angular momentum betwen the pre- and post-impact system assuming a simple bounce + ! We will start by assuming that kinetic energy gets partitioned such that the change in kinetic energy of body 1 is equal to the + ! change in kinetic energy between bodies 2 and all fragments. This will then be used to compute a torque on body/fragment 1. + ! All other fragments will be given a random velocity with a magnitude scaled by the change in the orbital system angular momentum mass_init = impactors%mass(2) mass_final = sum(fragments%mass(2:nfrag)) v_init = .mag.(impactors%vb(:,2) - impactors%vb(:,1)) - v_final = sqrt(mass_init / mass_final * v_init**2 - 2 * impactors%Qloss / mass_final) + KE_init = 0.5_DP * mass_init * v_init**2 + + ! Initialize fragment rotations and velocities to be pre-impact rotations in order to compute the energy. This will get adjusted later + fragments%rot(:,1) = impactors%rot(:,1) + fragments%vc(:,1) = impactors%vc(:,1) + do concurrent(i = 2:nfrag) + fragments%rot(:,i) = impactors%rot(:,2) + fragments%vc(:,i) = impactors%vc(:,2) + end do + call collider%get_energy_and_momentum(nbody_system, param, phase="after") + dKE = 0.5_DP * (collider%pe(2) - collider%pe(1) + collider%be(2) - collider%be(1) - impactors%Qloss) + KE_final = max(KE_init + dKE,0.0_DP) + + v_final = sqrt(2 * KE_final / mass_final) Lbefore(:) = mass_init * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1)) Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:)) L_spin(:) = impactors%L_spin(:,1) + (Lbefore(:) - Lafter(:)) fragments%rot(:,1) = L_spin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) - + fragments%rotmag(1) = .mag.fragments%rot(:,1) ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random + call random_number(fragments%rotmag(2:nfrag)) do concurrent(i = 2:nfrag) call random_number(rotdir) rotdir = rotdir - 0.5_DP rotdir = .unit. rotdir - call random_number(rotmag) - rotmag = random_scale_factor * rotmag * .mag. (Lbefore(:) - Lafter(:)) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) - fragments%rot(:,i) = rotmag * rotdir + rotmag = random_scale_factor * .mag. (Lbefore(:) - Lafter(:)) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) + fragments%rot(:,i) = rotmag * fragments%rotmag(i) * rotdir end do - end associate return @@ -416,10 +439,10 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP) :: vmag, vesc, dE, E_residual, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove - real(DP), parameter :: vmin_initial_factor = 1.5_DP ! For the initial "guess" of fragment velocities, this is the maximum velocity relative to escape velocity that the fragments will have - real(DP), parameter :: vmax_initial_factor = 5.0_DP ! For the initial "guess" of fragment velocities, this is the maximum velocity relative to escape velocity that the fragments will have - integer(I4B), parameter :: MAXLOOP = 200 - integer(I4B), parameter :: MAXTRY = 20 + real(DP), parameter :: vmin_initial_factor = 2.0_DP ! For the initial "guess" of fragment velocities, this is the maximum velocity relative to escape velocity that the fragments will have + real(DP) :: vmax_initial_factor = 8.0_DP + integer(I4B), parameter :: MAXLOOP = 20 + integer(I4B), parameter :: MAXTRY = 100 real(DP), parameter :: SUCCESS_METRIC = 1.0e-3_DP class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -440,6 +463,13 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu vsign(:) = 1 end where + ! Hit and run collisions should only affect the runner + if (lhitandrun) then + istart = 2 + else + istart = 1 + end if + ! The minimum fragment velocity will be set by the escape velocity if (lhitandrun) then vesc = sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) @@ -448,7 +478,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end if ! Scale the magnitude of the velocity by the distance from the impact point - ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct + ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct do concurrent(i = 1:nfrag) rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))) @@ -459,40 +489,35 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu fscale = log(vmax_initial_factor - vmin_initial_factor + 1.0_DP) / log(maxval(vscale(:))) vscale(:) = vscale(:)**fscale + vmin_initial_factor - 1.0_DP - ! Set the velocities of all fragments using all of the scale factors determined above - do concurrent(i = 1:nfrag) - j = fragments%origin_body(i) - vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) - if (lhitandrun) then - if (i == 1) then - fragments%vc(:,1) = impactors%vc(:,1) - else - vmag = .mag.impactors%vc(:,2) / maxval(vscale(:)) - fragments%vc(:,i) = vmag * vscale(i) * impactors%bounce_unit(:) * vsign(i) + vrot(:) - end if - else - ! Add more velocity dispersion to disruptions vs hit and runs. - vmag = vesc * vscale(i) - rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) - vimp_unit(:) = .unit. rimp(:) - fragments%vc(:,i) = vmag * vscale(i) * vimp_unit(:) + vrot(:) - end if - end do - - ! Hit and run collisions should only affect the runner - if (lhitandrun) then - istart = 2 - else - istart = 1 - end if - - ! Every time the collision-frame velocities are altered, we need to be sure to shift everything back to the center-of-mass frame - call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) - call fragments%set_coordinate_system() E_residual_best = huge(1.0_DP) lfailure = .false. dE_metric = huge(1.0_DP) + outer: do try = 1, MAXTRY + ! Set the velocities of all fragments using all of the scale factors determined above + do concurrent(i = 1:nfrag) + j = fragments%origin_body(i) + vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) + if (lhitandrun) then + if (i == 1) then + fragments%vc(:,1) = impactors%vc(:,1) + else + vmag = .mag.impactors%vc(:,2) / maxval(vscale(:)) + fragments%vc(:,i) = vmag * vscale(i) * impactors%bounce_unit(:) * vsign(i) + vrot(:) + end if + else + ! Add more velocity dispersion to disruptions vs hit and runs. + vmag = vesc * vscale(i) + rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) + vimp_unit(:) = .unit. rimp(:) + fragments%vc(:,i) = vmag * vscale(i) * vimp_unit(:) + vrot(:) + end if + end do + + ! Every time the collision-frame velocities are altered, we need to be sure to shift everything back to the center-of-mass frame + call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) + call fragments%set_coordinate_system() + do loop = 1, MAXLOOP call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body @@ -527,7 +552,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ke_remove = min(f_orbit * E_residual, ke_avail) f_orbit = ke_remove / E_residual - fscale = sqrt((fragments%ke_orbit_tot - ke_remove)/fragments%ke_orbit_tot) + fscale = sqrt((max(fragments%ke_orbit_tot - ke_remove, 0.0_DP))/fragments%ke_orbit_tot) fragments%vc(:,:) = fscale * fragments%vc(:,:) f_spin = 1.0_DP - f_orbit