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