1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
|
!
! wrap all possible functions together, they need unique ids
!
MODULE functions
USE mpfr
USE mpfr_ops
USE mpfr_cutoff_gamma
USE mpfr_yukawa
USE function_types
IMPLICIT NONE
INTEGER, SAVE :: functionid=-1
INTEGER, SAVE :: should_output=-1
CONTAINS
! evaluate the function, + nderiv other functions in the point X1,X2 in a reduced domain
! as a first step, this transforms X1,X2 into the natural variables of the original domain (e.g. R,T)
! write additionally to an open unit (should_output>0) connected to a file, for later checking the quality of the interpolation
SUBROUTINE f(res,nderiv,x1,x2)
INTEGER :: nderiv
TYPE(mpfr_type) :: x1,x2
TYPE(mpfr_type) :: res(0:nderiv)
TYPE(mpfr_type) :: T,r,upper,lower,zero
TYPE(mpfr_type) :: dummy(0:21)
INTEGER :: i
zero=.CONVERT."0"
SELECT CASE(functionid)
CASE(fun_trunc_coulomb_farfield,fun_trunc_coulomb_nearfield)
! 1 = farfield (R>11)
! 2 = nearfield (R<11)
SELECT CASE(functionid)
CASE (fun_trunc_coulomb_farfield)
! if R=+Infinity the result is zero
IF (mpfr_cmp(x2,zero)<=0) THEN
res=.CONVERT."0"
return
ENDIF
r=11/x2
upper=r**2 + 11*r + 50
lower=r**2 - 11*r
CASE (fun_trunc_coulomb_nearfield)
R=x2*11
upper=r**2 + 11*r + 50
lower=.CONVERT."0"
CASE DEFAULT
STOP "Function ID not implemented"
END SELECT
t=lower+x1*(upper-lower)
!t is zero, return the limiting expansion
IF (mpfr_cmp(t,zero)<=0) THEN
CALL cutoff_gamma_T0(21,R,dummy)
ELSE
CALL cutoff_gamma(21,t,r,dummy)
END IF
res(0:nderiv)=dummy(0:nderiv)
CASE(fun_yukawa)
! deal with infinite T,R locally
IF((X1==zero) .OR. (X2==zero)) THEN
res=zero
return
ELSE
T=(1/X1-1)**2
R=(1/X2-1)**2
CALL yukawa_gn_all(nderiv,T,R,dummy)
ENDIF
res(0:nderiv)=dummy(0:nderiv)
CASE DEFAULT
STOP "Function ID not implemented"
END SELECT
IF (should_output>0) THEN
!$OMP CRITICAL
WRITE(should_output,*) REAL(R),REAL(T),REAL(res(0:nderiv))
!$OMP END CRITICAL
ENDIF
END SUBROUTINE f
END MODULE functions
|
