Homec4science
Diffusion COSOlver (master) Browse

cosolver/src/coulomb_speed_functions.h
05a553ae087cmaster

Last Change
Raw File

coulomb_speed_functions.h
View Options

!======================================================================!
!======================================================================!
!
! Test-Particle Speed functions
!
!======================================================================!
!======================================================================!
!-----------------------------------------------------------
!> @brief
!> Computes \f$ \alpha_{ab1}^{pl}\f$
!
!> \f[ \alpha_{ab1}^{pl} =4 \frac{\sigma}{\tau}l(\tau -2)-p (p -2 \sigma +4\frac{ \sigma}{\tau} +1) \f]
!
!> @param[in] p
!> @param[in] l
!>
!> @param[out] XOUT
!-----------------------------------------------------------
FUNCTION aplab1(p,l) RESULT(XOUT)
!
! CHECK
IMPLICIT NONE
!
INTEGER, INTENT(in) :: p,l
TYPE(FM) :: XOUT
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM('0.0')
!
XOUT = 4*l*sigmafm_/taufm_*( taufm_ -2) - p*(p - 2*sigmafm_ + 4*sigmafm_/taufm_ + 1 )
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
end function aplab1
!======================================================================!
!-----------------------------------------------------------
!> @brief
!> Computes \f$ \alpha_{ab2}^{pl}\f$
!
!> \f[ \alpha_{ab2}^{pl} =\frac{\sigma}{\tau} \left(4 l^2+4 l (p-1)+\frac{3}{2} p(p-1) \right)\f]
!
!> @param[in] p
!> @param[in] l
!>
!> @param[out] XOUT
!-----------------------------------------------------------
function aplab2(p,l) result(XOUT)
!
! CHECK
implicit none
!
integer, intent(in) :: p,l
TYPE(FM):: XOUT
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM('0.0')
!
XOUT = sigmafm_/taufm_*(4*(TO_FM(l)**2) + 4*l*(p -1) + TO_FM(3*p*(p -1))/2)
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION aplab2
!======================================================================!
!-----------------------------------------------------------
!> @brief
!> Computes \f$ \alpha_{ab1}^{'pl}\f$
!
!> \f[ a^{'pl}_{ab1} &= -2 (p+2 l) (\frac{\sigma}{\tau} (\tau -2)-1),
!
!> @param[in] p
!> @param[in] l
!>
!> @param[out] XOUT
!-----------------------------------------------------------
FUNCTION applab1(p,l) RESULT(XOUT)
!
! CHECK
IMPLICIT NONE
!
INTEGER, intent(in) :: p,l
TYPE(FM) :: XOUT
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM('0.0')
!
XOUT = - 2*( 2*l + p)*(sigmafm_/taufm_*(taufm_ - 2) -1 )
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION applab1
!======================================================================!
!-----------------------------------------------------------
!> @brief
!> Computes \f$ \alpha_{ab2}^{'pl}\f$
!
!> \f[ a^{'pl}_{ab2}= - \frac{\sigma}{\tau} \left(4l^2+4 l (p-1)+\frac{3}{2}p (p-1) \right) \f]
!
!> @param[in] p
!> @param[in] l
!>
!> @param[out] XOUT
!-----------------------------------------------------------
FUNCTION applab2(p,l) RESULT(XOUT)
!
! CHECK
IMPLICIT NONE
!
INTEGER, intent(in) :: p,l
TYPE(FM):: XOUT
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM('0.0')
!
XOUT = - sigmafm_/taufm_*(4*TO_FM(l)**2 + 4*l*(p - 1) + TO_FM(3)/2*p*(p -1))
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION applab2
!======================================================================!
!======================================================================!
FUNCTION nusTabpjgd(p,j,g,d) RESULT(XOUT)
!
! computes the speed function of the test particle part.
!
! TODO: Tested against Mathematica
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: p,j,g,d
!
INTEGER :: j1,n
TYPE(FM), SAVE :: Lpjl_,ekab0_,Erfkab0_,ekab1_,Erfkab1_,ekab2_,Erfkab2_,aplab0,applab0
TYPE(FM), SAVE :: AX1, AX2,APX1,APX2
TYPE(FM) :: XOUT
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = TO_FM('0.0')
!
!
! ! Precompute some quantities
! ekab0_ = ekab(d + j +1 + g)
! Erfkab0_ = Erfkab(d + j +g)
! ekab1_ = ekab(d + j + g )
! Erfkab1_ = Erfkab(d + j + g - 1)
! ekab2_ = ekab(d + j -1 + g )
! Erfkab2_ = Erfkab(d + j + g-2)
aplab0 = -4*sigmafm_*(taufm_ - TO_FM(1))/taufm_
applab0 = 4*(taufm_ -TO_FM(1))*(sigmafm_ + taufm_)/taufm_
DO j1=0,j
! Compute error function and its derivatives
!_______ INTRINSIC FUNCTION ____________
! ekab0_ = ekab(d + j1 +1 + g)
! Erfkab0_ = Erfkab(d + j1 +g)
! ekab1_ = ekab(d + j1 + g )
! Erfkab1_ = Erfkab(d + j1 + g - 1)
! ekab2_ = ekab(d + j1 -1 + g )
! Erfkab2_ = Erfkab(d + j1 + g-2)
!_______ PRECOMPUTED VALUES ____________
ekab0_ = ekab_array(d + j1 +1 + g)
Erfkab0_ = Erfkab_array(d + j1 +g)
ekab1_ = ekab_array(d + j1 + g )
Erfkab1_ = Erfkab_array(d + j1 + g - 1)
ekab2_ = ekab_array(d + j1 -1 + g )
Erfkab2_ = Erfkab_array(d + j1 + g-2)
Lpjl_ = Lpjl(real(p,dp),real(j,dp),real(j1,dp))
AX1 = aplab1(p,j1)
APX1 = applab1(p,j1)
AX2 = aplab2(p,j1)
APX2 = applab2(p,j1)
XOUT = XOUT + Lpjl_*(aplab0*Erfkab0_ + chifm_*applab0*ekab0_ + &
AX1*Erfkab1_ + chifm_*APX1*ekab1_ &
+ AX2*Erfkab2_ + chifm_*APX2*ekab2_ )
ENDDO
!
!
XOUT = XOUT*4*TO_FM(nu)/SQRT(PIFM)
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION nusTabpjgd
!______________________________________________________________________
FUNCTION sumnusTabpjgd(t,p,j) RESULT(XOUT)
! cumputes \sum_d^t L^p_{td} \nu^{Tpjpd}_{sab}
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: t,p,j
TYPE(FM) :: XOUT
!
! Local vars.
INTEGER :: d
TYPE(FM),SAVE :: Lptd
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = TO_FM(0)
DO d=0,t
Lptd = Lpjl(real(p,dp),real(t,dp),real(d,dp))
XOUT = XOUT + Lptd*nusTabpjgd(p,j,p,d)
ENDDO
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION sumnusTabpjgd
!
!======================================================================!
!======================================================================!
!
! Back-reaction Speed functions
!
!======================================================================!
!======================================================================!
FUNCTION bpplab1(p) RESULT(XOUT)
!
! CHECKED
IMPLICIT NONE
!
INTEGER,intent(in) :: p
TYPE(FM) :: XOUT
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = TO_FM('0.0')
!
!
XOUT = -8*(TO_FM(p)**2 + p -1)/TO_FM(2*p -1)/TO_FM(2*p +3)
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION bpplab1
!======================================================================!
FUNCTION bppplab0(p,l) RESULT(XOUT)
!
! CHECKED
IMPLICIT NONE
!
INTEGER,intent(in) :: p,l
TYPE(FM) :: XOUT
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM('0.0')
!
!
XOUT = - 8*TO_FM((p*(p-1)*(l+1)))/TO_FM( (2*p+1)*(2*p-1) ) - taufm_*8*(TO_FM(1 + p) - sigmafm_*p)/sigmafm_/TO_FM((2*p +1))
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION bppplab0
!======================================================================!
FUNCTION bppplab1(p) RESULT(XOUT)
!
! CHECKED
IMPLICIT NONE
!
INTEGER,intent(in) :: p
TYPE(FM) :: XOUT
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = 8*TO_FM((p+1)*(p+2))/TO_FM( (2*p +1)*(2*p +3))
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION bppplab1
!======================================================================!
FUNCTION bmplab0(p,l) RESULT(XOUT)
!
! CHECKED
IMPLICIT NONE
!
INTEGER,intent(in) :: p,l
TYPE(FM) :: XOUT
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = 4*TO_FM((p+1)*(p+2)*(2*p + 2*l +3))/TO_FM( (2*p +1)*(2*p +3)) &
+ 8*chifm_**2*(TO_FM( p) -p*sigmafm_ - sigmafm_)/TO_FM(2*p +1)
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION bmplab0
!======================================================================!
FUNCTION bmplab1(p) RESULT(XOUT)
!
! CHECKED
IMPLICIT NONE
!
INTEGER,intent(in) :: p
TYPE(FM) :: XOUT
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = - 8*TO_FM( p*(p-1))/TO_FM( (2*p +1)*(2*p-1))
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION bmplab1
!======================================================================!
FUNCTION betaeabgpdnj1(g,p,d,n,j1) RESULT(XOUT)
!
! CHECKED
IMPLICIT NONE
!
INTEGER, INTENT(in) :: g,p,d,n,j1
TYPE(FM) :: XOUT,X1
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM(0)
X1 = TO_FM(0)
!
! X1= ekab( g + 1 + d +j1 +n)
X1= ekab_array( g + 1 + d +j1 +n)
XOUT = chifm_**( p + 2*j1 + 2*n -1)*X1
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION betaeabgpdnj1
!======================================================================!
FUNCTION betapeabgpdnj1(g,p,d,n,j1) RESULT(XOUT)
!
! CHECKED
IMPLICIT NONE
!
INTEGER, INTENT(in) :: g,p,d,n,j1
TYPE(FM) :: XOUT
TYPE(FM),SAVE :: X1
INTEGER :: j2
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM('0.0')
X1 = TO_FM('0.0')
!
DO j2=0,j1
! X1= ekab(p + n +1 + d +j2)
X1= ekab_array(p + n +1 + d +j2)
XOUT = XOUT + chifm_**(2*j2)*FACTORIAL(TO_FM(j1))*X1/FACTORIAL(TO_FM(j2))/2
ENDDO
!
XOUT = chifm_**( p + 2*n -1)*XOUT
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION betapeabgpdnj1
!======================================================================!
FUNCTION betamErfabgpdnj1(p,d,n,j1) RESULT(XOUT)
!
! CHECKED
IMPLICIT NONE
!
INTEGER, INTENT(in) :: p,d,n,j1
TYPE(FM) :: XOUT
TYPE(FM), SAVE :: X1,X2,FACT
INTEGER :: j2
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = TO_FM('0.0')
!
! Precompute coefficients
!
FACT =Factorial(TO_FM(2*p+2*j1+1)/2)
X1= Erfkab_array(d+n)
! X1= Erfkab(d+n)
X1 = FACT*X1/2/FACTORIAL(TO_FM(1)/2)
!
DO j2=0,p+j1
! X2= ekab(d+n+j2+1)
X2= ekab_array(d+n+j2+1)
X2 = X2*FACT/2/FACTORIAL(TO_FM(2*j2 + 1)/2)*chifm_**(2*j2 +1)
XOUT = XOUT + X2
ENDDO
!
XOUT = (X1 - XOUT)*chifm_**(2*n - p-2)
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION betamErfabgpdnj1
!======================================================================!
FUNCTION nusFabpjgd(p,j,g,d) RESULT(XOUT)
! Compute the back-reaction speed function
!
! Check against Mathematica
IMPLICIT NONE
!
INTEGER, intent(in) :: p,j,g,d
TYPE(FM) :: XOUT
!
INTEGER :: j1
TYPE(FM), SAVE :: Lpjj1_,X1,X2
TYPE(FM), SAVE :: bpplab0_, bpplab1_,bppplab0_,bppplab1_,bmplab1_,bmplab0_
TYPE(FM), SAVE :: B10,B20,B30,B11,B21,B31
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = TO_FM('0.0')
! Precompute quantities (independent of j1)
bpplab0_ = 4*taufm_*taufm_/sigmafm_
bpplab1_ = bpplab1(p)
bppplab1_ = bppplab1(p)
bmplab1_ = bmplab1(p)
! Perform Associated Laguerre summation
DO j1=0,j
Lpjj1_= Lpjl(real(p,dp),real(j,dp),real(j1,dp))
bppplab0_ = bppplab0(p,j1)
bmplab0_ = bmplab0(p,j1)
! n=0
B10 = betaeabgpdnj1(g,p,d,0,j1)
B20 = betapeabgpdnj1(g,p,d,0,j1)
B30 = betamErfabgpdnj1(p,d,0,j1)
! n=1
B11 = betaeabgpdnj1(g,p,d,1,j1)
B21 = betapeabgpdnj1(g,p,d,1,j1)
B31 = betamErfabgpdnj1(p,d,1,j1)
XOUT = XOUT + Lpjj1_*( bpplab0_*B10 + bppplab0_*B20 + bmplab0_*B30 & ! n=0
+ bpplab1_*B11 + bppplab1_*B21 + bmplab1_*B31) ! n=1
ENDDO
!
XOUT = 4*XOUT*TO_FM(nu)/SQRT(PIFM)
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION nusFabpjgd
c4science · Help