MODULE aux CONTAINS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FUNCTION up_bd(x,P2,Q2,Pint) !This function returns the value of the upper boundary !(pressure value) for a given energy density value. IMPLICIT NONE REAL(KIND=8) :: x, P2(:), Q2(:), Pint(:) REAL(KIND=8) :: up_bd IF( x .LT. Pint(1) ) THEN up_bd = x + P2(2) - P2(1) ELSE up_bd = Q2(2) END IF END FUNCTION up_bd !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FUNCTION lo_bd(x,P1,Q1,Pint) !This function returns the value of the lower boundary !(pressure value) for a given energy density value. IMPLICIT NONE REAL(KIND=8) :: x, P1(:), Q1(:), Pint(:) REAL(KIND=8) :: lo_bd IF( x .LT. Pint(1) ) THEN lo_bd = P1(2) ELSE lo_bd = x + Q1(2) - Q1(1) END IF END FUNCTION lo_bd !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FUNCTION check_constraints(pt1, pt2) !This function returns T if, for two given points p1 and p2, monotony and !causality constraints are satisfied. IMPLICIT NONE REAL(KIND=8) :: pt1(:), pt2(:) LOGICAL ::check_constraints IF( (pt2(2) .GE. pt1(2)) .AND. (pt2(2)-pt1(2))/(pt2(1)-pt1(1)) .LE. 1D0 ) THEN check_constraints = .TRUE. ELSE check_constraints = .FALSE. END IF !PRINT *, (pt2(2)-pt1(2))/(pt2(1)-pt1(1)) !PRINT *, pt1(:) !PRINT *, pt2(:) END FUNCTION check_constraints !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE spline(x, y, n, yp1, ypn, y2) !This routine has been taken from Numerical Recipes in Fortran 90, p109 !Given arrays x(1:n) and y(1:n) containing a tabulated function. i.e. yi=f(xi), !with x1