! Some useful snippets

!-----------------------------------------------------------------------

! Create double precision parameter 'pi' that is hardware and

! compiler independent [and two alternative formulations, but use first]

      DOUBLE PRECISION, PARAMETER :: pi  = 4.d0 * ATAN(1.d0)

      DOUBLE PRECISION, PARAMETER :: pi1 = ACOS(-1.d0)

      DOUBLE PRECISION, PARAMETER :: pi2 = 2.d0 * ASIN(1.d0)

!-----------------------------------------------------------------------

! Cube root function [vs. x**(1.d0/3.d0)]

      DOUBLE PRECISION FUNCTION CubeRt (x)

      IMPLICIT none

      DOUBLE PRECISION x

      INTEGER i

      IF (x .EQ. 0.d0) THEN ; CubeRt = 0.d0 ; RETURN ; END IF

      CubeRt = SIGN(1.d0, x) * EXP(LOG(ABS(x))/3.d0)

      DO i = 1, 2

        CubeRt = CubeRt - (CubeRt**3 - x) / (3.d0 * CubeRt**2)

      END DO

      END