Massive restructuring to help clarify the symba_frag_pos code. Split units of work into nested subroutines with clear names of their purpose. Simplified non-linear solver to only consider the radial component. Lots of other small improvements.

src/modules/module_interfaces.f90

@@ -726,7 +726,7 @@ function symba_casedisruption (symba_plA, family, nmergeadd, mergeadd_list, x, v
          real(DP), dimension(:,:), intent(inout)   :: x, v, L_spin, Ip
          real(DP), dimension(:),   intent(inout)   :: mass, radius, mass_res
          type(user_input_parameters),intent(inout) :: param
-         real(DP), intent(in)                      :: Qloss
+         real(DP), intent(inout)                   :: Qloss
          integer(I4B)                              :: status
       end function symba_casedisruption
 
@@ -744,7 +744,7 @@ function symba_casehitandrun (symba_plA, family, nmergeadd, mergeadd_list, name,
          real(DP), dimension(:,:), intent(inout)   :: x, v, Lspin, Ip
          real(DP), dimension(:), intent(inout)     :: mass, radius, mass_res
          type(user_input_parameters),intent(inout) :: param
-         real(DP), intent(in)                      :: Qloss
+         real(DP), intent(inout)                   :: Qloss
          integer(I4B)                              :: status
       end function symba_casehitandrun
 
@@ -776,7 +776,7 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis
          real(DP), dimension(:,:), intent(in)      :: x, v, lspin, Ip
          real(DP), dimension(:), intent(in)        :: mass, radius, mass_res
          type(user_input_parameters),intent(inout) :: param
-         real(DP), intent(in)                      :: Qloss
+         real(DP), intent(inout)                   :: Qloss
          integer(I4B)                              :: status
       end function symba_casesupercatastrophic
      end interface
@@ -929,7 +929,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
             type(user_input_parameters), intent(in)   :: param 
             type(symba_pl), intent(inout)             :: symba_plA
             integer(I4B), dimension(:), intent(in)    :: family
-            real(DP), intent(in)                      :: Qloss
+            real(DP), intent(inout)                   :: Qloss
             real(DP), dimension(:,:), intent(inout)   :: x, v, L_spin, Ip
             real(DP), dimension(:), intent(inout)     :: mass, radius, m_frag, rad_frag
             real(DP), dimension(:,:), intent(inout)   :: Ip_frag

src/modules/module_symba.f90

@@ -103,78 +103,4 @@ MODULE module_symba
 
    end type symba_merger 
 
-   type, public, extends(lambda_obj) :: symba_frag_lambda
-      procedure(abstract_objective_func), pointer, nopass :: ke_objective_func_ptr => null()
-      real(DP), dimension(:),   allocatable :: m_frag
-      real(DP), dimension(:,:), allocatable :: x_frag
-      real(DP), dimension(:,:), allocatable :: v_frag
-      real(DP), dimension(NDIM)             :: L_target
-      real(DP)                              :: ke_target
-   contains
-      generic   :: init => ke_objective_func_init
-      procedure :: eval => ke_objective_func_eval
-      procedure, nopass :: ke_objective_func_init
-      final     :: ke_objective_func_destroy
-   end type symba_frag_lambda
-   interface symba_frag_lambda
-      module procedure ke_objective_func_init
-   end interface
-
-   abstract interface
-      function abstract_objective_func(vflat, v_frag, x_frag, m_frag, L_target, ke_target) result(fnorm)
-         ! Template for the kinetic energy constraint function used for minimizing
-         import DP
-         real(DP), dimension(:),   intent(in) :: vflat          !! Unrolled unknown velocity vectors
-         real(DP), dimension(:,:), intent(in) :: v_frag, x_frag !! Velocity and position vectors
-         real(DP), dimension(:),   intent(in) :: m_frag         !! Fragment masses
-         real(DP), dimension(:),   intent(in) :: L_target       !! Target orbital momentum
-         real(DP),                 intent(in) :: ke_target      !! Target kinetic energ
-         real(DP)                             :: fnorm          !! The objective function result: norm of the vector composed of the tangential momentum and energy
-      end function
-   end interface
-
-   contains
-      type(symba_frag_lambda) function ke_objective_func_init(lambda, v_frag, x_frag, m_frag, L_target, ke_target)
-         implicit none
-         ! Arguments
-         procedure(abstract_objective_func)       :: lambda  !! The lambda function
-         real(DP), dimension(:,:), intent(in) :: v_frag, x_frag !! Velocity and position vectors
-         real(DP), dimension(:),   intent(in) :: m_frag         !! Fragment masses
-         real(DP), dimension(:),   intent(in) :: L_target       !! Target orbital momentum
-         real(DP),                 intent(in) :: ke_target      !! Target kinetic energ
-         ! Internals
-         associate(self => ke_objective_func_init)
-            self%ke_objective_func_ptr  => lambda
-            allocate(self%m_frag, source=m_frag)
-            allocate(self%v_frag, source=v_frag)
-            allocate(self%x_frag, source=x_frag)
-            self%L_target(:) = L_target(:)
-            self%ke_target = ke_target
-         end associate
-         return
-      end function ke_objective_func_init
-
-      subroutine ke_objective_func_destroy(self)
-         implicit none
-         type(symba_frag_lambda) :: self
-         if (allocated(self%m_frag)) deallocate(self%m_frag)
-         if (allocated(self%v_frag)) deallocate(self%v_frag)
-         if (allocated(self%x_frag)) deallocate(self%x_frag)
-         if (associated(self%ke_objective_func_ptr)) nullify(self%ke_objective_func_ptr)
-      end subroutine ke_objective_func_destroy 
-
-   function ke_objective_func_eval(self, x) result(fnorm)
-      implicit none
-      ! Arguments
-      class(symba_frag_lambda), intent(in) :: self
-      real(DP), dimension(:),   intent(in) :: x
-      ! Result
-      real(DP)                      :: fnorm
-
-      if (associated(self%ke_objective_func_ptr)) then
-         fnorm = self%ke_objective_func_ptr(x, self%v_frag, self%x_frag, self%m_frag, self%L_target, self%ke_target)
-      else
-         error stop "KE Objective function was not initialized."
-      end if
-   end function ke_objective_func_eval
 END MODULE module_symba