maths_library.f90 Source File


Contents

Source Code


Source Code

module maths_library

  !! Library of mathematical and numerical routines
  !! author: P J Knight, CCFE, Culham Science Centre
  !! N/A
  !! This module contains a large number of routines to enable
  !! PROCESS to perform a variety of numerical procedures, including
  !! linear algebra, zero finding, integration and minimisation.
  !! The module is an amalgamation of the contents of several
  !! different pre-existing PROCESS source files, which themselves
  !! were derived from a number of different numerical libraries
  !! including BLAS and MINPAC.
  !! http://en.wikipedia.org/wiki/Gamma_function
  !
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#ifndef dp
  use, intrinsic :: iso_fortran_env, only: dp=>real64
#endif

  !use process_output

  implicit none

contains

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function find_y_nonuniform_x(x0,x,y,n)

    !! Routine to find y0 such that y0 = y(x0) given a set of
    !! values x(1:n), y(1:n)
    !! author: P J Knight, CCFE, Culham Science Centre
    !! x0 : input real : x value at which we want to find y
    !! x(1:n) : input real array : monotonically de/increasing x values
    !! y(1:n) : input real array : y values at each x
    !! n : input integer : size of array
    !! This routine performs a simple linear interpolation method
    !! to find the y value at x = x0. If x0 lies outside the
    !! range [x(1),x(n)], the y value at the nearest 'end' of the data
    !! is returned.
    !! None
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    real(dp) :: find_y_nonuniform_x

    !  Arguments

    integer, intent(in) :: n
    real(dp), intent(in) :: x0
    real(dp), dimension(n), intent(in) :: x
    real(dp), dimension(n), intent(in) :: y

    !  Local variables
    integer :: i,j

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Step through arrays until x crosses the value of interest

    j = 0
    rough_search: do i = 1,n-1

       if (((x(i)-x0)*(x(i+1)-x0)) <= 0.0D0) then
          j = i
          exit rough_search
       end if

    end do rough_search

    if (j /= 0) then

       !  Simply do a linear interpolation between the two grid points
       !  spanning the point of interest

       find_y_nonuniform_x = y(j) + (y(j+1)-y(j))*(x0-x(j))/(x(j+1)-x(j))

    else  !  No points found, so return the 'nearest' y value

       if (x(n) > x(1)) then  !  values are monotonically increasing
          if (x0 > x(n)) then
             find_y_nonuniform_x = y(n)
          else
             find_y_nonuniform_x = y(1)
          end if

       else  !  values are monotonically decreasing
          if (x0 < x(n)) then
             find_y_nonuniform_x = y(n)
          else
             find_y_nonuniform_x = y(1)
          end if

       end if

    end if

  end function find_y_nonuniform_x

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ellipke(sqk,kk,ek)

    !! Routine that calculates the complete elliptic integral
    !! of the first and second kinds
    !! author: P J Knight, CCFE, Culham Science Centre
    !! sqk : input real : square of the elliptic modulus
    !! kk  : output real : complete elliptic integral of the first kind
    !! ek  : output real : complete elliptic integral of the second kind
    !! This routine calculates the complete elliptic integral
    !! of the first and second kinds.
    !! <P>The method used is that described in the reference, and
    !! the code is taken from the Culham maglib library routine FN02A.
    !! Approximations for Digital Computers, C. Hastings,
    !! Princeton University Press, 1955
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    !  Arguments

    real(dp), intent(in) :: sqk
    real(dp), intent(out) :: kk,ek

    !  Local variables

    real(dp) :: a,b,d,e

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    d = 1.0D0 - sqk
    e = log(d)

    !  Evaluate series for integral of first kind

    a = (((0.014511962D0*d + 0.037425637D0)*d + 0.035900924D0)*d &
         + 0.096663443D0)*d + 1.386294361D0
    b = (((0.004417870D0*d + 0.033283553D0)*d + 0.06880249D0)*d &
         + 0.12498594D0)*d + 0.5D0

    kk = a - b*e

    !  Evaluate series for integral of second kind

    a = (((0.017365065D0*d + 0.047573835D0)*d + 0.06260601D0)*d &
         + 0.44325141D0)*d + 1.0D0
    b = (((0.005264496D0*d + 0.040696975D0)*d + 0.09200180D0)*d &
         + 0.24998368D0)*d

    ek = a - b*e

  end subroutine ellipke

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  function binomial(n,k) result(coefficient)
    ! This outputs a real approximation to the coefficient
    ! http://en.wikipedia.org/wiki/Binomial_coefficient#Multiplicative_formula
    implicit none
    integer, intent(in) :: n, k
    real(dp) :: coefficient
    integer :: numerator, i
    if (k == 0) then
        coefficient = 1
    else
        coefficient = 1.0D0
        do i = 1, k
            numerator = n + 1 -i
            coefficient = coefficient * real(numerator)/real(i)
        end do
    end if
  end function binomial

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1)

    !! Singular Value Decomposition
    !! author: P J Knight, CCFE, Culham Science Centre
    !! author: B. S. Garbow, Applied Mathematics Division, Argonne National Laboratory
    !! nm : input integer : Max number of rows of arrays a, u, v; >= m,n
    !! m : input integer : Actual number of rows of arrays a, u
    !! n : input integer : Number of columns of arrays a, u, and the order of v
    !! a(nm,n) : input/output real array : On input matrix to be decomposed;
    !! on output, either unchanged or overwritten with u or v
    !! w(n) : output real array : The n (non-negative) singular values of a
    !! (the diagonal elements of s); unordered.  If an error exit
    !! is made, the singular values should be correct for indices
    !! ierr+1,ierr+2,...,n.
    !! matu : input logical : Set to .true. if the u matrix in the
    !! decomposition is desired, and to .false. otherwise.
    !! u(nm,n) : output real array : The matrix u (orthogonal column vectors)
    !! of the decomposition if matu has been set to .true., otherwise
    !! u is used as a temporary array.  u may coincide with a.
    !! If an error exit is made, the columns of u corresponding
    !! to indices of correct singular values should be correct.
    !! matv : input logical : Set to .true. if the v matrix in the
    !! decomposition is desired, and to .false. otherwise.
    !! v(nm,n) : output real array : The matrix v (orthogonal) of the
    !! decomposition if matv has been set to .true., otherwise
    !! v is not referenced.  v may also coincide with a if u is
    !! not needed.  If an error exit is made, the columns of v
    !! corresponding to indices of correct singular values
    !! should be correct.
    !! ierr : output integer :  zero for normal return, or <I>k</I> if the
    !! k-th singular value has not been determined after 30 iterations.
    !! rv1(n) : output real array : work array
    !! This subroutine is a translation of the algol procedure SVD,
    !! Num. Math. 14, 403-420(1970) by Golub and Reinsch,
    !! Handbook for Auto. Comp., vol II - Linear Algebra, 134-151(1971).
    !! <P>It determines the singular value decomposition
    !! <I>a=usv<SUP>t</SUP></I> of a real m by n rectangular matrix.
    !! Householder bidiagonalization and a variant of the QR
    !! algorithm are used.
    !! None
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    !  Arguments

    integer, intent(in) :: nm, m, n
    logical, intent(in) :: matu, matv
    real(dp), dimension(nm,n), intent(inout) :: a
    real(dp), dimension(nm,n), intent(out) :: u, v
    real(dp), dimension(n), intent(out) :: w, rv1
    integer, intent(out) :: ierr

    !  Local variables

    integer :: i,j,k,l,ii,i1,kk,k1,ll,l1,mn,its
    real(dp) :: c,f,g,h,s,x,y,z,scale,anorm

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ierr = 0

    u = a

    !  Householder reduction to bidiagonal form

    g = 0.0D0
    scale = 0.0D0
    anorm = 0.0D0

    do i = 1, n

       l = i + 1
       rv1(i) = scale * g
       g = 0.0D0
       s = 0.0D0
       scale = 0.0D0

       if (i <= m) then

          do k = i, m
             scale = scale + abs(u(k,i))
          end do

          if (scale /= 0.0D0) then

             do k = i, m
                u(k,i) = u(k,i) / scale
                s = s + u(k,i)**2
             end do

             f = u(i,i)
             g = -sign(sqrt(s),f)
             h = f * g - s
             u(i,i) = f - g

             if (i /= n) then
                do j = l, n
                   s = 0.0D0
                   do k = i, m
                      s = s + u(k,i) * u(k,j)
                   end do
                   f = s / h
                   do k = i, m
                      u(k,j) = u(k,j) + f * u(k,i)
                   end do
                end do
             end if

             do k = i, m
                u(k,i) = scale * u(k,i)
             end do

          end if

       end if

       w(i) = scale * g
       g = 0.0D0
       s = 0.0D0
       scale = 0.0D0

       if (.not.((i > m) .or. (i == n))) then

          do k = l, n
             scale = scale + abs(u(i,k))
          end do

          if (scale /= 0.0D0) then

             do k = l, n
                u(i,k) = u(i,k) / scale
                s = s + u(i,k)**2
             end do

             f = u(i,l)
             g = -sign(sqrt(s),f)
             h = f * g - s
             u(i,l) = f - g

             do k = l, n
                rv1(k) = u(i,k) / h
             end do

             if (i /= m) then
                do j = l, m
                   s = 0.0D0
                   do k = l, n
                      s = s + u(j,k) * u(i,k)
                   end do
                   do k = l, n
                      u(j,k) = u(j,k) + s * rv1(k)
                   end do
                end do
             end if

             do k = l, n
                u(i,k) = scale * u(i,k)
             end do

          end if

       end if

       anorm = max(anorm,abs(w(i))+abs(rv1(i)))

    end do  ! i

    !  Accumulation of right-hand transformations

    if (matv) then

       !  For i=n step -1 until 1 do
       do ii = 1, n
          i = n + 1 - ii
          if (i /= n) then

             if (g /= 0.0D0) then
                do j = l, n
                   !  Double division avoids possible underflow
                   v(j,i) = (u(i,j) / u(i,l)) / g
                end do
                do j = l, n
                   s = 0.0D0
                   do k = l, n
                      s = s + u(i,k) * v(k,j)
                   end do
                   do k = l, n
                      v(k,j) = v(k,j) + s * v(k,i)
                   end do
                end do
             end if

             do j = l, n
                v(i,j) = 0.0D0
                v(j,i) = 0.0D0
             end do

          end if

          v(i,i) = 1.0D0
          g = rv1(i)
          l = i
       end do

    end if

    !  Accumulation of left-hand transformations

    if (matu) then

       !  For i=min(m,n) step -1 until 1 do
       mn = n
       if (m < n) mn = m

       do ii = 1, mn
          i = mn + 1 - ii
          l = i + 1
          g = w(i)
          if (i /= n) then
             do j = l, n
                u(i,j) = 0.0D0
             end do
          end if

          if (g /= 0.0D0) then

             if (i /= mn) then
                do j = l, n
                   s = 0.0D0
                   do k = l, m
                      s = s + u(k,i) * u(k,j)
                   end do
                   f = (s / u(i,i)) / g  !  Double division avoids possible underflow
                   do k = i, m
                      u(k,j) = u(k,j) + f * u(k,i)
                   end do
                end do
             end if

             do j = i, m
                u(j,i) = u(j,i) / g
             end do

          else
             do j = i, m
                u(j,i) = 0.0D0
             end do
          end if

          u(i,i) = u(i,i) + 1.0D0

       end do

    end if

    !  Diagonalization of the bidiagonal form
    !  For k=n step -1 until 1 do

    do kk = 1, n
       k1 = n - kk
       k = k1 + 1
       its = 0

       !  Test for splitting.
       !  For l=k step -1 until 1 do

       do
          do ll = 1, k
             l1 = k - ll
             l = l1 + 1
             if ((abs(rv1(l)) + anorm) == anorm) goto 470

             !  rv1(1) is always zero, so there is no exit
             !  through the bottom of the loop

             !+**PJK 23/05/06 Prevent problems from the code getting here with l1=0
             if (l1 == 0) then
                write(*,*) 'SVD: Shouldn''t get here...'
                goto 470
             end if

             if ((abs(w(l1)) + anorm) == anorm) exit
          end do

          !  Cancellation of rv1(l) if l greater than 1

          c = 0.0D0
          s = 1.0D0

          do i = l, k
             f = s * rv1(i)
             rv1(i) = c * rv1(i)
             if ((abs(f) + anorm) == anorm) exit
             g = w(i)
             h = sqrt(f*f+g*g)
             w(i) = h
             c = g / h
             s = -f / h
             if (.not. matu) cycle

             do j = 1, m
                y = u(j,l1)
                z = u(j,i)
                u(j,l1) = y * c + z * s
                u(j,i) = -y * s + z * c
             end do
          end do

