maths.f90

Go to the documentation of this file.

MODULE maths

  USE base_routines
  USE constants
  USE kinds
  USE iso_varying_string
  USE strings

#include "macros.h"  
  
  IMPLICIT NONE

  PRIVATE

  !Module parameters

  !Module types

  !Interfaces

  INTERFACE cross_product
    MODULE PROCEDURE crossproductintg
    MODULE PROCEDURE crossproductsp
    MODULE PROCEDURE crossproductdp
  END INTERFACE cross_product

  INTERFACE crossproduct
    MODULE PROCEDURE crossproductintg
    MODULE PROCEDURE crossproductsp
    MODULE PROCEDURE crossproductdp
  END INTERFACE crossproduct

  INTERFACE d_cross_product
    MODULE PROCEDURE dcrossproductintg
    MODULE PROCEDURE dcrossproductsp
    MODULE PROCEDURE dcrossproductdp
  END INTERFACE d_cross_product

  INTERFACE dcrossproduct
    MODULE PROCEDURE dcrossproductintg
    MODULE PROCEDURE dcrossproductsp
    MODULE PROCEDURE dcrossproductdp
  END INTERFACE dcrossproduct

  INTERFACE determinant
    MODULE PROCEDURE determinantfullintg
    MODULE PROCEDURE determinantfullsp
    MODULE PROCEDURE determinantfulldp
  END INTERFACE determinant
        
  INTERFACE edp
    MODULE PROCEDURE edpsp
    MODULE PROCEDURE edpdp
  END INTERFACE edp

  INTERFACE eigenvalue
    MODULE PROCEDURE eigenvaluefullsp
    MODULE PROCEDURE eigenvaluefulldp
  END INTERFACE eigenvalue

  INTERFACE eigenvector
    MODULE PROCEDURE eigenvectorfullsp
    MODULE PROCEDURE eigenvectorfulldp
  END INTERFACE eigenvector

  INTERFACE i0
    MODULE PROCEDURE i0dp
    MODULE PROCEDURE i0sp
  END INTERFACE i0

  INTERFACE i1
    MODULE PROCEDURE i1dp
    MODULE PROCEDURE i1sp
  END INTERFACE i1

  INTERFACE invert
    MODULE PROCEDURE invertfullsp
    MODULE PROCEDURE invertfulldp
  END INTERFACE invert

  INTERFACE k0
    MODULE PROCEDURE k0dp
    MODULE PROCEDURE k0sp
  END INTERFACE k0

  INTERFACE k1
    MODULE PROCEDURE k1dp
    MODULE PROCEDURE k1sp
  END INTERFACE k1

  INTERFACE kdp
    MODULE PROCEDURE kdpdp
    MODULE PROCEDURE kdpsp
  END INTERFACE kdp

  INTERFACE identitymatrix
    MODULE PROCEDURE identitymatrixsp
    MODULE PROCEDURE identitymatrixdp
  END INTERFACE identitymatrix

  INTERFACE l2norm
    MODULE PROCEDURE l2normsp
    MODULE PROCEDURE l2normdp
  END INTERFACE l2norm
    
  INTERFACE matrix_product
    MODULE PROCEDURE matrixproductsp
    MODULE PROCEDURE matrixproductdp
  END INTERFACE matrix_product
  
  INTERFACE matrixproduct
    MODULE PROCEDURE matrixproductsp
    MODULE PROCEDURE matrixproductdp
  END INTERFACE matrixproduct
  
  INTERFACE matrixtransposeproduct
    MODULE PROCEDURE matrixtransposeproductsp
    MODULE PROCEDURE matrixtransposeproductdp
  END INTERFACE matrixtransposeproduct
  
  INTERFACE matrixproducttranspose
    MODULE PROCEDURE matrixproducttransposesp
    MODULE PROCEDURE matrixproducttransposedp
  END INTERFACE matrixproducttranspose
  
  INTERFACE matrix_transpose
    MODULE PROCEDURE matrixtransposesp
    MODULE PROCEDURE matrixtransposedp
  END INTERFACE matrix_transpose

  INTERFACE matrixtranspose
    MODULE PROCEDURE matrixtransposesp
    MODULE PROCEDURE matrixtransposedp
  END INTERFACE matrixtranspose

  INTERFACE matrix_vector_product
    MODULE PROCEDURE matrixvectorproductsp
    MODULE PROCEDURE matrixvectorproductdp
  END INTERFACE matrix_vector_product

  INTERFACE matrixvectorproduct
    MODULE PROCEDURE matrixvectorproductsp
    MODULE PROCEDURE matrixvectorproductdp
  END INTERFACE matrixvectorproduct
  
  INTERFACE matrixtransposevectorproduct
    MODULE PROCEDURE matrixtransposevectorproductsp
    MODULE PROCEDURE matrixtransposevectorproductdp
  END INTERFACE matrixtransposevectorproduct
  
  INTERFACE normalise
    MODULE PROCEDURE normalisesp
    MODULE PROCEDURE normalisedp
  END INTERFACE normalise

  INTERFACE norm_cross_product
    MODULE PROCEDURE normcrossproductsp
    MODULE PROCEDURE normcrossproductdp
  END INTERFACE norm_cross_product

  INTERFACE normcrossproduct
    MODULE PROCEDURE normcrossproductsp
    MODULE PROCEDURE normcrossproductdp
  END INTERFACE normcrossproduct

  INTERFACE solve_small_linear_system
    MODULE PROCEDURE solvesmalllinearsystemsp
    MODULE PROCEDURE solvesmalllinearsystemdp
  END INTERFACE solve_small_linear_system

  INTERFACE solvesmalllinearsystem
    MODULE PROCEDURE solvesmalllinearsystemsp
    MODULE PROCEDURE solvesmalllinearsystemdp
  END INTERFACE solvesmalllinearsystem

  INTERFACE coth
    MODULE PROCEDURE cothsp
    MODULE PROCEDURE cothdp
  END INTERFACE coth 

  PUBLIC cross_product,crossproduct,d_cross_product,dcrossproduct,determinant,eigenvalue,eigenvector,identitymatrix,invert, &
    & l2norm,matrix_product,matrixproduct,matrixtransposeproduct,matrixproducttranspose,matrix_transpose,matrixtranspose, &
    & normalise,norm_cross_product,normcrossproduct,solve_small_linear_system,solvesmalllinearsystem,coth,spline_cubic_set, &
    & s3_fs,spline_cubic_val,matrix_vector_product,matrixvectorproduct,matrixtransposevectorproduct
  
  
