libMath Module

Procedures to manipulate matrices and vectors


Uses


Variables

Type Visibility Attributes Name Initial
integer, public, parameter :: dp = kind(1.d0)

Double precision setting
real(kind=dp), public, parameter :: pi = atan(1._dp)*4._dp

Value of pi
real(kind=dp), public, parameter :: eps = epsilon(1._dp)

Value of machine epsilon
real(kind=dp), public, parameter :: degToRad = pi/180._dp

For converting degree to radian
real(kind=dp), public, parameter :: radToDeg = 180._dp/pi

For converting radian to degree
real(kind=dp), public, parameter, dimension(3) :: xAxis = [1._dp, 0._dp, 0._dp]

Global x-axis
real(kind=dp), public, parameter, dimension(3) :: yAxis = [0._dp, 1._dp, 0._dp]

Global y-axis
real(kind=dp), public, parameter, dimension(3) :: zAxis = [0._dp, 0._dp, 1._dp]

Global z-axis

Interfaces

public interface lsq2

  • public function lsq2_scalar(xQuery, xData, yData)

    Linear Least Squares fitting (2nd order)

    Arguments

    Type IntentOptional Attributes Name
    real(kind=dp), intent(in) :: xQuery
    real(kind=dp), intent(in), dimension(:) :: xData
    real(kind=dp), intent(in), dimension(:) :: yData

    Return Value real(kind=dp)

  • public function lsq2_array(xQuery, xData, yData)

    Arguments

    Type IntentOptional Attributes Name
    real(kind=dp), intent(in), dimension(:) :: xQuery
    real(kind=dp), intent(in), dimension(:) :: xData
    real(kind=dp), intent(in), dimension(:) :: yData

    Return Value real(kind=dp), dimension(size(xQuery))


Functions

public function isFloatEqual(a, b, tol)

Checks if a == b within a tolerance

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: a

Real number a
real(kind=dp), intent(in) :: b

Real number b
real(kind=dp), optional :: tol

Tolerance is epsilon if not specified

Return Value logical

public function inv2(A) result(Ainv)

Inverse of a matrix calculated by finding the LU

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(:, :) :: A

Square matrix A

Return Value real(kind=dp), dimension(size(A,1),size(A,2))

Inverse of matrix A

public function matmul2(A, B) result(AB)

Matrix multiplication

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(:, :) :: A

Matrix A
real(kind=dp), intent(in), dimension(:, :) :: B

Matrix B

Return Value real(kind=dp), dimension(size(A, 1), size(B, 2))

public function matmulAX(A, X) result(AX)

Matrix multiplication with vector

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(:, :) :: A

Matrix A
real(kind=dp), intent(in), dimension(:) :: X

Vector X

Return Value real(kind=dp), dimension(size(A, 1))

Product of matrix A and vector X

public function length3d(P1, P2) result(length)

Compute length of linesegment between points P1 and P2

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(3) :: P1

Coordinates of point P1
real(kind=dp), intent(in), dimension(3) :: P2

Coordinates of point P2

Return Value real(kind=dp)

Length of linesegment

public function linspace(xstart, xend, nx) result(xout)

Return linearly-spaced array over a specified interval

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: xstart

Start value of sequence
real(kind=dp), intent(in) :: xend

End value of sequence
integer, intent(in) :: nx

Number of points sequence

Return Value real(kind=dp), dimension(nx)

Array of real numbers

public function cosspace(xstart, xend, nx) result(xout)

Return cossine-spaced array over a specified interval

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: xstart
real(kind=dp), intent(in) :: xend
integer, intent(in) :: nx

Return Value real(kind=dp), dimension(nx)

public function halfsinspace(xstart, xend, nx) result(xout)

Return half sine-spaced array over a specified interval

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: xstart
real(kind=dp), intent(in) :: xend
integer, intent(in) :: nx

Return Value real(kind=dp), dimension(nx)

public function tanspace(xstart, xend, nx) result(xout)

Return tan-spaced array over a specified interval

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: xstart
real(kind=dp), intent(in) :: xend
integer, intent(in) :: nx

Return Value real(kind=dp), dimension(nx)

public pure function cross_product(aVec, bVec)

Compute cross product between two 3d vectors

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(3) :: aVec
real(kind=dp), intent(in), dimension(3) :: bVec

Return Value real(kind=dp), dimension(3)

