Started building an objective function for the KE constraint using th…

…e lambda function class.
MintonGroup · May 13, 2021 · b643779 · b643779
1 parent 32dbbcc
commit b643779
Show file tree

Hide file tree

Showing 4 changed files with 22 additions and 6 deletions.
diff --git a/Makefile b/Makefile
@@ -47,14 +47,14 @@
 SWIFTEST_MODULES =   swiftest_globals.f90 \
                      user.f90 \
                      swiftest_data_structures.f90 \
+							lambda_function.f90\
                      module_swifter.f90 \
                      module_helio.f90 \
                      module_nrutil.f90 \
                      module_symba.f90 \
                      module_swiftestalloc.f90 \
                      io.f90 \
                      module_interfaces.f90 \
-							lambda_function.f90\
                      swiftest.f90 
 
 include Makefile.Defines

diff --git a/src/modules/lambda_function.f90 b/src/modules/lambda_function.f90
@@ -13,11 +13,11 @@ module lambda_function
    end type
 
    abstract interface
-      function lambda0(x)
+      function lambda0(x) result(y)
          ! Template for a 0 argument function
          import DP
          real(DP), intent(in) :: x
-         real(DP) :: lambda
+         real(DP) :: y
       end function
    end interface
 

diff --git a/src/modules/module_symba.f90 b/src/modules/module_symba.f90
@@ -27,6 +27,7 @@ MODULE module_symba
 
    USE swiftest_globals
    USE module_helio
+   use lambda_function
    IMPLICIT NONE
 
    INTEGER(I4B), PARAMETER :: NENMAX = 32767
@@ -101,6 +102,21 @@ MODULE module_symba
       type(swiftest_particle_info), dimension(:), allocatable :: info
 
    end type symba_merger 
+
+   type, public, extends(lambda_obj) :: symba_vel_lambda_obj
+      procedure(abstract_objective_func), pointer, nopass :: ke_objective_func_ptr => null()
+   end type symba_vel_lambda_obj
+
+   abstract interface
+      function abstract_objective_func(v, m, rhat, t, G, B, L) result(fnorm)
+         ! Template for the kinetic energy constraint function used for minimizing
+         import DP
+         real(DP), dimension(:),   intent(in) :: v, m, t, G
+         real(DP), dimension(:,:), intent(in) :: rhat
+         real(DP),                 intent(in) :: B, L      
+         real(DP)                             :: fnorm
+      end function
+   end interface
 END MODULE module_symba
 !**********************************************************************************************************************************
 !

diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90
@@ -370,7 +370,7 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr
       end do
       if (Lambda > 0.0_DP) then
          lmerge = .false.
-         v_r_mag(:) = symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambda, tau)
+         v_r_mag(:) = symba_frag_pos_fragment_velocity(nfrag, m_frag, v_r_unit, Lambda, tau)
          do i = 1, nfrag
             v_frag(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_t(:, i)
          end do
@@ -385,12 +385,12 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr
       return
    end subroutine symba_frag_pos_kinetic_energy
 
-   function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambda, tau) result(v_r_mag)
+
+   function symba_frag_pos_fragment_velocity(nfrag, m_frag, v_r_unit, Lambda, tau) result(v_r_mag)
       implicit none
       ! Arguments
       integer(I4B),             intent(in)  :: nfrag
       real(DP), dimension(:),   intent(in)  :: m_frag                   !! Fragment masses
-      real(DP), dimension(:,:), intent(in)  :: x_frag                   !! Fragment position in the center of mass frame 
       real(DP), dimension(:,:), intent(in)  :: v_r_unit                 !! Radial velocity unit vector for each fragment
       real(DP), intent(in)                  :: Lambda                   !! Sum of the radial kinetic energy of all fragments 
       real(DP), dimension(:), intent(in)    :: tau                      !! Sum of the tangential momentum vector of all fragments