CONTAINS

  !
  !================================================================================================================================
  !

  SUBROUTINE crossproductintg(a,b,c,err,error,*)
      
    !Argument variables
    INTEGER(INTG), INTENT(IN) :: a(:)
    INTEGER(INTG), INTENT(IN) :: b(:)
    INTEGER(INTG), INTENT(OUT) :: c(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("CrossProductIntg",err,error,*999)

    IF(SIZE(a,1)==SIZE(b,1)) THEN
      IF(SIZE(c,1)==3) THEN
        SELECT CASE(SIZE(a,1))
        CASE(3)
          c(1)=a(2)*b(3)-a(3)*b(2)
          c(2)=a(3)*b(1)-a(1)*b(3)
          c(3)=a(1)*b(2)-a(2)*b(1)
        CASE DEFAULT
          CALL flagerror("Invalid vector size.",err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("The vector c is not the correct size.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("The vectors a and b are not the same size.",err,error,*999)
    ENDIF

    exits("CrossProductIntg")
    RETURN
999 errorsexits("CrossProductIntg",err,error)
    RETURN 1
    
  END SUBROUTINE crossproductintg
  
  !
  !================================================================================================================================
  !

  SUBROUTINE crossproductsp(a,b,c,err,error,*)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: a(:)
    REAL(SP), INTENT(IN) :: b(:)
    REAL(SP), INTENT(OUT) :: c(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("CrossProductSP",err,error,*999)

    IF(SIZE(a,1)==SIZE(b,1)) THEN
      IF(SIZE(c,1)==3) THEN
        SELECT CASE(SIZE(a,1))
        CASE(3)
          c(1)=a(2)*b(3)-a(3)*b(2)
          c(2)=a(3)*b(1)-a(1)*b(3)
          c(3)=a(1)*b(2)-a(2)*b(1)
        CASE DEFAULT
          CALL flagerror("Invalid vector size.",err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("The vector c is not the correct size.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("The vectors a and b are not the same size.",err,error,*999)
    ENDIF

    exits("CrossProductSP")
    RETURN
999 errorsexits("CrossProductSP",err,error)
    RETURN 1
    
  END SUBROUTINE crossproductsp
  
  !
  !================================================================================================================================
  !

  SUBROUTINE crossproductdp(a,b,c,err,error,*)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: a(:)
    REAL(DP), INTENT(IN) :: b(:)
    REAL(DP), INTENT(OUT) :: c(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("CrossProductDP",err,error,*999)

    IF(SIZE(a,1)==SIZE(b,1)) THEN
      IF(SIZE(c,1)==3) THEN
        SELECT CASE(SIZE(a,1))
        CASE(3)
          c(1)=a(2)*b(3)-a(3)*b(2)
          c(2)=a(3)*b(1)-a(1)*b(3)
          c(3)=a(1)*b(2)-a(2)*b(1)
        CASE DEFAULT
          CALL flagerror("Invalid vector size.",err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("The vector C is not the correct size.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("The vectors A and B are not the same size.",err,error,*999)
    ENDIF

    exits("CrossProductDP")
    RETURN
999 errorsexits("CrossProductDP",err,error)
    RETURN 1
    
  END SUBROUTINE crossproductdp
  
  !
  !================================================================================================================================
  !
  
  SUBROUTINE dcrossproductintg(n,a,b,c,da,db,dc,err,error,*)
    
    !Argument variables
    INTEGER(INTG), INTENT(IN) :: n
    INTEGER(INTG), INTENT(IN) :: a(:)
    INTEGER(INTG), INTENT(IN) :: b(:)
    INTEGER(INTG), INTENT(OUT) :: c(:)
    INTEGER(INTG), INTENT(IN) :: da(:,:)
    INTEGER(INTG), INTENT(IN) :: db(:,:)
    INTEGER(INTG), INTENT(OUT) :: dc(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: derivIdx
    
    enters("dCrossProductIntg",err,error,*999)

    CALL crossproduct(a,b,c,err,error,*999)
    IF(SIZE(da,1)==SIZE(db,1).AND.SIZE(a,1)==SIZE(da,1).AND.SIZE(b,1)==SIZE(db,1)) THEN
      IF(SIZE(da,2)==n.AND.SIZE(db,2)==n) THEN
        IF(SIZE(c,1)==3) THEN
          SELECT CASE(SIZE(a,1))
          CASE(3)
            DO derividx=1,n
              dc(1,derividx)=da(2,derividx)*b(3)-da(3,derividx)*b(2)+a(2)*db(3,derividx)-a(3)*db(2,derividx)
              dc(2,derividx)=da(3,derividx)*b(1)-da(1,derividx)*b(3)+a(3)*db(1,derividx)-a(1)*db(3,derividx)
              dc(3,derividx)=da(1,derividx)*b(2)-da(2,derividx)*b(1)+a(1)*db(2,derividx)-a(2)*db(1,derividx)
            ENDDO !derivIdx
          CASE DEFAULT
            CALL flagerror("Invalid vector size.",err,error,*999)
          END SELECT
        ELSE
          CALL flagerror("The vector c is not the correct size.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The number of derivative vectors is too small.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("The vectors for da and db are not the same size.",err,error,*999)
    ENDIF

    exits("dCrossProductIntg")
    RETURN
999 errorsexits("dCrossProductIntg",err,error)
    RETURN 1
    
  END SUBROUTINE dcrossproductintg
  
  !
  !================================================================================================================================
  !
  
  SUBROUTINE dcrossproductsp(n,a,b,c,da,db,dc,err,error,*)
  
    !Argument variables
    INTEGER(INTG), INTENT(IN) :: n
    REAL(SP), INTENT(IN) :: a(:)
    REAL(SP), INTENT(IN) :: b(:)
    REAL(SP), INTENT(OUT) :: c(:)
    REAL(SP), INTENT(IN) :: da(:,:)
    REAL(SP), INTENT(IN) :: db(:,:)
    REAL(SP), INTENT(OUT) :: dc(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: derivIdx
    
    enters("dCrossProductSP",err,error,*999)

    CALL crossproduct(a,b,c,err,error,*999)
    IF(SIZE(da,1)==SIZE(db,1).AND.SIZE(a,1)==SIZE(da,1).AND.SIZE(b,1)==SIZE(db,1)) THEN
      IF(SIZE(da,2)==n.AND.SIZE(db,2)==n) THEN
        IF(SIZE(c,1)==3) THEN
          SELECT CASE(SIZE(a,1))
          CASE(3)
            DO derividx=1,n
              dc(1,derividx)=da(2,derividx)*b(3)-da(3,derividx)*b(2)+a(2)*db(3,derividx)-a(3)*db(2,derividx)
              dc(2,derividx)=da(3,derividx)*b(1)-da(1,derividx)*b(3)+a(3)*db(1,derividx)-a(1)*db(3,derividx)
              dc(3,derividx)=da(1,derividx)*b(2)-da(2,derividx)*b(1)+a(1)*db(2,derividx)-a(2)*db(1,derividx)
            ENDDO !derivIdx
          CASE DEFAULT
            CALL flagerror("Invalid vector size.",err,error,*999)
          END SELECT
        ELSE
          CALL flagerror("The vector c is not the correct size.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The number of derivative vectors is too small.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("The vectors for da and db are not the same size.",err,error,*999)
    ENDIF

    exits("dCrossProductSP")
    RETURN
999 errorsexits("dCrossProductSP",err,error)
    RETURN 1
    
  END SUBROUTINE dcrossproductsp
  
  !
  !================================================================================================================================
  !

  SUBROUTINE dcrossproductdp(n,a,b,c,da,db,dc,err,error,*)
  
    !Argument variables
    INTEGER(INTG), INTENT(IN) :: n
    REAL(DP), INTENT(IN) :: a(:)
    REAL(DP), INTENT(IN) :: b(:)
    REAL(DP), INTENT(OUT) :: c(:)
    REAL(DP), INTENT(IN) :: da(:,:)
    REAL(DP), INTENT(IN) :: db(:,:)
    REAL(DP), INTENT(OUT) :: dc(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: derivIdx
    
    enters("dCrossProductDP",err,error,*999)

    CALL crossproduct(a,b,c,err,error,*999)
    IF(SIZE(da,1)==SIZE(db,1).AND.SIZE(a,1)==SIZE(da,1).AND.SIZE(b,1)==SIZE(db,1)) THEN
      IF(SIZE(da,2)==n.AND.SIZE(db,2)==n) THEN
        IF(SIZE(c,1)==3) THEN
          SELECT CASE(SIZE(a,1))
          CASE(3)
            DO derividx=1,n
              dc(1,derividx)=da(2,derividx)*b(3)-da(3,derividx)*b(2)+a(2)*db(3,derividx)-a(3)*db(2,derividx)
              dc(2,derividx)=da(3,derividx)*b(1)-da(1,derividx)*b(3)+a(3)*db(1,derividx)-a(1)*db(3,derividx)
              dc(3,derividx)=da(1,derividx)*b(2)-da(2,derividx)*b(1)+a(1)*db(2,derividx)-a(2)*db(1,derividx)
            ENDDO !derivIdx
          CASE DEFAULT
            CALL flagerror("Invalid vector size.",err,error,*999)
          END SELECT
        ELSE
          CALL flagerror("The vector c is not the correct size.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The number of derivative vectors is too small.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("The vectors for da and db are not the same size.",err,error,*999)
    ENDIF

    exits("dCrossProductDP")
    RETURN
999 errorsexits("dCrossProductDP",err,error)
    RETURN 1
    
  END SUBROUTINE dcrossproductdp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixvectorproductsp(A,b,c,err,error,*)

    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(IN) :: b(:)
    REAL(SP), INTENT(OUT) :: c(:)
    INTEGER(INTG) :: err
    TYPE(varying_string), INTENT(OUT) :: error

    enters("MatrixVectorProductSP",err,error,*999)

    IF(SIZE(a,2)==SIZE(b,1).AND.SIZE(a,1)==SIZE(c,1)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        c(1)=a(1,1)*b(1)
      CASE(2)
        c(1)=a(1,1)*b(1)+a(1,2)*b(2)
        c(2)=a(2,1)*b(1)+a(2,2)*b(2)
      CASE(3)
        c(1)=a(1,1)*b(1)+a(1,2)*b(2)+a(1,3)*b(3)
        c(2)=a(2,1)*b(1)+a(2,2)*b(2)+a(2,3)*b(3)
        c(3)=a(3,1)*b(1)+a(3,2)*b(2)+a(3,3)*b(3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix and/or vector size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix sizes.",err,error,*999)
    ENDIF
    
    exits("MatrixVectorProductSP")
    RETURN
999 errorsexits("MatrixVectorProductSP",err,error)
    RETURN 1
    
  END SUBROUTINE matrixvectorproductsp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixvectorproductdp(A,b,c,err,error,*)

    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(IN) :: B(:)
    REAL(DP), INTENT(OUT) :: C(:)
    INTEGER(INTG) :: err
    TYPE(varying_string), INTENT(OUT) :: error

    enters("MatrixVectorProductDP",err,error,*999)

    IF(SIZE(a,2)==SIZE(b,1).AND.SIZE(a,1)==SIZE(c,1)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        c(1)=a(1,1)*b(1)
      CASE(2)
        c(1)=a(1,1)*b(1)+a(1,2)*b(2)
        c(2)=a(2,1)*b(1)+a(2,2)*b(2)
      CASE(3)
        c(1)=a(1,1)*b(1)+a(1,2)*b(2)+a(1,3)*b(3)
        c(2)=a(2,1)*b(1)+a(2,2)*b(2)+a(2,3)*b(3)
        c(3)=a(3,1)*b(1)+a(3,2)*b(2)+a(3,3)*b(3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix and/or vector size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix sizes.",err,error,*999)
    ENDIF
    
    exits("MatrixVectorProductDP")
    RETURN
999 errorsexits("MatrixVectorProductDP",err,error)
    RETURN 1
     
  END SUBROUTINE matrixvectorproductdp
  
  !
  !================================================================================================================================
  !

  SUBROUTINE matrixtransposevectorproductsp(A,b,c,err,error,*)

    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(IN) :: b(:)
    REAL(SP), INTENT(OUT) :: c(:)
    INTEGER(INTG) :: err
    TYPE(varying_string), INTENT(OUT) :: error

    enters("MatrixTransposeVectorProductSP",err,error,*999)

    IF(SIZE(a,1)==SIZE(b,1).AND.SIZE(a,2)==SIZE(c,1)) THEN
       SELECT CASE(SIZE(a,2))
       CASE(1)
         c(1)=a(1,1)*b(1)
       CASE(2)
         c(1)=a(1,1)*b(1)+a(2,1)*b(2)
         c(2)=a(1,2)*b(1)+a(2,2)*b(2)
       CASE(3)
         c(1)=a(1,1)*b(1)+a(2,1)*b(2)+a(3,1)*b(3)
         c(2)=a(1,2)*b(1)+a(2,2)*b(2)+a(3,2)*b(3)
         c(3)=a(1,3)*b(1)+a(2,3)*b(2)+a(3,3)*b(3)
       CASE DEFAULT
         CALL flagerror("Invalid matrix and/or vector size.",err,error,*999)
       END SELECT
     ELSE
       CALL flagerror("Invalid matrix sizes.",err,error,*999)
     ENDIF

     exits("MatrixTransposeVectorProductSP")
     RETURN
 999 errorsexits("MatrixTransposeVectorProductSP",err,error)
     RETURN 1
     
  END SUBROUTINE matrixtransposevectorproductsp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixtransposevectorproductdp(A,b,c,err,error,*)

    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(IN) :: b(:)
    REAL(DP), INTENT(OUT) :: c(:)
    INTEGER(INTG) :: err
    TYPE(varying_string), INTENT(OUT) :: error

    enters("MatrixTransposeVectorProductDP",err,error,*999)

    IF(SIZE(a,1)==SIZE(b,1).AND.SIZE(a,2)==SIZE(c,1)) THEN
       SELECT CASE(SIZE(a,2))
       CASE(1)
         c(1)=a(1,1)*b(1)
       CASE(2)
         c(1)=a(1,1)*b(1)+a(2,1)*b(2)
         c(2)=a(1,2)*b(1)+a(2,2)*b(2)
       CASE(3)
         c(1)=a(1,1)*b(1)+a(2,1)*b(2)+a(3,1)*b(3)
         c(2)=a(1,2)*b(1)+a(2,2)*b(2)+a(3,2)*b(3)
         c(3)=a(1,3)*b(1)+a(2,3)*b(2)+a(3,3)*b(3)
       CASE DEFAULT
         CALL flagerror("Invalid matrix and/or vector size.",err,error,*999)
       END SELECT
     ELSE
       CALL flagerror("Invalid matrix sizes.",err,error,*999)
     ENDIF

     exits("MatrixTransposeVectorProductDP")
     RETURN
 999 errorsexits("MatrixTransposeVectorProductDP",err,error)
     RETURN 1
     
  END SUBROUTINE matrixtransposevectorproductdp

  !
  !================================================================================================================================
  !

  FUNCTION determinantfullintg(A,err,error)
  
    !Argument variables
    INTEGER(INTG), INTENT(IN) :: A(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Function variable
    INTEGER(INTG) :: DeterminantFullIntg
    
    enters("DeterminantFullIntg",err,error,*999)

    determinantfullintg=0
    
    IF(SIZE(a,1)==SIZE(a,2)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        determinantfullintg=a(1,1)
      CASE(2)
        determinantfullintg=a(1,1)*a(2,2)-a(2,1)*a(1,2)
      CASE(3)
        determinantfullintg=a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+a(1,3)*a(2,1)*a(3,2)-a(1,1)*a(3,2)*a(2,3)- &
          a(2,1)*a(1,2)*a(3,3)-a(3,1)*a(2,2)*a(1,3)
      CASE DEFAULT
        CALL flagerror("Matrix size not implemented.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Matrix is not square.",err,error,*999)
    ENDIF

    exits("DeterminantFullIntg")
    RETURN
999 errorsexits("DeterminantFullIntg",err,error)
    RETURN
    
  END FUNCTION determinantfullintg
  
  !
  !================================================================================================================================
  !

  FUNCTION determinantfullsp(A,err,error)
    
    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Function variable
    REAL(SP) :: DeterminantFullSP

    enters("DeterminantFullSP",err,error,*999)

    determinantfullsp=0.0_sp
    
    IF(SIZE(a,1)==SIZE(a,2)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        determinantfullsp=a(1,1)
      CASE(2)
        determinantfullsp=a(1,1)*a(2,2)-a(2,1)*a(1,2)
      CASE(3)
        determinantfullsp=a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+a(1,3)*a(2,1)*a(3,2)-a(1,1)*a(3,2)*a(2,3)- &
          a(2,1)*a(1,2)*a(3,3)-a(3,1)*a(2,2)*a(1,3)
      CASE DEFAULT
        CALL flagerror("Matrix size not implemented.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Matrix is not square.",err,error,*999)
    ENDIF

    exits("DeterminantFullSP")
    RETURN
999 errorsexits("DeterminantFullSP",err,error)
    RETURN
    
  END FUNCTION determinantfullsp
  
  !
  !================================================================================================================================
  !

  FUNCTION determinantfulldp(A,err,error)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Function variable
    REAL(DP) :: DeterminantFullDP
    
    enters("DeterminantFullDP",err,error,*999)

    determinantfulldp=0.0_dp
    
    IF(SIZE(a,1)==SIZE(a,2)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        determinantfulldp=a(1,1)
      CASE(2)
        determinantfulldp=a(1,1)*a(2,2)-a(2,1)*a(1,2)
      CASE(3)
        determinantfulldp=a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+a(1,3)*a(2,1)*a(3,2)-a(1,1)*a(3,2)*a(2,3)- &
          a(2,1)*a(1,2)*a(3,3)-a(3,1)*a(2,2)*a(1,3)
      CASE DEFAULT
        CALL flagerror("Matrix size not implemented.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Matrix is not square.",err,error,*999)
    ENDIF

    exits("DeterminantFullDP")
    RETURN
999 errorsexits("DeterminantFullDP",err,error)
    RETURN
    
  END FUNCTION determinantfulldp

  !
  !================================================================================================================================
  !

  PURE FUNCTION edpdp(x)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: x
    !Function variable
    REAL(DP) :: EdpDP
    !Local variables
    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    REAL(DP), PARAMETER :: a1=0.44325141463_dp
    REAL(DP), PARAMETER :: a2=0.06260601220_dp
    REAL(DP), PARAMETER :: a3=0.04757383546_dp
    REAL(DP), PARAMETER :: a4=0.01736506451_dp
    REAL(DP), PARAMETER :: b1=0.24998368310_dp
    REAL(DP), PARAMETER :: b2=0.09200180037_dp
    REAL(DP), PARAMETER :: b3=0.04069697526_dp
    REAL(DP), PARAMETER :: b4=0.00526449639_dp
    REAL(DP) :: term1,term2,x1
    
    x1=1.0_dp-x
    term1=1.0_dp+(a1+(a2+(a3+a4*x1)*x1)*x1)*x1
    term2=(b1+(b2+(b3+b4*x1)*x1)*x1)*x1
    edpdp=term1+term2*log(1.0_dp/x1)
    
    RETURN 
  END FUNCTION edpdp
  
  !
  !================================================================================================================================
  !

  PURE FUNCTION edpsp(x)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: x
    !Function variable
    REAL(SP) :: EdpSP
    !Local variables
    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    REAL(SP), PARAMETER :: a1=0.44325141463_sp
    REAL(SP), PARAMETER :: a2=0.06260601220_sp
    REAL(SP), PARAMETER :: a3=0.04757383546_sp
    REAL(SP), PARAMETER :: a4=0.01736506451_sp
    REAL(SP), PARAMETER :: b1=0.24998368310_sp
    REAL(SP), PARAMETER :: b2=0.09200180037_sp
    REAL(SP), PARAMETER :: b3=0.04069697526_sp
    REAL(SP), PARAMETER :: b4=0.00526449639_sp
    REAL(SP) :: term1,term2,x1
    
    x1=1.0_sp-x
    term1=1.0_sp+(a1+(a2+(a3+a4*x1)*x1)*x1)*x1
    term2=(b1+(b2+(b3+b4*x1)*x1)*x1)*x1
    edpsp=term1+term2*log(1.0_sp/x1)
        
    RETURN
    
  END FUNCTION edpsp
  
  !
  !================================================================================================================================
  !

  SUBROUTINE eigenvaluefullsp(A,eValues,err,error,*)
    
    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(OUT) :: eValues(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: i
    REAL(SP) :: angle,b2,b3,c1,c2,d,q,q3,r,ri1,ri2,ri3,ri4,rq,temp,theta
    
    enters("EigenvalueFullSP",err,error,*999)

    IF(SIZE(a,1)==SIZE(a,2)) THEN
      IF(SIZE(a,1)<=SIZE(evalues,1)) THEN
        SELECT CASE(SIZE(a,1))
        CASE(1)
          evalues(1)=a(1,1)
        CASE(2)
          IF(abs(a(1,2))>zero_tolerance_sp) THEN
            ri1=a(1,1)+a(2,2)
            ri2=a(1,1)*a(2,2)-a(1,2)**2
            b2=ri1/2.0_sp
            c1=ri1*ri1
            c2=4.0_sp*ri2
            IF(c2>c1) CALL flagerror("Complex roots found in quadratic equation.",err,error,*999)
            b3=sqrt(c1-c2)/2.0_sp
            evalues(1)=b2+b3
            evalues(2)=b2-b3
          ELSE
            evalues(1)=a(1,1)
            evalues(2)=a(2,2)
          ENDIF
          IF(abs(evalues(2))>abs(evalues(1))) THEN
            temp=evalues(1)
            evalues(1)=evalues(2)
            evalues(2)=temp
          ENDIF
        CASE(3)
          ri1=a(1,1)+a(2,2)+a(3,3)
          ri2=a(1,1)*a(2,2)+a(2,2)*a(3,3)+a(3,3)*a(1,1)-(a(1,2)**2+a(2,3)**2+a(3,1)**2)
          ri3=determinant(a,err,error)
          IF(err/=0) GOTO 999
          ri4=ri1/3.0_sp
          q=ri4*ri4-ri2/3.0_sp   
          r=ri4*(ri4*ri4-ri2/2.0_sp)+ri3/2.0_sp
          q3=q*q*q
          d=r*r-q3
          IF(abs(d)>zero_tolerance_sp) CALL flagerror("Complex roots found in solution of cubic equation.",err,error,*999)
          rq=sqrt(abs(q))
          IF(abs(q)<zero_tolerance_sp) THEN
            theta=0.0_sp
          ELSE
            theta=acos(r/sqrt(abs(q3)))/3.0_sp
          ENDIF
          angle=2.0_sp*REAL(pi,sp)/3.0_SP
          evalues(1)=2.0_sp*rq*cos(theta)+ri4
          evalues(2)=2.0_sp*rq*cos(theta+angle)+ri4
          evalues(3)=2.0_sp*rq*cos(theta+2.0_sp*angle)+ri4
          DO i=1,2
            IF(abs(evalues(3))>abs(evalues(i))) THEN
              temp=evalues(i)
              evalues(i)=evalues(3)
              evalues(3)=temp
            ENDIF
          ENDDO !i
        CASE DEFAULT
          CALL flagerror("Matrix size not implemented.",err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Evalues is too small.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not square.",err,error,*999)
    ENDIF

    exits("EigenvalueFullSP")
    RETURN
999 errorsexits("EigenvalueFullSP",err,error)
    RETURN 1
    
  END SUBROUTINE eigenvaluefullsp

  !
  !================================================================================================================================
  !

  SUBROUTINE eigenvaluefulldp(A,eValues,err,error,*)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(OUT) :: eValues(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: i
    REAL(DP) :: angle,b2,b3,c1,c2,d,q,q3,r,ri1,ri2,ri3,ri4,rq,temp,theta
    
    enters("EigenvalueFullDP",err,error,*999)

    IF(SIZE(a,1)==SIZE(a,2)) THEN
      IF(SIZE(a,1)<= SIZE(evalues,1)) THEN
        SELECT CASE(SIZE(a,1))
        CASE(1)
          evalues(1)=a(1,1)
        CASE(2)
          IF(abs(a(1,2))>zero_tolerance_dp) THEN
            ri1=a(1,1)+a(2,2)
            ri2=a(1,1)*a(2,2)-a(1,2)**2
            b2=ri1/2.0_dp
            c1=ri1*ri1
            c2=4.0_dp*ri2
            IF(c2>c1) CALL flagerror("Complex roots found in quadratic equation.",err,error,*999)
            b3=sqrt(c1-c2)/2.0_dp
            evalues(1)=b2+b3
            evalues(2)=b2-b3
          ELSE
            evalues(1)=a(1,1)
            evalues(2)=a(2,2)
          ENDIF
          IF(abs(evalues(2))>abs(evalues(1))) THEN
            temp=evalues(1)
            evalues(1)=evalues(2)
            evalues(2)=temp
          ENDIF
        CASE(3)
          ri1=a(1,1)+a(2,2)+a(3,3)
          ri2=a(1,1)*a(2,2)+a(2,2)*a(3,3)+a(3,3)*a(1,1)-(a(1,2)**2+a(2,3)**2+a(3,1)**2)
          ri3=determinant(a,err,error)
          IF(err/=0) GOTO 999
          ri4=ri1/3.0_dp
          q=ri4*ri4-ri2/3.0_dp   
          r=ri4*(ri4*ri4-ri2/2.0_dp)+ri3/2.0_dp
          q3=q*q*q
          d=r*r-q3
          IF(abs(d)>zero_tolerance_dp) CALL flagerror("Complex roots found in solution of cubic equation.",err,error,*999)
          rq=sqrt(abs(q))
          IF(abs(q)<zero_tolerance_dp) THEN
            theta=0.0_dp
          ELSE
            theta=acos(r/sqrt(abs(q3)))/3.0_dp
          ENDIF
          angle=2.0_dp*pi/3.0_dp
          evalues(1)=2.0_dp*rq*cos(theta)+ri4
          evalues(2)=2.0_dp*rq*cos(theta+angle)+ri4
          evalues(3)=2.0_dp*rq*cos(theta+2.0_dp*angle)+ri4
          DO i=1,2
            IF(abs(evalues(3))>abs(evalues(i))) THEN
              temp=evalues(i)
              evalues(i)=evalues(3)
              evalues(3)=temp
            ENDIF
          ENDDO !i
        CASE DEFAULT
          CALL flagerror("Matrix size not implemented.",err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Evalues is too small.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not square.",err,error,*999)
    ENDIF

    exits("EigenvalueFullDP")
    RETURN
999 errorsexits("EigenvalueFullDP",err,error)
    RETURN 1
    
  END SUBROUTINE eigenvaluefulldp

  !
  !================================================================================================================================
  !

  SUBROUTINE eigenvectorfullsp(A,eValue,eVector,err,error,*)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(IN) :: eValue
    REAL(SP), INTENT(OUT) :: eVector(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: i,i1,i2,i3,iCycle(3,3)
    REAL(SP) :: al,b(size(a,1)),sum,U(size(a,1),size(a,2)),x(size(a,1))

    DATA icycle /2,1,1,3,1,2,1,2,1/
    
    enters("EigenvectorFullSP",err,error,*999)

!!THIS NEEDS TO BE CHECKED
    
    IF(SIZE(a,1)==SIZE(a,2)) THEN
      IF(SIZE(a,1)<=SIZE(evector,1)) THEN
        SELECT CASE(SIZE(a,1))
        CASE(1)
          evector(1)=1.0_sp
        CASE(2)
          IF(abs(a(1,2))>zero_tolerance_sp) THEN
            IF(abs(a(1,1)-evalue)>abs(a(2,2)-evalue)) THEN
              al=sqrt(a(1,2)**2+(a(1,1)-evalue)**2)
              evector(1)=a(1,2)/al
              evector(2)=(evalue-a(1,1))/al
            ELSE
              al=sqrt(a(1,2)**2+(a(2,2)-evalue)**2)
              evector(1)=(evalue-a(2,2))/al
              evector(2)=a(1,2)/al
            ENDIF
          ELSE IF(abs(evalue-a(1,1))<zero_tolerance_sp) THEN
            evector(1)=1.0_sp
            evector(2)=0.0_sp
          ELSE IF(abs(evalue-a(2,2))<zero_tolerance_sp) THEN
            evector(1)=0.0_sp
            evector(2)=1.0_dp
          ENDIF
        CASE(3)
          IF(abs(a(1,2))<zero_tolerance_sp.AND.abs(a(1,3))<zero_tolerance_sp.AND.abs(a(2,3))<zero_tolerance_sp) THEN
            evector=0.0_sp
            CALL flagerror("Zero matrix?? Eigenvectors undetermined.",err,error,*999)
          ELSE
            DO i=1,3
              u(i,:)=a(i,:)
              u(i,i)=u(i,i)-evalue
            ENDDO !i
            DO i=1,3
              x(i)=1.0_sp
              i1=icycle(i,1)
              i2=icycle(i,2)
              i3=icycle(i,3)
              b(1)=-1.0_sp*u(i1,i)
              b(2)=-1.0_sp*u(i2,i)
              CALL solvesmalllinearsystem(u(i1:i2:i3,i1:i2:i3),x,b,err,error,*999)
              sum=dot_product(u(i,:),x)
              IF(abs(sum)<zero_tolerance_sp) THEN
                evector=normalise(x,err,error)
                IF(err /= 0) GOTO 999
              ENDIF
            ENDDO !i
          ENDIF
        CASE DEFAULT
          CALL flagerror("Matrix size not implemented.",err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Evector is too small.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not square.",err,error,*999)
    ENDIF

    exits("EigenvectorFullSP")
    RETURN
999 errorsexits("EigenvectorFullSP",err,error)
    RETURN 1
    
  END SUBROUTINE eigenvectorfullsp

  !
  !================================================================================================================================
  !

  SUBROUTINE eigenvectorfulldp(A,eValue,eVector,err,error,*)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(IN) :: eValue
    REAL(DP), INTENT(OUT) :: eVector(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: i,i1,i2,i3,iCycle(3,3)
    REAL(DP) :: al,b(size(a,1)),sum,U(size(a,1),size(a,2)),x(size(a,1))

    DATA icycle /2,1,1,3,1,2,1,2,1/
    
    enters("EigenvectorFullDP",err,error,*999)

!!THIS NEEDS TO BE CHECKED

    IF(SIZE(a,1)==SIZE(a,2)) THEN
      IF(SIZE(a,1)<=SIZE(evector,1)) THEN
        SELECT CASE(SIZE(a,1))
        CASE(1)
          evector(1)=1.0_dp
        CASE(2)
          IF(abs(a(1,2))>zero_tolerance_dp) THEN
            IF(abs(a(1,1)-evalue)>abs(a(2,2)-evalue)) THEN
              al=sqrt(a(1,2)**2+(a(1,1)-evalue)**2)
              evector(1)=a(1,2)/al
              evector(2)=(evalue-a(1,1))/al
            ELSE
              al=sqrt(a(1,2)**2+(a(2,2)-evalue)**2)
              evector(1)=(evalue-a(2,2))/al
              evector(2)=a(1,2)/al
            ENDIF
          ELSE IF(abs(evalue-a(1,1))<zero_tolerance_dp) THEN
            evector(1)=1.0_dp
            evector(2)=0.0_dp
          ELSE IF(abs(evalue-a(2,2))<zero_tolerance_dp) THEN
            evector(1)=0.0_dp
            evector(2)=1.0_dp
          ENDIF
        CASE(3)
          IF(abs(a(1,2))<zero_tolerance_dp.AND.abs(a(1,3))<zero_tolerance_dp.AND.abs(a(2,3))<zero_tolerance_dp) THEN
            evector=0.0_dp
            CALL flagerror("Zero matrix?? Eigenvectors undetermined.",err,error,*999)
          ELSE
            DO i=1,3
              u(i,:)=a(i,:)
              u(i,i)=u(i,i)-evalue
            ENDDO !i
            DO i=1,3
              x(i)=1.0_dp
              i1=icycle(i,1)
              i2=icycle(i,2)
              i3=icycle(i,3)
              b(1)=-1.0_dp*u(i1,i)
              b(2)=-1.0_dp*u(i2,i)
              CALL solvesmalllinearsystem(u(i1:i2:i3,i1:i2:i3),x,b,err,error,*999)
              sum=dot_product(u(i,:),x)
              IF(abs(sum)<zero_tolerance_dp) THEN
                evector=normalise(x,err,error)
                IF(err /= 0) GOTO 999
              ENDIF
            ENDDO !i
          ENDIF
        CASE DEFAULT
          CALL flagerror("Matrix size not implemented.",err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Evector is too small.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not square.",err,error,*999)
    ENDIF

    exits("EigenvectorFullDP")
    RETURN
999 errorsexits("EigenvectorFullDP",err,error)
    RETURN 1
    
  END SUBROUTINE eigenvectorfulldp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION i0dp(x)
      
    !Argument variables
    REAL(DP), INTENT(IN) :: x
    !Function variable
    REAL(DP) :: I0DP
    !Local variables
    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    REAL(DP), PARAMETER :: a1=1.0_dp
    REAL(DP), PARAMETER :: a2=3.5156229_dp
    REAL(DP), PARAMETER :: a3=3.0899424_dp
    REAL(DP), PARAMETER :: a4=1.2067492_dp
    REAL(DP), PARAMETER :: a5=0.2659732_dp
    REAL(DP), PARAMETER :: a6=0.0360768_dp
    REAL(DP), PARAMETER :: a7=0.0045813_dp
    REAL(DP) :: t

    !Calculate I0(x) for x < 3.75
    t=x*x/(3.75_dp*3.75_dp)
    i0dp=a1+(a2+(a3+(a4+(a5+(a6+a7*t)*t)*t)*t)*t)*t
    
    RETURN
    
  END FUNCTION i0dp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION i0sp(x)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: x
    !Function variable
    REAL(SP) :: I0SP
    !Local variables
    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    REAL(SP), PARAMETER :: a1=1.0_sp
    REAL(SP), PARAMETER :: a2=3.5156229_sp
    REAL(SP), PARAMETER :: a3=3.0899424_sp
    REAL(SP), PARAMETER :: a4=1.2067492_sp
    REAL(SP), PARAMETER :: a5=0.2659732_sp
    REAL(SP), PARAMETER :: a6=0.0360768_sp
    REAL(SP), PARAMETER :: a7=0.0045813_sp
    REAL(SP) :: t

    !Calculate I0(x) for x < 3.75
    t=x*x/(3.75_sp*3.75_sp)
    i0sp=a1+(a2+(a3+(a4+(a5+(a6+a7*t)*t)*t)*t)*t)*t
    
    RETURN
    
  END FUNCTION i0sp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION i1dp(x)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: x
    !Function variable
    REAL(DP) :: I1DP
    !Local variables
    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    REAL(DP), PARAMETER :: a1=0.5_dp
    REAL(DP), PARAMETER :: a2=0.87890594_dp
    REAL(DP), PARAMETER :: a3=0.51498869_dp
    REAL(DP), PARAMETER :: a4=0.15084934_dp
    REAL(DP), PARAMETER :: a5=0.02658733_dp
    REAL(DP), PARAMETER :: a6=0.00301532_dp
    REAL(DP), PARAMETER :: a7=0.00032411_dp
    REAL(DP) :: t

    !Calculate I1(x)
    t=(x/3.75_dp)*(x/3.75_dp)
    i1dp=(a1+(a2+(a3+(a4+(a5+(a6+a7*t)*t)*t)*t)*t)*t)*x
    
    RETURN
    
  END FUNCTION i1dp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION i1sp(x)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: x
    !Function variable
    REAL(SP) :: I1SP
    !Local variables
    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    REAL(SP), PARAMETER :: a1=0.5_sp
    REAL(SP), PARAMETER :: a2=0.87890594_sp
    REAL(SP), PARAMETER :: a3=0.51498869_sp
    REAL(SP), PARAMETER :: a4=0.15084934_sp
    REAL(SP), PARAMETER :: a5=0.02658733_sp
    REAL(SP), PARAMETER :: a6=0.00301532_sp
    REAL(SP), PARAMETER :: a7=0.00032411_sp
    REAL(SP) :: t

    !Calculate I1(x)
    t=(x/3.75_sp)*(x/3.75_sp)
    i1sp=(a1+(a2+(a3+(a4+(a5+(a6+7*t)*t)*t)*t)*t)*t)*x
    
    RETURN
    
  END FUNCTION i1sp
  
  !
  !================================================================================================================================
  !

  SUBROUTINE identitymatrixsp(A,err,error,*)
    
    !Argument variables
    REAL(SP), INTENT(OUT) :: A(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: i
    
    enters("IdentityMatrixSP",err,error,*999)

    IF(SIZE(a,1)==SIZE(a,2)) THEN
      SELECT CASE(SIZE(a,1)) 
      CASE(1)
        a(1,1)=1.0_sp
      CASE(2)
        a(1,1)=1.0_sp
        a(2,1)=0.0_sp
        a(1,2)=0.0_sp
        a(2,2)=1.0_sp
      CASE(3)
        a(1,1)=1.0_sp
        a(2,1)=0.0_sp
        a(3,1)=0.0_sp
        a(1,2)=0.0_sp
        a(2,2)=1.0_sp
        a(3,2)=0.0_sp
        a(1,3)=0.0_sp
        a(2,3)=0.0_sp
        a(3,3)=1.0_sp
      CASE DEFAULT
        a=0.0_dp
        DO i=1,SIZE(a,1)
          a(i,i)=1.0_sp
        ENDDO !i
      END SELECT
    ELSE
      CALL flagerror("Matrix A is not square.",err,error,*999)
    ENDIF

    exits("IdentityMatrixSP")
    RETURN
999 errorsexits("IdentityMatrixSP",err,error)
    RETURN 1
    
  END SUBROUTINE identitymatrixsp

  !
  !================================================================================================================================
  !

  SUBROUTINE identitymatrixdp(A,err,error,*)
    
    !Argument variables
    REAL(DP), INTENT(OUT) :: A(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: i
    
    enters("IdentityMatrixDP",err,error,*999)

    IF(SIZE(a,1)==SIZE(a,2)) THEN
      SELECT CASE(SIZE(a,1)) 
      CASE(1)
        a(1,1)=1.0_dp
      CASE(2)
        a(1,1)=1.0_dp
        a(2,1)=0.0_dp
        a(1,2)=0.0_dp
        a(2,2)=1.0_dp
      CASE(3)
        a(1,1)=1.0_dp
        a(2,1)=0.0_dp
        a(3,1)=0.0_dp
        a(1,2)=0.0_dp
        a(2,2)=1.0_dp
        a(3,2)=0.0_dp
        a(1,3)=0.0_dp
        a(2,3)=0.0_dp
        a(3,3)=1.0_dp
      CASE DEFAULT
        a=0.0_dp
        DO i=1,SIZE(a,1)
          a(i,i)=1.0_dp
        ENDDO !i
      END SELECT
    ELSE
      CALL flagerror("Matrix A is not square.",err,error,*999)
    ENDIF

    exits("IdentityMatrixDP")
    RETURN
999 errorsexits("IdentityMatrixDP",err,error)
    RETURN 1
    
  END SUBROUTINE identitymatrixdp

  !
  !================================================================================================================================
  !

  SUBROUTINE invertfullsp(A,B,det,err,error,*)
    
    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(OUT) :: B(:,:)
    REAL(SP), INTENT(OUT) :: det
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables

    enters("InvertFullSP",err,error,*999)

    IF(SIZE(a,1)==SIZE(a,2)) THEN
      IF(SIZE(b,1)==SIZE(a,1).AND.SIZE(b,2)==SIZE(a,2)) THEN
        SELECT CASE(SIZE(a,1)) 
        CASE(1)
          det=a(1,1)
          IF(abs(det)>zero_tolerance_sp) THEN
            b(1,1)=1.0_sp/a(1,1)
          ELSE
            CALL flagwarning("Matrix A is zero and cannot be inverted.",err,error,*999)
            b(1,1)=0.0_dp
          ENDIF
        CASE(2)
          det=a(1,1)*a(2,2)-a(1,2)*a(2,1)
          IF(abs(det)>zero_tolerance_sp) THEN
            b(1,1)=a(2,2)/det
            b(1,2)=-a(1,2)/det
            b(2,1)=-a(2,1)/det
            b(2,2)=a(1,1)/det
          ELSE
            CALL flagwarning("Zero Determinant for matrix A.",err,error,*999)
            b=0.0_dp
          ENDIF
        CASE(3)
          det=a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+a(1,3)*a(2,1)*a(3,2)-a(1,1)*a(3,2)*a(2,3)-a(2,1)*a(1,2)*a(3,3)- &
            & a(3,1)*a(2,2)*a(1,3)
          IF(abs(det)>zero_tolerance_sp) THEN
            b(1,1)=(a(2,2)*a(3,3)-a(3,2)*a(2,3))/det
            b(2,1)=(a(2,3)*a(3,1)-a(3,3)*a(2,1))/det
            b(3,1)=(a(2,1)*a(3,2)-a(3,1)*a(2,2))/det
            b(1,2)=(a(3,2)*a(1,3)-a(1,2)*a(3,3))/det
            b(2,2)=(a(3,3)*a(1,1)-a(1,3)*a(3,1))/det
            b(3,2)=(a(3,1)*a(1,2)-a(1,1)*a(3,2))/det
            b(1,3)=(a(1,2)*a(2,3)-a(2,2)*a(1,3))/det
            b(2,3)=(a(1,3)*a(2,1)-a(2,3)*a(1,1))/det
            b(3,3)=(a(1,1)*a(2,2)-a(2,1)*a(1,2))/det
          ELSE
            CALL flagwarning("Zero Determinant for matrix A.",err,error,*999)
            b=0.0_dp
          ENDIF
        CASE DEFAULT
          CALL flagerror("Matrix size is not implemented.",err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Matrix B is not the same size as matrix A.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix A is not square.",err,error,*999)
    ENDIF

    exits("InvertFullSP")
    RETURN
999 errorsexits("InvertFullSP",err,error)
    RETURN 1
    
  END SUBROUTINE invertfullsp

  !
  !================================================================================================================================
  !

  SUBROUTINE invertfulldp(A,B,det,err,error,*)
    
    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(OUT) :: B(:,:)
    REAL(DP), INTENT(OUT) :: det
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables

    enters("InvertFullDP",err,error,*999)

    IF(SIZE(a,1)==SIZE(a,2)) THEN
      IF(SIZE(b,1)==SIZE(a,1).AND.SIZE(b,2)==SIZE(a,2)) THEN
        SELECT CASE(SIZE(a,1)) 
        CASE(1)
          det=a(1,1)
          IF(abs(det)>zero_tolerance_dp) THEN
            b(1,1)=1.0_dp/a(1,1)
          ELSE
            CALL flagwarning("Matrix A is zero and cannot be inverted",err,error,*999)
            b(1,1)=0.0_dp
          ENDIF
        CASE(2)
          det=a(1,1)*a(2,2)-a(2,1)*a(1,2)
          IF(abs(det)>zero_tolerance_dp) THEN
            b(1,1)=a(2,2)/det
            b(1,2)=-a(1,2)/det
            b(2,1)=-a(2,1)/det
            b(2,2)=a(1,1)/det
          ELSE
            CALL flagwarning("Zero Determinant for matrix A.",err,error,*999)
            b=0.0_dp
          ENDIF
        CASE(3)
          det=a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+a(1,3)*a(2,1)*a(3,2)-a(1,1)*a(3,2)*a(2,3)-a(2,1)*a(1,2)*a(3,3)- &
            & a(3,1)*a(2,2)*a(1,3)
          IF(abs(det)>zero_tolerance_dp) THEN
            b(1,1)=(a(2,2)*a(3,3)-a(3,2)*a(2,3))/det
            b(2,1)=(a(2,3)*a(3,1)-a(3,3)*a(2,1))/det
            b(3,1)=(a(2,1)*a(3,2)-a(3,1)*a(2,2))/det
            b(1,2)=(a(3,2)*a(1,3)-a(1,2)*a(3,3))/det
            b(2,2)=(a(3,3)*a(1,1)-a(1,3)*a(3,1))/det
            b(3,2)=(a(3,1)*a(1,2)-a(1,1)*a(3,2))/det
            b(1,3)=(a(1,2)*a(2,3)-a(2,2)*a(1,3))/det
            b(2,3)=(a(1,3)*a(2,1)-a(2,3)*a(1,1))/det
            b(3,3)=(a(1,1)*a(2,2)-a(2,1)*a(1,2))/det
          ELSE
            CALL flagwarning("Zero Determinant for matrix A.",err,error,*999)
            b=0.0_dp
          ENDIF
        CASE DEFAULT
          CALL flagerror("Matrix size is not implemented.",err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Matrix B is not the same size as matrix A.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix A is not square.",err,error,*999)
    ENDIF

    exits("InvertFullDP")
    RETURN
999 errorsexits("InvertFullDP",err,error)
    RETURN 1
    
  END SUBROUTINE invertfulldp
  
  !
  !================================================================================================================================
  !

  PURE FUNCTION k0dp(x)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: x
    !Function variable
    REAL(DP) :: K0DP
    !Local variables
    REAL(DP) :: a1,a2,a3,a4,a5,a6,a7,t

    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    !Calculate K0(x)
    IF(x<=2.0_dp) THEN
      t=x*x/4.0_dp
      a1=-0.57721566_dp
      a2=0.42278420_dp
      a3=0.23069756_dp
      a4=0.03488590_dp
      a5=0.00262698_dp
      a6=0.00010750_dp
      a7=0.00000740_dp
      k0dp=-log(x/2.0_dp)*i0(x)+(a1+(a2+(a3+(a4+(a5+(a6+a7*t)*t)*t)*t)*t)*t)
    ELSE
      IF(x>174.0_dp) THEN
        k0dp=0.0_dp
      ELSE
        t=2.0_dp/x
        a1=1.25331414_dp
        a2=-0.07832358_dp
        a3=0.02189568_dp
        a4=-0.01062446_dp
        a5=0.00587872_dp
        a6=-0.00251540_dp
        a7=0.00053208_dp
        k0dp=(a1+(a2+(a3+(a4+(a5+(a6+a7*t)*t)*t)*t)*t)*t)*exp(-x)/sqrt(x)
      ENDIF
    ENDIF
    
    RETURN
    
  END FUNCTION k0dp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION k0sp(x)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: x
    !Function variable
    REAL(SP) :: K0SP
    !Local variables
    REAL(SP) :: a1,a2,a3,a4,a5,a6,a7,t

    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    !Calculate K0(x)
    IF(x<=2.0_sp) THEN
      t=x*x/4.0_sp
      a1=-0.57721566_sp
      a2=0.42278420_sp
      a3=0.23069756_sp
      a4=0.03488590_sp
      a5=0.00262698_sp
      a6=0.00010750_sp
      a7=0.00000740_sp
      k0sp=-log(x/2.0_sp)*i0(x)+(a1+(a2+(a3+(a4+(a5+(a6+a7*t)*t)*t)*t)*t)*t)
    ELSE
      IF(x>174.0_sp) THEN
        k0sp=0.0_sp
      ELSE
        t=2.0_sp/x
        a1=1.25331414_sp
        a2=-0.07832358_sp
        a3=0.02189568_sp
        a4=-0.01062446_sp
        a5=0.00587872_sp
        a6=-0.00251540_sp
        a7=0.00053208_sp
        k0sp=(a1+(a2+(a3+(a4+(a5+(a6+a7*t)*t)*t)*t)*t)*t)*exp(-x)/sqrt(x)
      ENDIF
    ENDIF
    
    RETURN
    
  END FUNCTION k0sp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION k1dp(x)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: x
    !Function variable
    REAL(DP) :: K1DP
    !Local variables
    REAL(DP) :: a1,a2,a3,a4,a5,a6,a7,a8,t

    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    !Calculate K1(x)
    IF(x<=2.0_dp) THEN
      t=(x/2.0_dp)*(x/2.0_dp)
      a1=log(x/2.0_dp)*i1(x)
      a2=1.0_dp
      a3=0.15443144_dp
      a4=-0.67278579_dp
      a5=-0.18156897_dp
      a6=-0.01919402_dp
      a7=-0.00110404_dp
      a8=-0.00004686_dp
      k1dp=a1+a2/x+(a3+(a4+(a5+(a6+(a7+a8*t)*t)*t)*t)*t)*x/4.0_dp
    ELSE
      IF (x>174.0_dp) THEN
        k1dp=0.0_dp
      ELSE
        t=2.0_dp/x
        a1=1.25331414_dp
        a2=0.23498619_dp
        a3=-0.03655620_dp
        a4=0.01504268_dp
        a5=-0.00780355_dp
        a6=0.00325614_dp
        a7=-0.00068245_dp
        k1dp=(a1+(a2+(a3+(a4+(a5+(a6+a7*t)*t)*t)*t)*t)*t)*exp(-x)/sqrt(x)
      ENDIF
    ENDIF
    
    RETURN
    
  END FUNCTION k1dp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION k1sp(x)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: x
    !Function variable
    REAL(SP) :: K1SP
    !Local variables
    REAL(SP) :: a1,a2,a3,a4,a5,a6,a7,a8,t

    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    !Calculate K1(x)
    IF(x<=2.0_sp) THEN
      t=(x/2.0_sp)*(x/2.0_sp)
      a1=log(x/2.0_sp)*i1(x)
      a2=1.0_sp
      a3=0.15443144_sp
      a4=-0.67278579_sp
      a5=-0.18156897_sp
      a6=-0.01919402_sp
      a7=-0.00110404_sp
      a8=-0.00004686_sp
      k1sp=a1+a2/x+(a3+(a4+(a5+(a6+(a7+a8*t)*t)*t)*t)*t)*x/4.0_sp
    ELSE
      IF (x>174.0_sp) THEN
        k1sp=0.0_sp
      ELSE
        t=2.0_sp/x
        a1=1.25331414_sp
        a2=0.23498619_sp
        a3=-0.03655620_sp
        a4=0.01504268_sp
        a5=-0.00780355_sp
        a6=0.00325614_sp
        a7=-0.00068245_sp
        k1sp=(a1+(a2+(a3+(a4+(a5+(a6+a7*t)*t)*t)*t)*t)*t)*exp(-x)/sqrt(x)
      ENDIF
    ENDIF
    
    RETURN
    
  END FUNCTION k1sp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION kdpdp(x)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: x
    !Function variable
    REAL(DP) :: KdpDP
    !Local variables
    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    REAL(DP), PARAMETER :: a0=1.38629436112_dp
    REAL(DP), PARAMETER :: a1=0.09666344259_dp
    REAL(DP), PARAMETER :: a2=0.03590092383_dp
    REAL(DP), PARAMETER :: a3=0.03742563713_dp
    REAL(DP), PARAMETER :: a4=0.01451196212_dp
    REAL(DP), PARAMETER :: b0=0.5_dp
    REAL(DP), PARAMETER :: b1=0.12498593597_dp
    REAL(DP), PARAMETER :: b2=0.06880248576_dp
    REAL(DP), PARAMETER :: b3=0.03328355346_dp
    REAL(DP), PARAMETER :: b4=0.00441787012_dp
    REAL(DP) :: term1,term2,x1

    x1=1.0_dp-x
    term1=a0+(a1+(a2+(a3+a4*x)*x)*x)*x
    term2=b0+(b1+(b2+(b3+b4*x)*x)*x)*x
    kdpdp=term1+term2*log(1.0_dp/x)    
    
    RETURN
    
  END FUNCTION kdpdp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION kdpsp(x)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: x
    !Function variable
    REAL(SP) :: KdpSP
    !Local variables
    !!TODO:: CHECK PRECISION ON THESE CONSTANTS
    REAL(SP), PARAMETER :: a0=1.38629436112_sp
    REAL(SP), PARAMETER :: a1=0.09666344259_sp
    REAL(SP), PARAMETER :: a2=0.03590092383_sp
    REAL(SP), PARAMETER :: a3=0.03742563713_sp
    REAL(SP), PARAMETER :: a4=0.01451196212_sp
    REAL(SP), PARAMETER :: b0=0.5_sp
    REAL(SP), PARAMETER :: b1=0.12498593597_sp
    REAL(SP), PARAMETER :: b2=0.06880248576_sp
    REAL(SP), PARAMETER :: b3=0.03328355346_sp
    REAL(SP), PARAMETER :: b4=0.00441787012_sp
    REAL(SP) :: term1,term2,x1

    x1=1.0_sp-x
    term1=a0+(a1+(a2+(a3+a4*x)*x)*x)*x
    term2=b0+(b1+(b2+(b3+b4*x)*x)*x)*x
    kdpsp=term1+term2*log(1.0_sp/x)    
    
    RETURN
    
  END FUNCTION kdpsp
  
  !
  !================================================================================================================================
  !
  
  PURE FUNCTION l2normsp(a)

    !Argument variables
    REAL(SP), INTENT(IN) :: a(:)
    !Function variable
    REAL(SP) :: L2NormSP
    !Local variables
    REAL(SP) :: asum
    
    asum=sum(a*a,1)
    l2normsp=sqrt(asum)

    RETURN
    
  END FUNCTION l2normsp

  !
  !================================================================================================================================
  !

  FUNCTION l2normdp(A)

    !Argument variables
    REAL(DP), INTENT(IN) :: a(:)
    !Function variable
    REAL(DP) :: L2NormDP
    !Local variables
    REAL(DP) :: asum

    asum=sum(a*a,1)
    l2normdp=sqrt(asum)

    RETURN
    
  END FUNCTION l2normdp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixproductsp(A,B,C,err,error,*)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(IN) :: B(:,:)
    REAL(SP), INTENT(OUT) :: C(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("MatrixProductSP",err,error,*999)

    IF(SIZE(a,2)==SIZE(b,1).AND.SIZE(a,1)==SIZE(c,1).AND.SIZE(b,2)==SIZE(c,2)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        c(1,1)=a(1,1)*b(1,1)
      CASE(2)
        c(1,1)=a(1,1)*b(1,1)+a(1,2)*b(2,1)
        c(1,2)=a(1,1)*b(1,2)+a(1,2)*b(2,2)
        c(2,1)=a(2,1)*b(1,1)+a(2,2)*b(2,1)
        c(2,2)=a(2,1)*b(1,2)+a(2,2)*b(2,2)
      CASE(3)
        c(1,1)=a(1,1)*b(1,1)+a(1,2)*b(2,1)+a(1,3)*b(3,1)
        c(1,2)=a(1,1)*b(1,2)+a(1,2)*b(2,2)+a(1,3)*b(3,2)
        c(1,3)=a(1,1)*b(1,3)+a(1,2)*b(2,3)+a(1,3)*b(3,3)        
        c(2,1)=a(2,1)*b(1,1)+a(2,2)*b(2,1)+a(2,3)*b(3,1)
        c(2,2)=a(2,1)*b(1,2)+a(2,2)*b(2,2)+a(2,3)*b(3,2)
        c(2,3)=a(2,1)*b(1,3)+a(2,2)*b(2,3)+a(2,3)*b(3,3)
        c(3,1)=a(3,1)*b(1,1)+a(3,2)*b(2,1)+a(3,3)*b(3,1)
        c(3,2)=a(3,1)*b(1,2)+a(3,2)*b(2,2)+a(3,3)*b(3,2)
        c(3,3)=a(3,1)*b(1,3)+a(3,2)*b(2,3)+a(3,3)*b(3,3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix sizes.",err,error,*999)
    ENDIF

    exits("MatrixProductSP")
    RETURN
999 errorsexits("MatrixProductSP",err,error)
    RETURN 1
    
  END SUBROUTINE matrixproductsp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixproductdp(A,B,C,err,error,*)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(IN) :: B(:,:)
    REAL(DP), INTENT(OUT) :: C(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
        
    enters("MatrixProductDP",err,error,*999)
    
   IF(SIZE(a,2)==SIZE(b,1).AND.SIZE(a,1)==SIZE(c,1).AND.SIZE(b,2)==SIZE(c,2)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        c(1,1)=a(1,1)*b(1,1)
      CASE(2)
        c(1,1)=a(1,1)*b(1,1)+a(1,2)*b(2,1)
        c(1,2)=a(1,1)*b(1,2)+a(1,2)*b(2,2)
        c(2,1)=a(2,1)*b(1,1)+a(2,2)*b(2,1)
        c(2,2)=a(2,1)*b(1,2)+a(2,2)*b(2,2)
      CASE(3)
        c(1,1)=a(1,1)*b(1,1)+a(1,2)*b(2,1)+a(1,3)*b(3,1)
        c(1,2)=a(1,1)*b(1,2)+a(1,2)*b(2,2)+a(1,3)*b(3,2)
        c(1,3)=a(1,1)*b(1,3)+a(1,2)*b(2,3)+a(1,3)*b(3,3)        
        c(2,1)=a(2,1)*b(1,1)+a(2,2)*b(2,1)+a(2,3)*b(3,1)
        c(2,2)=a(2,1)*b(1,2)+a(2,2)*b(2,2)+a(2,3)*b(3,2)
        c(2,3)=a(2,1)*b(1,3)+a(2,2)*b(2,3)+a(2,3)*b(3,3)
        c(3,1)=a(3,1)*b(1,1)+a(3,2)*b(2,1)+a(3,3)*b(3,1)
        c(3,2)=a(3,1)*b(1,2)+a(3,2)*b(2,2)+a(3,3)*b(3,2)
        c(3,3)=a(3,1)*b(1,3)+a(3,2)*b(2,3)+a(3,3)*b(3,3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix sizes.",err,error,*999)
    ENDIF

    exits("MatrixProductDP")
    RETURN
999 errorsexits("MatrixProductDP",err,error)
    RETURN 1
    
  END SUBROUTINE matrixproductdp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixtransposeproductsp(A,B,C,err,error,*)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(IN) :: B(:,:)
    REAL(SP), INTENT(OUT) :: C(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("MatrixTransposeProductSP",err,error,*999)

    IF(SIZE(a,2)==SIZE(b,1).AND.SIZE(a,1)==SIZE(c,1).AND.SIZE(b,2)==SIZE(c,2)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        c(1,1)=a(1,1)*b(1,1)
      CASE(2)
        c(1,1)=a(1,1)*b(1,1)+a(2,1)*b(2,1)
        c(1,2)=a(1,1)*b(1,2)+a(2,1)*b(2,2)
        c(2,1)=a(1,2)*b(1,1)+a(2,2)*b(2,1)
        c(2,2)=a(1,2)*b(1,2)+a(2,2)*b(2,2)
      CASE(3)
        c(1,1)=a(1,1)*b(1,1)+a(2,1)*b(2,1)+a(3,1)*b(3,1)
        c(1,2)=a(1,1)*b(1,2)+a(2,1)*b(2,2)+a(3,1)*b(3,2)
        c(1,3)=a(1,1)*b(1,3)+a(2,1)*b(2,3)+a(3,1)*b(3,3)        
        c(2,1)=a(1,2)*b(1,1)+a(2,2)*b(2,1)+a(3,2)*b(3,1)
        c(2,2)=a(1,2)*b(1,2)+a(2,2)*b(2,2)+a(3,2)*b(3,2)
        c(2,3)=a(1,2)*b(1,3)+a(2,2)*b(2,3)+a(3,2)*b(3,3)
        c(3,1)=a(1,3)*b(1,1)+a(2,3)*b(2,1)+a(3,3)*b(3,1)
        c(3,2)=a(1,3)*b(1,2)+a(2,3)*b(2,2)+a(3,3)*b(3,2)
        c(3,3)=a(1,3)*b(1,3)+a(2,3)*b(2,3)+a(3,3)*b(3,3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix sizes.",err,error,*999)
    ENDIF

    exits("MatrixTransposeProductSP")
    RETURN
999 errorsexits("MatrixTransposeProductSP",err,error)
    RETURN 1
    
  END SUBROUTINE matrixtransposeproductsp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixtransposeproductdp(A,B,C,err,error,*)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(IN) :: B(:,:)
    REAL(DP), INTENT(OUT) :: C(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("MatrixTransposeProductDP",err,error,*999)

    IF(SIZE(a,2)==SIZE(b,1).AND.SIZE(a,1)==SIZE(c,1).AND.SIZE(b,2)==SIZE(c,2)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        c(1,1)=a(1,1)*b(1,1)
      CASE(2)
        c(1,1)=a(1,1)*b(1,1)+a(2,1)*b(2,1)
        c(1,2)=a(1,1)*b(1,2)+a(2,1)*b(2,2)
        c(2,1)=a(1,2)*b(1,1)+a(2,2)*b(2,1)
        c(2,2)=a(1,2)*b(1,2)+a(2,2)*b(2,2)
      CASE(3)
        c(1,1)=a(1,1)*b(1,1)+a(2,1)*b(2,1)+a(3,1)*b(3,1)
        c(1,2)=a(1,1)*b(1,2)+a(2,1)*b(2,2)+a(3,1)*b(3,2)
        c(1,3)=a(1,1)*b(1,3)+a(2,1)*b(2,3)+a(3,1)*b(3,3)        
        c(2,1)=a(1,2)*b(1,1)+a(2,2)*b(2,1)+a(3,2)*b(3,1)
        c(2,2)=a(1,2)*b(1,2)+a(2,2)*b(2,2)+a(3,2)*b(3,2)
        c(2,3)=a(1,2)*b(1,3)+a(2,2)*b(2,3)+a(3,2)*b(3,3)
        c(3,1)=a(1,3)*b(1,1)+a(2,3)*b(2,1)+a(3,3)*b(3,1)
        c(3,2)=a(1,3)*b(1,2)+a(2,3)*b(2,2)+a(3,3)*b(3,2)
        c(3,3)=a(1,3)*b(1,3)+a(2,3)*b(2,3)+a(3,3)*b(3,3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix sizes.",err,error,*999)
    ENDIF

    exits("MatrixTransposeProductDP")
    RETURN
999 errorsexits("MatrixTransposeProductDP",err,error)
    RETURN 1
    
  END SUBROUTINE matrixtransposeproductdp

  !
  !================================================================================================================================
  !
  SUBROUTINE matrixproducttransposesp(A,B,C,err,error,*)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(IN) :: B(:,:)
    REAL(SP), INTENT(OUT) :: C(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("MatrixProductTransposeSP",err,error,*999)

    IF(SIZE(a,2)==SIZE(b,2).AND.SIZE(a,1)==SIZE(c,1).AND.SIZE(b,1)==SIZE(c,2)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        c(1,1)=a(1,1)*b(1,1)
      CASE(2)
        c(1,1)=a(1,1)*b(1,1)+a(1,2)*b(1,2)
        c(1,2)=a(1,1)*b(2,1)+a(1,2)*b(2,2)
        c(2,1)=a(2,1)*b(1,1)+a(2,2)*b(1,2)
        c(2,2)=a(2,1)*b(2,1)+a(2,2)*b(2,2)
      CASE(3)
        c(1,1)=a(1,1)*b(1,1)+a(1,2)*b(1,2)+a(1,3)*b(1,3)
        c(1,2)=a(1,1)*b(2,1)+a(1,2)*b(2,2)+a(1,3)*b(2,3)
        c(1,3)=a(1,1)*b(3,1)+a(1,2)*b(3,2)+a(1,3)*b(3,3)        
        c(2,1)=a(2,1)*b(1,1)+a(2,2)*b(1,2)+a(2,3)*b(1,3)
        c(2,2)=a(2,1)*b(2,1)+a(2,2)*b(2,2)+a(2,3)*b(2,3)
        c(2,3)=a(2,1)*b(3,1)+a(2,2)*b(3,2)+a(2,3)*b(3,3)
        c(3,1)=a(3,1)*b(1,1)+a(3,2)*b(1,2)+a(3,3)*b(1,3)
        c(3,2)=a(3,1)*b(2,1)+a(3,2)*b(2,2)+a(3,3)*b(2,3)
        c(3,3)=a(3,1)*b(3,1)+a(3,2)*b(3,2)+a(3,3)*b(3,3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix sizes.",err,error,*999)
    ENDIF

    exits("MatrixProductTransposeSP")
    RETURN
999 errorsexits("MatrixProductTransposeSP",err,error)
    RETURN 1
    
  END SUBROUTINE matrixproducttransposesp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixproducttransposedp(A,B,C,err,error,*)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(IN) :: B(:,:)
    REAL(DP), INTENT(OUT) :: C(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
        
    enters("MatrixProductTranposeDP",err,error,*999)
    
   IF(SIZE(a,2)==SIZE(b,2).AND.SIZE(a,1)==SIZE(c,1).AND.SIZE(b,1)==SIZE(c,2)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        c(1,1)=a(1,1)*b(1,1)
      CASE(2)
        c(1,1)=a(1,1)*b(1,1)+a(1,2)*b(1,2)
        c(1,2)=a(1,1)*b(2,1)+a(1,2)*b(2,2)
        c(2,1)=a(2,1)*b(1,1)+a(2,2)*b(1,2)
        c(2,2)=a(2,1)*b(2,1)+a(2,2)*b(2,2)
      CASE(3)
        c(1,1)=a(1,1)*b(1,1)+a(1,2)*b(1,2)+a(1,3)*b(1,3)
        c(1,2)=a(1,1)*b(2,1)+a(1,2)*b(2,2)+a(1,3)*b(2,3)
        c(1,3)=a(1,1)*b(3,1)+a(1,2)*b(3,2)+a(1,3)*b(3,3)        
        c(2,1)=a(2,1)*b(1,1)+a(2,2)*b(1,2)+a(2,3)*b(1,3)
        c(2,2)=a(2,1)*b(2,1)+a(2,2)*b(2,2)+a(2,3)*b(2,3)
        c(2,3)=a(2,1)*b(3,1)+a(2,2)*b(3,2)+a(2,3)*b(3,3)
        c(3,1)=a(3,1)*b(1,1)+a(3,2)*b(1,2)+a(3,3)*b(1,3)
        c(3,2)=a(3,1)*b(2,1)+a(3,2)*b(2,2)+a(3,3)*b(2,3)
        c(3,3)=a(3,1)*b(3,1)+a(3,2)*b(3,2)+a(3,3)*b(3,3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix sizes.",err,error,*999)
    ENDIF

    exits("MatrixProductTransposeDP")
    RETURN
999 errorsexits("MatrixProductTransposeDP",err,error)
    RETURN 1
    
  END SUBROUTINE matrixproducttransposedp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixtransposesp(A,AT,err,error,*)
          
    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(OUT) :: AT(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("MatrixTransposeSP",err,error,*999)

    IF(SIZE(a,1)==SIZE(at,2).AND.SIZE(a,2)==SIZE(at,1)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        at(1,1)=a(1,1)
      CASE(2)
        at(1,1)=a(1,1)
        at(1,2)=a(2,1)
        at(2,1)=a(1,2)
        at(2,2)=a(2,2)
      CASE(3)
        at(1,1)=a(1,1)
        at(1,2)=a(2,1)
        at(1,3)=a(3,1)
        at(2,1)=a(1,2)
        at(2,2)=a(2,2)
        at(2,3)=a(3,2)
        at(3,1)=a(1,3)
        at(3,2)=a(2,3)
        at(3,3)=a(3,3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix size.",err,error,*999)
    ENDIF
 
    exits("MatrixTransposeSP")
    RETURN
999 errorsexits("MatrixTransposeSP",err,error)
    RETURN 1
    
  END SUBROUTINE matrixtransposesp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrixtransposedp(A,AT,err,error,*)
    
    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(OUT) :: AT(:,:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
        
    enters("MatrixTransposeDP",err,error,*999)

    IF(SIZE(a,1)==SIZE(at,2).AND.SIZE(a,2)==SIZE(at,1)) THEN
      SELECT CASE(SIZE(a,1))
      CASE(1)
        at(1,1)=a(1,1)
      CASE(2)
        at(1,1)=a(1,1)
        at(1,2)=a(2,1)
        at(2,1)=a(1,2)
        at(2,2)=a(2,2)
      CASE(3)
        at(1,1)=a(1,1)
        at(1,2)=a(2,1)
        at(1,3)=a(3,1)
        at(2,1)=a(1,2)
        at(2,2)=a(2,2)
        at(2,3)=a(3,2)
        at(3,1)=a(1,3)
        at(3,2)=a(2,3)
        at(3,3)=a(3,3)
      CASE DEFAULT
        CALL flagerror("Invalid matrix size.",err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Invalid matrix size.",err,error,*999)
    ENDIF
    
    exits("MatrixTransposeDP")
    RETURN
999 errorsexits("MatrixTransposeDP",err,error)
    RETURN 1
    
  END SUBROUTINE matrixtransposedp

  !
  !================================================================================================================================
  !
  
  FUNCTION normalisesp(a,err,error)

    !Argument variables
    REAL(SP), INTENT(IN) :: a(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Function variable
    REAL(SP) :: NormaliseSP(size(a,1))
    !Local variables
    REAL(SP) :: length
    
    enters("",err,error,*999)

    length=l2norm(a)
    IF(abs(length)<zero_tolerance_sp) THEN
      normalisesp=a
      CALL flagerror("Length of the vector to normalise is zero.",err,error,*999)
    ELSE
      normalisesp=a/length
    ENDIF

    exits("NormaliseSP")
    RETURN
999 errorsexits("NormaliseSP",err,error)
    RETURN
    
  END FUNCTION normalisesp

  !
  !================================================================================================================================
  !

  FUNCTION normalisedp(a,err,error)

    !Argument variables
    REAL(DP), INTENT(IN) :: a(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Function variable
    REAL(DP) :: NormaliseDP(size(a,1))
    !Local variables
    REAL(DP) :: length
    
    enters("NormaliseDP",err,error,*999)

    length=l2norm(a)
    IF(abs(length)<zero_tolerance_dp) THEN
      normalisedp=a
      CALL flagerror("Length of the vector to normalise is zero.",err,error,*999)
    ELSE
      normalisedp=a/length
    ENDIF

    exits("NormaliseDP")
    RETURN
999 errorsexits("NormaliseDP",err,error)
    RETURN
    
  END FUNCTION normalisedp

  !
  !================================================================================================================================
  !

  SUBROUTINE normcrossproductsp(a,b,c,err,error,*)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: a(:)
    REAL(SP), INTENT(IN) :: b(:)
    REAL(SP), INTENT(OUT) :: c(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("NormCrossProductSP",err,error,*999)

    CALL crossproduct(a,b,c,err,error,*999)
    c=normalise(c,err,error)
    IF(err/=0) GOTO 999

    exits("NormCrossProductSP")
    RETURN
999 errorsexits("NormCrossProductSP",err,error)
    RETURN 1
    
  END SUBROUTINE normcrossproductsp
  
  !
  !================================================================================================================================
  !

  SUBROUTINE normcrossproductdp(a,b,c,err,error,*)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: a(:)
    REAL(DP), INTENT(IN) :: b(:)
    REAL(DP), INTENT(OUT) :: c(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    
    enters("NormCrossProductDP",err,error,*999)

    CALL crossproduct(a,b,c,err,error,*999)
    c=normalise(c,err,error)
    IF(err/=0) GOTO 999

    exits("NormCrossProductDP")
    RETURN
999 errorsexits("NormCrossProductDP",err,error)
    RETURN 1
    
  END SUBROUTINE normcrossproductdp
  
  !
  !================================================================================================================================
  !

  SUBROUTINE solvesmalllinearsystemsp(A,x,b,err,error,*)
  
    !Argument variables
    REAL(SP), INTENT(IN) :: A(:,:)
    REAL(SP), INTENT(OUT) :: x(:)
    REAL(SP), INTENT(IN) :: b(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    REAL(SP) :: AInv(size(a,1),size(a,2)),det
    
    enters("SolveSmallLinearSystemSP",err,error,*999)

    IF(SIZE(a,1)==SIZE(a,2)) THEN
      IF(SIZE(a,1)==SIZE(b,1)) THEN
        IF(SIZE(a,1)<=SIZE(x,1)) THEN
          SELECT CASE(SIZE(a,1))
          CASE(1:3)
            CALL invert(a,ainv,det,err,error,*999)
            x=matmul(ainv,b)
          CASE DEFAULT
            CALL flagerror("Matrix size not implemented.",err,error,*999)
          END SELECT
        ELSE
          CALL flagerror("x is too small.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Size of b is not the same as the number of rows in A.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not square.",err,error,*999)
    ENDIF

    exits("SolveSmallLinearSystemSP")
    RETURN
999 errorsexits("SolveSmallLinearSystemSP",err,error)
    RETURN 1
    
  END SUBROUTINE solvesmalllinearsystemsp

  !
  !================================================================================================================================
  !

  SUBROUTINE solvesmalllinearsystemdp(A,x,b,err,error,*)
  
    !Argument variables
    REAL(DP), INTENT(IN) :: A(:,:)
    REAL(DP), INTENT(OUT) :: x(:)
    REAL(DP), INTENT(IN) :: b(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    REAL(DP) :: AInv(size(a,1),size(a,2)),det
    
    enters("SolveSmallLinearSystemDP",err,error,*999)

    IF(SIZE(a,1)==SIZE(a,2)) THEN
      IF(SIZE(a,1)==SIZE(b,1)) THEN
        IF(SIZE(a,1)<=SIZE(x,1)) THEN
          SELECT CASE(SIZE(a,1))
          CASE(1:3)
            CALL invert(a,ainv,det,err,error,*999)
            x=matmul(ainv,b)
          CASE DEFAULT
            CALL flagerror("Matrix size not implemented.",err,error,*999)
          END SELECT
        ELSE
          CALL flagerror("x is too small.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Size of b is not the same as the number of rows in A.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not square.",err,error,*999)
    ENDIF

    exits("SolveSmallLinearSystemDP")
    RETURN
999 errorsexits("SolveSmallLinearSystemDP",err,error)
    RETURN 1
    
  END SUBROUTINE solvesmalllinearsystemdp

  !
  !================================================================================================================================
  !
  
  FUNCTION cothsp(a)

    !Argument variables
    REAL(SP), INTENT(IN) :: a
    !Function variable
    REAL(SP) :: CothSP
    
    cothsp=(exp(a)+exp(-1.0_sp*a))/(exp(a)-exp(-1.0_sp*a))

    RETURN
    
  END FUNCTION cothsp

  !
  !================================================================================================================================
  !

  FUNCTION cothdp(a)

    !Argument variables
    REAL(DP), INTENT(IN) :: a
    !Function variable
    REAL(DP) :: CothDP

    cothdp=(exp(a)+exp(-1.0_dp*a))/(exp(a)-exp(-1.0_dp*a))

    RETURN
    
  END FUNCTION cothdp

  !
  !================================================================================================================================
  !

  SUBROUTINE spline_cubic_set (n, t, y, ibcbeg, ybcbeg, ibcend, ybcend, ypp, err, error, *)

    !Argument variables
    INTEGER(INTG), INTENT(IN) :: n
    REAL(DP), INTENT(IN) :: t(n)
    REAL(DP), INTENT(IN) :: y(n)
    INTEGER(INTG), INTENT(IN) :: ibcbeg
    REAL(DP), INTENT(IN) :: ybcbeg
    INTEGER(INTG), INTENT(IN) :: ibcend
    REAL(DP), INTENT(IN) :: ybcend
    REAL(DP), INTENT(OUT) :: ypp(n)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    REAL(DP) :: diag(n)
    REAL(DP) :: sub(2:n)
    REAL(DP) :: sup(1:n-1)
    INTEGER(INTG) :: i
    TYPE(varying_string) :: localError

    enters("spline_cubic_set",err,error,*999)

    ! Sanity checks
    IF ( n <= 1 ) then
      localerror="spline interpolation requires at least 2 knots- user supplied "//trim(numbertovstring(n,"*",err,error))
      CALL flagerror(localerror,err,error,*999)
    ENDIF
    DO i = 1, n-1
      IF ( t(i) >= t(i+1) ) then
        localerror="Non-increasing knots supplied for cubic spline interpolation."
        CALL flagerror(localerror,err,error,*999)
      ENDIF
    ENDDO

    !  Set the first equation.
    IF ( ibcbeg == 0 ) then
      ypp(1) = 0.0e+00_dp
      diag(1) = 1.0e+00_dp
      sup(1) = -1.0e+00_dp
    ELSE IF ( ibcbeg == 1 ) then
      ypp(1) = ( y(2) - y(1) ) / ( t(2) - t(1) ) - ybcbeg
      diag(1) = ( t(2) - t(1) ) / 3.0e+00_dp
      sup(1) = ( t(2) - t(1) ) / 6.0e+00_dp
    ELSE IF ( ibcbeg == 2 ) then
      ypp(1) = ybcbeg
      diag(1) = 1.0e+00_dp
      sup(1) = 0.0e+00_dp
    ELSE
      localerror="The boundary flag IBCBEG must be 0, 1 or 2."
      CALL flagerror(localerror,err,error,*999)
    ENDIF

    !  Set the intermediate equations.
    DO i = 2, n-1
      ypp(i) = ( y(i+1) - y(i) ) / ( t(i+1) - t(i) ) &
       &     - ( y(i) - y(i-1) ) / ( t(i) - t(i-1) )
      sub(i) = ( t(i) - t(i-1) ) / 6.0e+00_dp
      diag(i) = ( t(i+1) - t(i-1) ) / 3.0e+00_dp
      sup(i) = ( t(i+1) - t(i) ) / 6.0e+00_dp
    ENDDO

    !  Set the last equation.
    IF ( ibcend == 0 ) then
      ypp(n) = 0.0e+00_dp
      sub(n) = -1.0e+00_dp
      diag(n) = 1.0e+00_dp
    ELSE IF ( ibcend == 1 ) then
      ypp(n) = ybcend - ( y(n) - y(n-1) ) / ( t(n) - t(n-1) )
      sub(n) = ( t(n) - t(n-1) ) / 6.0e+00_dp
      diag(n) = ( t(n) - t(n-1) ) / 3.0e+00_dp
    ELSE IF ( ibcend == 2 ) then
      ypp(n) = ybcend
      sub(n) = 0.0e+00_dp
      diag(n) = 1.0e+00_dp
    ELSE
      localerror="The boundary flag IBCEND must be 0, 1 or 2."
      CALL flagerror(localerror,err,error,*999)
    ENDIF

    !  Special case:
    !    N = 2, IBCBEG = IBCEND = 0.
    IF ( n == 2 .and. ibcbeg == 0 .and. ibcend == 0 ) then
      ypp(1) = 0.0e+00_dp
      ypp(2) = 0.0e+00_dp

    !  Solve the linear system.
    ELSE
      CALL s3_fs ( sub, diag, sup, n, ypp, ypp, err, error, *999 )
    ENDIF

    exits("spline_cubic_set")
    RETURN
999 errorsexits("spline_cubic_set",err,error)
    RETURN 1
    
  END SUBROUTINE spline_cubic_set

  !
  !================================================================================================================================
  !

  SUBROUTINE s3_fs ( a1, a2, a3, n, b, x, err, error, *)

    !Argument variables
    REAL(DP), INTENT(INOUT) :: a1(2:n)
    REAL(DP), INTENT(INOUT) :: a2(1:n)
    REAL(DP), INTENT(INOUT) :: a3(1:n-1)
    INTEGER(INTG), INTENT(IN) :: n
    REAL(DP), INTENT(INOUT) :: b(n)
    REAL(DP), INTENT(OUT) :: x(n)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    INTEGER(INTG) :: i
    REAL(DP) :: xmult
    TYPE(varying_string) :: localError

    enters("s3_fs",err,error,*999)

    !  The diagonal entries can't be zero.
    DO i = 1, n
      IF ( abs(a2(i)) < zero_tolerance ) then
        localerror="Zero diagonal entry in tridiagonal linear system."
        CALL flagerror(localerror,err,error,*999)
      ENDIF
    ENDDO

    DO i = 2, n-1
      xmult = a1(i) / a2(i-1)
      a2(i) = a2(i) - xmult * a3(i-1)
      b(i) = b(i) - xmult * b(i-1)
    ENDDO

    xmult = a1(n) / a2(n-1)
    a2(n) = a2(n) - xmult * a3(n-1)
    x(n) = ( b(n) - xmult * b(n-1) ) / a2(n)
    DO i = n-1, 1, -1
      x(i) = ( b(i) - a3(i) * x(i+1) ) / a2(i)
    ENDDO

    exits("s3_fs")
    RETURN
999 errorsexits("s3_fs",err,error)
    RETURN 1
    
  END SUBROUTINE s3_fs

  !
  !================================================================================================================================
  !
  
  SUBROUTINE spline_cubic_val (n, t, y, ypp, tval, yval, ypval, yppval, err, error, *)

    !Argument variables
    INTEGER(INTG), INTENT(IN) :: n
    REAL(DP), INTENT(IN) :: t(n)
    REAL(DP), INTENT(IN) :: y(n)
    REAL(DP), INTENT(IN) :: ypp(n)
    REAL(DP), INTENT(IN) :: tval
    REAL(DP), INTENT(OUT) :: yval
    REAL(DP), INTENT(OUT) :: ypval
    REAL(DP), INTENT(OUT) :: yppval
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    REAL(DP) :: dt
    REAL(DP) :: h
    INTEGER(INTG) :: i
    INTEGER(INTG) :: left
    INTEGER(INTG) :: right
    LOGICAL :: foundInterval

    enters("spline_cubic_val",err,error,*999)

    !  Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL.
    !  Values below T(1) or above T(N) use extrapolation.
    foundinterval = .false.
    DO i = 2, n - 1
      IF ( tval < t(i) ) THEN
        foundinterval=.true.
        left = i - 1
        right = i
        EXIT
      ENDIF
    ENDDO
    IF (foundinterval .EQV. .false.) THEN
      left = n - 1
      right = n
    ENDIF

    !  Evaluate the polynomial.
    dt = tval - t(left)
    h = t(right) - t(left)

    yval = y(left) &
     &   + dt * ( ( y(right) - y(left) ) / h &
     &          - ( ypp(right) / 6.0e+00_dp + ypp(left) / 3.0e+00_dp ) * h &
     &   + dt * ( 0.5e+00 * ypp(left) &
     &   + dt * ( ( ypp(right) - ypp(left) ) / ( 6.0e+00_dp * h ) ) ) )

    ypval = ( y(right) - y(left) ) / h &
     &   - ( ypp(right) / 6.0e+00_dp + ypp(left) / 3.0e+00_dp ) * h &
     &   + dt * ( ypp(left) &
     &   + dt * ( 0.5e+00_dp * ( ypp(right) - ypp(left) ) / h ) )

    yppval = ypp(left) + dt * ( ypp(right) - ypp(left) ) / h 

    exits("spline_cubic_val")
    RETURN
999 errorsexits("spline_cubic_val",err,error)
    RETURN 1
    
  END SUBROUTINE spline_cubic_val

  !
  !================================================================================================================================
  !

END MODULE maths