public pure function outer_product(aVec, bVec)

Compute outer product between two 3d vectors

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(3) :: aVec
real(kind=dp), intent(in), dimension(3) :: bVec

Return Value real(kind=dp), dimension(3, 3)

public pure function getAngleTan(aVec, bVec)

Angle between two 3d vectors using tan formula Result will be from -pi to pi

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(3) :: aVec
real(kind=dp), intent(in), dimension(3) :: bVec

Return Value real(kind=dp)

public pure function getAngleCos(aVec, bVec)

Angle between two 3d vectors using cos formula Assumes no angle > 180 deg exists

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(3) :: aVec
real(kind=dp), intent(in), dimension(3) :: bVec

Return Value real(kind=dp)

public pure function unitVec(aVec)

Normalizes a non-zero 3d vector

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(3) :: aVec

Return Value real(kind=dp), dimension(3)

public pure function projVec(aVec, dirVec)

Returns 3d vector aVec projected along dirVec

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(3) :: aVec
real(kind=dp), intent(in), dimension(3) :: dirVec

Return Value real(kind=dp), dimension(3)

public pure function noProjVec(aVec, dirVec)

Removes component along dirVec from a vector aVec

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(3) :: aVec
real(kind=dp), intent(in), dimension(3) :: dirVec

Return Value real(kind=dp), dimension(3)

public function inv(A) result(Ainv)

Native implementation of matrix inverse based on Doolittle method

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(:, :) :: A

Return Value real(kind=dp), dimension(size(A, 1), size(A, 1))

public function isInverse(A, Ainv)

Check if Ainv is inverse of matrix A by multiplication

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(:, :) :: A
real(kind=dp), intent(in), dimension(:, :) :: Ainv

Return Value logical

public function pwl_interp1d(x, y, q)

Piecewise linear 1d interpolation

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(:) :: x
real(kind=dp), intent(in), dimension(:) :: y
real(kind=dp), intent(in) :: q

Return Value real(kind=dp)

public function interp1d(xq, x, y, order)

1-d Interpolation using 1st and 2nd order Lagrange polynomials

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: xq
real(kind=dp), intent(in), dimension(:) :: x
real(kind=dp), intent(in), dimension(:) :: y
integer, intent(in) :: order

Return Value real(kind=dp)

public function lsq2_scalar(xQuery, xData, yData)

Linear Least Squares fitting (2nd order)

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: xQuery
real(kind=dp), intent(in), dimension(:) :: xData
real(kind=dp), intent(in), dimension(:) :: yData

Return Value real(kind=dp)

public function lsq2_array(xQuery, xData, yData)

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(:) :: xQuery
real(kind=dp), intent(in), dimension(:) :: xData
real(kind=dp), intent(in), dimension(:) :: yData

Return Value real(kind=dp), dimension(size(xQuery))

public function Tbg(cs_phi, cs_theta, cs_psi)

Transformation matrix (body to global frame)

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(2) :: cs_phi

cos(phi), sin(phi)
real(kind=dp), intent(in), dimension(2) :: cs_theta

cos(theta), sin(theta)
real(kind=dp), intent(in), dimension(2) :: cs_psi

cos(psi), sin(psi)

Return Value real(kind=dp), dimension(3, 3)

3x3 transformation matrix

public function Tgb(cs_phi, cs_theta, cs_psi)

Transformation matrix (global to body frame)

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(2) :: cs_phi

cos(phi), sin(phi)
real(kind=dp), intent(in), dimension(2) :: cs_theta

cos(theta), sin(theta)
real(kind=dp), intent(in), dimension(2) :: cs_psi

cos(psi), sin(psi)

Return Value real(kind=dp), dimension(3, 3)

3x3 transformation matrix

public function getTransformAxis(theta, axisVec) result(TMat)

Transformation matrix for theta angular rotation about a 3d axis

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: theta

theta angle in radians
real(kind=dp), intent(in), dimension(3) :: axisVec

3d axis vector

Return Value real(kind=dp), dimension(3, 3)

3x3 transformation matrix

public function trapz(y, x)

Trapezoid integration with unequal intervals

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(:) :: y
real(kind=dp), intent(in), dimension(:) :: x

Return Value real(kind=dp)


Subroutines

public subroutine print_mat(M)

Display in matrix format

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(:, :) :: M