How call specific functions within a function?

Suppose that I have some functions

function V(r,Param1,Param2,Param3)result(result_V)
...
end function

function V2(r,Param1,Param2,Param3)result(result_V2)
...
end function

function V3(r,Param1,Param2,Param3)result(result_V3)
...
end function

Where the functions have different expressions. I want to invoke all of this one by one in a function,

Function delta(r,delta0,kcm,Param1,Param2,Param3)result(result_delta)!,Param1,Param2,Param3)
real*8, intent(in) :: r,delta0,kcm
real*8 :: result_delta
real*8 :: Param1,Param2,Param3

result_delta = (-V(r,Param1,Param2,Param3)*(sin(kcm*r + delta0))**2)/(kcm*frac) !,Param1,Param2,Param3)

end function

Is it possible to call the functions V, V2 and V3 one by one while I call delta function in main program? Like add some variable or method that modifies the V and change it for V2.

I first make one function with all different expressions, comment and uncomment for all cases, but I want to print results without modifying the code many times.

Welcome to the forum :slight_smile:

It is not clear to me what you want to do exactly, but I can think of various ways of achieving calls to these functions without changing the code. Probably the simplest is to use an extra variable in your function delta:

function delta( ..., choice ) result(...)
    ....
    select case (choice )
    case ( 1 ) 
         result_delta = V( ... )
    case ( 2 )
         result_delta = V2( ... )
    case ( ... )
         ...
    end select

The variable “choice” now comes in as an argument, but you can use a module variable to as well.

Does this help?

3 Likes

Another option is to use a procedure interface, assuming they all have the same interface:

abstract interface
    function potential_func(r, Param1, Param2, Param3) result(Vr)
        real(8), intent(in) :: r, Param1, Param2, Param3
        real(8) :: Vr
    end function potential_func
end interface

and then

function delta(r, delta0, kcm, Param1, Param2, Param3, V) result(result_delta)
    real(8), intent(in) :: r, delta0, kcm
    real(8), intent(in) :: Param1, Param2, Param3
    procedure(potential_func) :: V
    real(8) :: result_delta

    real(8), parameter :: frac = 1.0d0 ! Just assuming it's defined somewhere

    result_delta = (-V(r, Param1, Param2, Param3) * (sin(kcm*r + delta0))**2) / (kcm * frac)
end function delta
2 Likes

Thank you for the warm welcome.

Yes, you are right. What I try to do is obtain scattering phase shift (delta), putting some differential equation (result_delta) into a RK4 method for different potentials. With your suggest, I think to make one V function with all potentials as cases and choose them into delta. I do this:

!                                            POTENTIALS 
Function V(r,Param1,Param2,Param3,choice) result(result_V)           
Real*8, intent(in) :: r
Real*8 :: result_V
Real*8 :: Vr, Va, Mua
Real*8 :: Param1, Param2, Param3
Integer :: choice

    select case ( choice )
        case ( 1 )
            result_V = -Param2*exp(-Param3*r)/r + Param1*exp(-2*Param3*r)/r !-3.704*r)/r                       !MALFLIET-TJON
        case ( 2 )
            result_V = ( Param3*(Param3-1)*exp(-2*r/Param2)/( (1-exp(-r/Param2))**2 ) - Param1*exp(-r/Param2)/( 1-exp(-r/Param2) ) )/(Param2**2) !MANNING-Rosen
        case ( 3 )
            result_V = Param1*(exp(-2.0d0*(r-Param2)/Param3)-2.0d0*exp(-1.0d0*(r-Param2)/Param3))               !MORSE
        case ( 4 )
            result_V = ( -(Param1**2-Param2**2)*exp(-(Param1-Param2)*r)/(1-exp(-(Param1-Param2)*r)) ) + ( (Param1-Param2)**2*exp(-(Param1-Param2)*r)/((1-exp(Param1-Param2)*r)**2) )!Hulthen
    end select

end function

with delta function

Function deltaTrial(r,delta0,kcm,Param1,Param2,Param3,choice) result(dTrial)     !Trial delta
Real*8, intent(in) :: r, delta0,kcm
Real*8 :: dTrial,Param1,Param2,Param3
integer :: choice

