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