470       continue

          !  Test for convergence

          z = w(k)
          if (l == k) exit

          !  Shift from bottom 2 by 2 minor

          if (its == 30) then
             !  Set error - no convergence to a
             !  singular value after 30 iterations
             ierr = k
             return
          end if

          its = its + 1
          x = w(l)
          y = w(k1)
          g = rv1(k1)
          h = rv1(k)
          f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.D0 * h * y)
          g = sqrt(f*f+1.D0)
          f = ((x - z) * (x + z) + h * (y / (f + sign(g,f)) - h)) / x

          !  Next QR transformation

          c = 1.0D0
          s = 1.0D0

          do i1 = l, k1
             i = i1 + 1
             g = rv1(i)
             y = w(i)
             h = s * g
             g = c * g
             z = sqrt(f*f+h*h)
             rv1(i1) = z
             c = f / z
             s = h / z
             f = x * c + g * s
             g = -x * s + g * c
             h = y * s
             y = y * c

             if (matv) then
                do j = 1, n
                   x = v(j,i1)
                   z = v(j,i)
                   v(j,i1) = x * c + z * s
                   v(j,i) = -x * s + z * c
                end do
             end if

             z = sqrt(f*f+h*h)
             w(i1) = z

             !  Rotation can be arbitrary if z is zero

             if (z /= 0.0D0) then
                c = f / z
                s = h / z
             end if

             f = c * g + s * y
             x = -s * g + c * y
             if (.not. matu) cycle

             do j = 1, m
                y = u(j,i1)
                z = u(j,i)
                u(j,i1) = y * c + z * s
                u(j,i) = -y * s + z * c
             end do

          end do

          rv1(l) = 0.0D0
          rv1(k) = f
          w(k) = x

       end do

       !  Convergence

       if (z >= 0.0D0) cycle

       !  w(k) is made non-negative
       w(k) = -z
       if (.not. matv) cycle

       do j = 1, n
          v(j,k) = -v(j,k)
       end do

    end do

  end subroutine svd

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine eshellvol(rshell,rmini,rmino,zminor,drin,drout,dz,vin,vout,vtot)

    !! Routine to calculate the inboard, outboard and total volumes
    !! of a toroidal shell comprising two elliptical sections
    !! author: P J Knight, CCFE, Culham Science Centre
    !! rshell : input real : major radius of centre of both ellipses (m)
    !! rmini  : input real : horizontal distance from rshell to outer edge
    !! of inboard elliptical shell (m)
    !! rmino  : input real : horizontal distance from rshell to inner edge
    !! of outboard elliptical shell (m)
    !! zminor : input real : vertical internal half-height of shell (m)
    !! drin   : input real : horiz. thickness of inboard shell at midplane (m)
    !! drout  : input real : horiz. thickness of outboard shell at midplane (m)
    !! dz     : input real : vertical thickness of shell at top/bottom (m)
    !! vin    : output real : volume of inboard section (m3)
    !! vout   : output real : volume of outboard section (m3)
    !! vtot   : output real : total volume of shell (m3)
    !! This routine calculates the volume of the inboard and outboard sections
    !! of a toroidal shell defined by two co-centred semi-ellipses.
    !! Each section's internal and external surfaces are in turn defined
    !! by two semi-ellipses. The volumes of each section are calculated as
    !! the difference in those of the volumes of revolution enclosed by their
    !! inner and outer surfaces.
    !! <P>See also <A HREF="eshellarea.html"><CODE>eshellarea</CODE></A>
    !! Internal CCFE note T&amp;M/PKNIGHT/PROCESS/009, P J Knight:
    !! Surface Area and Volume Calculations for Toroidal Shells
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    use constants, only: pi, twopi
    implicit none

    !  Arguments
    real(kind=dp), intent(in) :: rshell, rmini, rmino, zminor, drin, drout, dz
    real(kind=dp), intent(out) :: vin, vout, vtot

    !  Local variables
    real(kind=dp) :: a, b, elong, v1, v2

    !  Global shared variables
    !  Input: pi,twopi
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! #TODO - Review both equations containing dz and attempt to separate
    !         top and bottom of vacuum vessel thickness
    !         See issue #433 for explanation

    !  Inboard section

    !  Volume enclosed by outer (higher R) surface of elliptical section
    !  and the vertical straight line joining its ends
    a = rmini ; b = zminor ; elong = b/a
    v1 = twopi * elong * (0.5D0*pi*rshell*a*a - 2.0D0/3.0D0*a*a*a)

    !  Volume enclosed by inner (lower R) surface of elliptical section
    !  and the vertical straight line joining its ends
    a = rmini+drin ; b = zminor+dz ; elong = b/a
    v2 = twopi * elong * (0.5D0*pi*rshell*a*a - 2.0D0/3.0D0*a*a*a)

    !  Volume of inboard section of shell
    vin = v2 - v1

    !  Outboard section

    !  Volume enclosed by inner (lower R) surface of elliptical section
    !  and the vertical straight line joining its ends
    a = rmino ; b = zminor ; elong = b/a
    v1 = twopi * elong * (0.5D0*pi*rshell*a*a + 2.0D0/3.0D0*a*a*a)

    !  Volume enclosed by outer (higher R) surface of elliptical section
    !  and the vertical straight line joining its ends
    a = rmino+drout ; b = zminor+dz ; elong = b/a
    v2 = twopi * elong * (0.5D0*pi*rshell*a*a + 2.0D0/3.0D0*a*a*a)

    !  Volume of outboard section of shell
    vout = v2 - v1

    !  Total shell volume
    vtot = vin + vout

  end subroutine eshellvol

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine dshellvol(rmajor,rminor,zminor,drin,drout,dz,vin,vout,vtot)

    !! Routine to calculate the inboard, outboard and total volumes
    !! of a D-shaped toroidal shell
    !! author: P J Knight, CCFE, Culham Science Centre
    !! rmajor : input real : major radius to outer point of inboard
    !! straight section of shell (m)
    !! rminor : input real : horizontal internal width of shell (m)
    !! zminor : input real : vertical internal half-height of shell (m)
    !! drin   : input real : horiz. thickness of inboard shell at midplane (m)
    !! drout  : input real : horiz. thickness of outboard shell at midplane (m)
    !! dz     : input real : vertical thickness of shell at top/bottom (m)
    !! vin    : output real : volume of inboard straight section (m3)
    !! vout   : output real : volume of outboard curved section (m3)
    !! vtot   : output real : total volume of shell (m3)
    !! This routine calculates the volume of the inboard and outboard sections
    !! of a D-shaped toroidal shell defined by the above input parameters.
    !! The inboard section is assumed to be a cylinder of uniform thickness.
    !! The outboard section's internal and external surfaces are defined
    !! by two semi-ellipses, centred on the outer edge of the inboard section;
    !! its volume is calculated as the difference in those of the volumes of
    !! revolution enclosed by the two surfaces.
    !! <P>See also <A HREF="dshellarea.html"><CODE>dshellarea</CODE></A>
    !! Internal CCFE note T&amp;M/PKNIGHT/PROCESS/009, P J Knight:
    !! Surface Area and Volume Calculations for Toroidal Shells
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    use constants, only: pi, twopi
    implicit none

    !  Arguments
    real(kind=dp), intent(in) :: rmajor, rminor, zminor, drin, drout, dz
    real(kind=dp), intent(out) :: vin, vout, vtot

    !  Local variables
    real(kind=dp) :: a, b, elong, v1, v2

    !  Global shared variables
    !  Input: pi,twopi
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! #TODO - Review both equations containing dz and attempt to separate
    !         top and bottom of vacuum vessel thickness
    !         See issue #433 for explanation

    !  Volume of inboard cylindrical shell
    vin = 2.0D0*(zminor+dz) * pi*(rmajor**2 - (rmajor-drin)**2)

    !  Volume enclosed by inner surface of elliptical outboard section
    !  and the vertical straight line joining its ends
    a = rminor ; b = zminor ; elong = b/a
    v1 = twopi * elong * (0.5D0*pi*rmajor*a*a + 2.0D0/3.0D0*a*a*a)

    !  Volume enclosed by outer surface of elliptical outboard section
    !  and the vertical straight line joining its ends
    a = rminor+drout ; b = zminor+dz ; elong = b/a
    v2 = twopi * elong * (0.5D0*pi*rmajor*a*a + 2.0D0/3.0D0*a*a*a)

    !  Volume of elliptical outboard shell
    vout = v2 - v1

    !  Total shell volume
    vtot = vin + vout

  end subroutine dshellvol

! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine dshellarea(rmajor,rminor,zminor,ain,aout,atot)

    !! Routine to calculate the inboard, outboard and total surface areas
    !! of a D-shaped toroidal shell
    !! author: P J Knight, CCFE, Culham Science Centre
    !! rmajor : input real : major radius of inboard straight section (m)
    !! rminor : input real : horizontal width of shell (m)
    !! zminor : input real : vertical half-height of shell (m)
    !! ain    : output real : surface area of inboard straight section (m3)
    !! aout   : output real : surface area of outboard curved section (m3)
    !! atot   : output real : total surface area of shell (m3)
    !! This routine calculates the surface area of the inboard and outboard
    !! sections of a D-shaped toroidal shell defined by the above input
    !! parameters.
    !! The inboard section is assumed to be a cylinder.
    !! The outboard section is defined by a semi-ellipse, centred on the
    !! major radius of the inboard section.
    !! <P>See also <A HREF="dshellvol.html"><CODE>dshellvol</CODE></A>
    !! Internal CCFE note T&amp;M/PKNIGHT/PROCESS/009, P J Knight:
    !! Surface Area and Volume Calculations for Toroidal Shells
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    use constants, only: pi, twopi
    implicit none

    !  Arguments
    real(dp), intent(in) :: rmajor,rminor,zminor
    real(dp), intent(out) :: ain,aout,atot

    !  Local variables
    real(dp) :: elong

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Area of inboard cylindrical shell
    ain = 4.0D0*zminor*pi*rmajor

    !  Area of elliptical outboard section
    elong = zminor/rminor
    aout = twopi * elong * (pi*rmajor*rminor + 2.0D0*rminor*rminor)

    !  Total surface area
    atot = ain + aout

  end subroutine dshellarea

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine eshellarea(rshell,rmini,rmino,zminor,ain,aout,atot)

    !! Routine to calculate the inboard, outboard and total surface areas
    !! of a toroidal shell comprising two elliptical sections
    !! author: P J Knight, CCFE, Culham Science Centre
    !! rshell : input real : major radius of centre of both ellipses (m)
    !! rmini  : input real : horizontal distance from rshell to
    !! inboard elliptical shell (m)
    !! rmino  : input real : horizontal distance from rshell to
    !! outboard elliptical shell (m)
    !! zminor : input real : vertical internal half-height of shell (m)
    !! ain    : output real : surface area of inboard section (m3)
    !! aout   : output real : surface area of outboard section (m3)
    !! atot   : output real : total surface area of shell (m3)
    !! This routine calculates the surface area of the inboard and outboard
    !! sections of a toroidal shell defined by two co-centred semi-ellipses.
    !! <P>See also <A HREF="eshellvol.html"><CODE>eshellvol</CODE></A>
    !! Internal CCFE note T&amp;M/PKNIGHT/PROCESS/009, P J Knight:
    !! Surface Area and Volume Calculations for Toroidal Shells
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    use constants, only: pi, twopi
    implicit none

    !  Arguments
    real(dp), intent(in) :: rshell,rmini,rmino,zminor
    real(dp), intent(out) :: ain,aout,atot

    !  Local variables
    real(dp) :: elong

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Inboard section
    elong = zminor/rmini
    ain = twopi * elong * (pi*rshell*rmini - 2.0D0*rmini*rmini)

    !  Outboard section
    elong = zminor/rmino
    aout = twopi * elong * (pi*rshell*rmino + 2.0D0*rmino*rmino)

    !  Total surface area
    atot = ain + aout

  end subroutine eshellarea

  ! ------------------------------------------------------------------------
  pure function variable_error(variable)
      real(dp), intent(in) ::variable
      logical::variable_error

      if((variable/=variable).or.(variable<-9.99D99).or.(variable>9.99D99))then
          variable_error = .TRUE.
      else
          variable_error = .FALSE.
      end if

  end function variable_error

  ! ------------------------------------------------------------------------
  pure function integer2string(value)
      ! Convert an integer value to a 2-digit string with leading zero if required.
      integer, intent(in) :: value
      character(len=2) integer2string
      write (integer2string,'(I2.2)') value
  end function integer2string

  pure function integer3string(value)
      ! Convert an integer value to a 3-digit string with leading zero if required.
      integer, intent(in) :: value
      character(len=3) integer3string
      write (integer3string,'(I3.3)') value
  end function integer3string

end module maths_library