dTrial = (-V(r,Param1,Param2,Param3,choice)*(sin(kcm*r + delta0))**2)/(kcm*frac)

end function

with RK4

!                                          RK4
function PFMRK4(kcom,Param1,Param2,Param3,choice)result(result_PFMRK4)
integer :: i,choice
real*8 :: r,y0,y1
real*8 :: a,b,h,N
real*8 :: result_PFMRK4
real*8 :: k1,k2,k3,k4
real*8, intent(in) :: kcom,Param1,Param2,Param3
N = 500000
a = 0.01d0
b = 5.0d0
h = (b - a)/N
            !! Initial conditions
r = a
y0 = 0.0d0

k1 = 0;k2 = 0;k3 = 0;k4 = 0
    do i=1, N
        k1 = deltaTrial( r,y0,kcom,Param1,Param2,Param3,choice)
        k2 = deltaTrial( r+0.5d0*h,y0+0.5d0*k1*h,kcom,Param1,Param2,Param3,choice)
        k3 = deltaTrial( r+0.5d0*h,y0+0.5d0*k2*h,kcom,Param1,Param2,Param3,choice)
        k4 = deltaTrial( r+h,y0+k3*h,kcom,Param1,Param2,Param3,choice)
        y1 = y0 + h*(k1 + 2.0d0*k2 + 2.0d0*k3 + k4)/6.0d0
    r = r + h
    y0 = y1
    end do
        result_PFMRK4 = y0*180.0d0/pi
end function

Thanks

I don’t understand very well. I should put all functions inside the abstract interface or put all expressions as cases inside the potential_func?

Unrelated to your question, but real*8 is not standard Fortran, although most compilers will accept it. real(8) is standard but not portable to all compilers, and for better approaches to declaring real variables see Floating Point Numbers — Fortran Programming Language.

2 Likes

Here is one possible refactor that removes the choice flag altogether and replaces it with a procedure argument that conforms to the shared abstract interface I showed in the previous post.
Everything lives in three small modules so you can mix-and-match any potential with the same RK4 integrator.

! =====================  potentials_mod.f90  =====================
module potentials_mod
   implicit none
   private       ! keep symbols local unless explicitly PUBLIC

   !---- Public API ----
   public :: potential_func                  ! the abstract interface name
   public :: V_MT, V_ManningRosen, V_Morse, V_Hulthén  ! concrete potentials

   !-----------------------------------------------------------------
   !  Generic “shape” that every concrete potential must match
   !-----------------------------------------------------------------
   abstract interface
      function potential_func(r, p1, p2, p3) result(Vr)
         real(8), intent(in) :: r, p1, p2, p3
         real(8)            :: Vr
      end function potential_func
   end interface

contains
!-------------------------------------------------------------------
!  1) Malfliet–Tjon (MT)
!-------------------------------------------------------------------
   function V_MT(r, p1, p2, p3) result(Vr)
      real(8), intent(in) :: r, p1, p2, p3
      real(8)             :: Vr
      Vr = -p2*exp(-p3*r)/r + p1*exp(-2*p3*r)/r
   end function V_MT

!-------------------------------------------------------------------
!  2) Manning–Rosen
!-------------------------------------------------------------------
   function V_ManningRosen(r, p1, p2, p3) result(Vr)
      real(8), intent(in) :: r, p1, p2, p3
      real(8)             :: Vr
      Vr = ( p3*(p3-1d0)*exp(-2d0*r/p2)        &
           / (1d0-exp(-r/p2))**2               &
           - p1*exp(-r/p2)/(1d0-exp(-r/p2)) )  &
           / (p2**2)
   end function V_ManningRosen

!-------------------------------------------------------------------
!  3) Morse
!-------------------------------------------------------------------
   function V_Morse(r, p1, p2, p3) result(Vr)
      real(8), intent(in) :: r, p1, p2, p3
      real(8)             :: Vr
      Vr = p1*(exp(-2d0*(r-p2)/p3) - 2d0*exp(-(r-p2)/p3))
   end function V_Morse

