FUNCTION hermite(npt, func, ier) RESULT(fn_val) ! HERMITE INTEGRATION, I.E. EVALUATION OF THE INTEGRAL OF ! exp(-x**2)*f(x) FROM -INFINITY TO +INFINITY ! NPT = INPUT, NO. OF POINTS AT WHICH f(x) IS TO BE EVALUATED. ! NPT MUST BE ONE OF 2, 4, 6, 8, 10, 12, 16, 20. ! FUNC = INPUT, NAME OF THE USER'S FUNCTION TO SUPPLY VALUES ! OF F(X). THE FUNCTION MAY HAVE ANY NAME BUT ONLY THE ! ONE ARGUMENT, X. ! IER = OUTPUT, ERROR INDICATOR ! = 0 NO ERROR DETECTED ! = 1 ILLEGAL VALUE OF NPT ! LATEST REVISION - 5 December 1999 !************************************************************************ IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) INTEGER, INTENT(IN) :: npt INTEGER, INTENT(OUT) :: ier REAL (dp) :: fn_val INTERFACE FUNCTION func(x) RESULT(fn_val) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: x REAL (dp) :: fn_val END FUNCTION func END INTERFACE ! Local variables REAL (dp) :: xpt(39) = (/ 0.707106781186548_dp, & 0.524647623275290_dp, 1.650680123885785_dp, & 0.436077411927617_dp, 1.335849074013697_dp, 2.350604973674492_dp, & 0.381186990207322_dp, 1.157193712446780_dp, 1.981656756695843_dp, & 2.930637420257244_dp, & 0.342901327223705_dp, 1.036610829789514_dp, 1.756683649299882_dp, & 2.532731674232790_dp, 3.436159118837738_dp, & 0.314240376254359_dp, 0.947788391240164_dp, 1.597682635152605_dp, & 2.279507080501060_dp, 3.020637025120890_dp, 3.889724897869782_dp, & 0.27348104613815_dp, 0.82295144914466_dp, 1.38025853919888_dp, & 1.95178799091625_dp, 2.54620215784748_dp, 3.17699916197996_dp, & 3.86944790486012_dp, 4.68873893930582_dp, & 0.2453407083009_dp, 0.7374737285454_dp, 1.2340762153953_dp, & 1.7385377121166_dp, 2.2549740020893_dp, 2.7888060584281_dp, & 3.3478545673832_dp, 3.9447640401156_dp, 4.6036824495507_dp, & 5.3874808900112_dp /), & wt(39) = (/ 8.862269254528D-1, & 8.049140900055D-1, 8.131283544725D-2, & 7.246295952244D-1, 1.570673203229D-1, 4.530009905509D-3, & 6.611470125582D-1, 2.078023258149D-1, 1.707798300741D-2, & 1.996040722114D-4, & 6.108626337353D-1, 2.401386110823D-1, 3.387439445548D-2, & 1.343645746781D-3, 7.640432855233D-6, & 5.701352362625D-1, 2.604923102642D-1, 5.160798561588D-2, & 3.905390584629D-3, 8.573687043588D-5, 2.658551684356D-7, & 5.079294790166D-1, 2.806474585285D-1, 8.381004139899D-2, & 1.288031153551D-2, 9.322840086242D-4, 2.711860092538D-5, & 2.320980844865D-7, 2.654807474011D-10, & 4.622436696006D-1, 2.866755053628D-1, 1.090172060200D-1, & 2.481052088746D-2, 3.243773342238D-3, 2.283386360163D-4, & 7.802556478532D-6, 1.086069370769D-7, 4.399340992273D-10, & 2.229393645534D-13 /), zero = 0.0_dp, x INTEGER :: npts(8) = (/ 2, 4, 6, 8, 10, 12, 16, 20 /), ipos, i, nby2 ! CHECK FOR PERMISSIBLE VALUE OF NPT. fn_val = zero ipos = 1 DO i = 1, 8 IF(npts(i).EQ.npt) GO TO 20 ipos = ipos + (npts(i) + 1)/2 END DO ier = 1 RETURN ! EVALUATE SUM OF WT(I) * FUNC(X(I)) 20 nby2 = (npt + 1)/2 DO i = 1, nby2 x = xpt(ipos) fn_val = fn_val + wt(ipos) * (func(x) + func(-x)) ipos = ipos + 1 END DO ier = 0 RETURN END FUNCTION hermite