!-------------------------------------------------------------------
!  4) Hulthén
!-------------------------------------------------------------------
   function V_Hulthén(r, p1, p2, p3) result(Vr)
      real(8), intent(in) :: r, p1, p2, p3
      real(8)             :: Vr
      real(8) :: a
      a  = p1 - p2
      Vr = -(p1**2 - p2**2)*exp(-a*r)/(1d0 - exp(-a*r))           &
           + a**2 * exp(-a*r) / (1d0 - exp(-a*r))**2
   end function V_Hulthén

end module potentials_mod

! =====================  delta_mod.f90  ===========================
module delta_mod
   use potentials_mod          ! brings in the abstract interface
   implicit none
   private
   public :: deltaTrial

   real(8), parameter :: frac = 1.0d0   !<< you had this in the original code

contains
!-------------------------------------------------------------------
!  Trial phase-shift integrand
!-------------------------------------------------------------------
   function deltaTrial(r, delta0, kcm, p1, p2, p3, V) result(dTrial)
      real(8), intent(in)           :: r, delta0, kcm, p1, p2, p3
      procedure(potential_func)     :: V
      real(8)                       :: dTrial
      dTrial = -V(r, p1, p2, p3) * ( sin(kcm*r + delta0) )**2  &
               / (kcm*frac)
   end function deltaTrial
end module delta_mod
! =====================  rk4_mod.f90  =============================
module rk4_mod
   use delta_mod
   implicit none
   private
   public :: PFMRK4

contains
!-------------------------------------------------------------------
!  Simple RK4 driver that works with ANY potential obeying
!  the potential_func interface
!-------------------------------------------------------------------
   function PFMRK4(kcom, p1, p2, p3, V) result(delta_deg)
      real(8), intent(in)       :: kcom, p1, p2, p3
      procedure(potential_func) :: V         !<< key change
      real(8)                   :: delta_deg

      integer, parameter :: N = 500000
      real(8), parameter :: a = 0.01d0, b = 5.0d0
      real(8), parameter :: h = (b-a)/N
      real(8)            :: r, y, k1, k2, k3, k4
      integer            :: i

      r = a;  y = 0.0d0
      do i = 1, N
         k1 = deltaTrial(r,           y,              kcom, p1, p2, p3, V)
         k2 = deltaTrial(r+0.5d0*h,   y+0.5d0*k1*h,   kcom, p1, p2, p3, V)
         k3 = deltaTrial(r+0.5d0*h,   y+0.5d0*k2*h,   kcom, p1, p2, p3, V)
         k4 = deltaTrial(r+h,         y+k3*h,         kcom, p1, p2, p3, V)

         y  = y + h*(k1 + 2d0*k2 + 2d0*k3 + k4)/6d0
         r  = r + h
      end do

      delta_deg = y * 180d0 / acos(-1d0)     ! rad -> deg
   end function PFMRK4
end module rk4_mod
! =====================  main.f90  ================================
program demo
   use rk4_mod
   use potentials_mod   ! to get the specific potentials
   implicit none

   real(8) :: kcm, p1, p2, p3, delta

   kcm = 1.0d0
   p1  = 1.0d0
   p2  = 1.0d0
   p3  = 1.0d0

   ! ---- Choose any concrete potential simply by passing its name ----
   delta = PFMRK4(kcm, p1, p2, p3, V_MT)      ! Malfliet–Tjon
   write(*,'("Δ(MT) = ",f8.4,"°")') delta

   delta = PFMRK4(kcm, p1, p2, p3, V_Morse)   ! Morse
   write(*,'("Δ(Morse) = ",f8.4,"°")') delta
end program demo

5 Likes

Omg, you really open my mind. I have many ideas to apply, thanks!

You may find it useful to look up the documentations and examples for some of the routines in software libraries such as IMSL and NAG. For example, to find the root of a trancendental function, one of the arguments to pass to the library routine would be a procedure argument. The abstract interface may be used by the compiler to check that the actual argument has the correct interface.

I also found @FedericoPerini’s program instructive but it worked only after two changes: from V_Hulthén to V_Hulthen because the Fortran character set does not include é, and I needed

   public :: deltaTrial, potential_func

in the module delta_mod to make the abstract interface available in
PFMRK4.

1 Like