fun

package
v1.2.12 Latest Latest
Warning

This package is not in the latest version of its module.

 Go to latest Published: May 25, 2023 License: BSD-3-Clause Imports: 8 Imported by: 39

README

Gosl. fun. Special functions, DFT, FFT, Bessel, elliptical integrals, orthogonal polynomials, interpolators

Go Reference

This package implements special functions such as orthogonal polynomials and elliptical functions of first, second and third kind.

Routines to interpolate and/or assist on spectral methods are also available; e.g. FourierInterp, ChebyInterp.

API

Please see the documentation here

Documentation

Overview

Package fun (functions) implements special functions such as elliptical, orthogonal polynomials, Bessel, discrete Fourier transform, polynomial interpolators, and more.

Index

Constants

This section is empty.

Variables

This section is empty.

Functions

func Atan2p 

func Atan2p(y, x float64) (αrad float64)

Atan2p implements a positive version of atan2, in such a way that: 0 ≤ α ≤ 2π

func Atan2pDeg 

func Atan2pDeg(y, x float64) (αdeg float64)

Atan2pDeg implements a positive version of atan2, in such a way that: 0 ≤ α ≤ 360

func Beta 

func Beta(a, b float64) float64

Beta computes the beta function by calling the Lgamma function

func Binomial 

func Binomial(n, k int) float64

Binomial comptues the binomial coefficient (n k)^T

func Boxcar 

func Boxcar(x, a, b float64) float64

Boxcar implements the boxcar function 

Boxcar(x;a,b) = Heav(x-a) - Heav(x-b)

                │ 0    if x < a or  x > b
Boxcar(x;a,b) = ┤ 1/2  if x = a or  x = b
                │ 1    if x > a and x < b

Note: a ≤ x ≤ b; i.e. b ≥ a (not checked)

func CarlsonRc 

func CarlsonRc(x, y float64) float64

CarlsonRc computes Carlson’s degenerate elliptic integral according to [1] Computes Rc(x,y) where x must be nonnegative and y must be nonzero. If y < 0, the Cauchy principal value is returned.

func CarlsonRd 

func CarlsonRd(x, y, z float64) float64

CarlsonRd computes Carlson’s elliptic integral of the second kind according to [1] Computes Rf(x,y,z) where x,y must be non-negative and at most one can be zero. z must be positive. 

References:
[1] Press WH, Teukolsky SA, Vetterling WT, Fnannery BP (2007) Numerical Recipes: The Art of
    Scientific Computing. Third Edition. Cambridge University Press. 1235p.

func CarlsonRf 

func CarlsonRf(x, y, z float64) float64

CarlsonRf computes Carlson's elliptic integral of the first kind according to [1]. See also [2] Computes Rf(x,y,z) where x,y,z must be non-negative and at most one can be zero. 

References:
[1] Press WH, Teukolsky SA, Vetterling WT, Fnannery BP (2007) Numerical Recipes: The Art of
    Scientific Computing. Third Edition. Cambridge University Press. 1235p.
[2] Carlson BC (1977) Elliptic Integrals of the First Kind, SIAM Journal on Mathematical
    Analysis, vol. 8, pp. 231-242.

func CarlsonRj 

func CarlsonRj(x, y, z, p float64) float64

CarlsonRj computes Carlson’s elliptic integral of the third kind according to [1] Computes Rj(x,y,z,p) where x,y,z must be nonnegative, and at most one can be zero. p must be nonzero. If p < 0, the Cauchy principal value is returned.

func ChebyshevT 

func ChebyshevT(n int, x float64) float64

ChebyshevT directly computes the Chebyshev polynomial of first kind Tn(x) using the trigonometric functions. 

        │ (-1)ⁿ cosh[n⋅acosh(-x)]   if x < -1
Tn(x) = ┤       cosh[n⋅acosh( x)]   if x > 1
        │       cos [n⋅acos ( x)]   if |x| ≤ 1

func ChebyshevTdiff1 added in v1.0.1 

func ChebyshevTdiff1(n int, x float64) float64

ChebyshevTdiff1 computes the first derivative of the Chebyshev function Tn(x) 

dTn
————
 dx

func ChebyshevTdiff2 added in v1.0.1 

func ChebyshevTdiff2(n int, x float64) float64

ChebyshevTdiff2 computes the second derivative of the Chebyshev function Tn(x) 

d²Tn
—————
 dx²

func ChebyshevXgauss added in v1.0.1 

func ChebyshevXgauss(N int) (X []float64)

ChebyshevXgauss computes Chebyshev-Gauss roots considering symmetry 

            /  (2i+1)⋅π  \
X[i] = -cos | —————————— |       i = 0 ... N
            \   2N + 2   /

func ChebyshevXlob added in v1.0.1 

func ChebyshevXlob(N int) (X []float64)

ChebyshevXlob computes Chebyshev-Gauss-Lobatto points using the sin function and considering symmetry 

                    /  π⋅(N-2i)  \
        X[i] = -sin | —————————— |       i = 0 ... N
                    \    2⋅N     /
  or
                    /  i⋅π  \
        X[i] = -cos | ————— |       i = 0 ... N
                    \   N   /

Reference:
[1] Baltensperger R and Trummer MR (2003) Spectral differencing with a twist, SIAM J. Sci.
    Comput. 24(5):1465-1487

func Dft1d 

func Dft1d(data []complex128, inverse bool)

Dft1d computes the discrete Fourier transform (DFT) in 1D. It replaces data by its discrete Fourier transform, if inverse==false or replaces data by its inverse discrete Fourier transform, if inverse==true 

Computes:
                   N-1         -i 2 π j k / N                 __
  forward:  X[k] =  Σ  x[j] ⋅ e                     with i = √-1
                   j=0

                   N-1         +i 2 π j k / N
  inverse:  Y[k] =  Σ  y[j] ⋅ e                     thus x[k] = Y[k] / N
                   j=0

NOTE: (1) the inverse operation does not divide by N
      (2) ideally, N=len(data) is an integer power of 2.
      (3) using FFTW: http://fftw.org/fftw3_doc/What-FFTW-Really-Computes.html

func Elliptic1 

func Elliptic1(φ, k float64) float64

Elliptic1 computes Legendre elliptic integral of the first kind F(φ,k), evaluated using Carlson’s function Rf [1]. The argument ranges are 0 ≤ φ ≤ π/2 and 0 ≤ k·sin(φ) ≤ 1 

Computes:
                   φ
                  ⌠          dt
      F(φ, k)  =  │  ___________________
                  │     _______________
                  ⌡   \╱ 1 - k² sin²(t)
                 0
where:
         0 ≤ φ ≤ π/2
         0 ≤ k·sin(φ) ≤ 1

References:
[1] Press WH, Teukolsky SA, Vetterling WT, Fnannery BP (2007) Numerical Recipes: The Art of
    Scientific Computing. Third Edition. Cambridge University Press. 1235p.

func Elliptic2 

func Elliptic2(φ, k float64) float64

Elliptic2 computes Legendre elliptic integral of the second kind E(φ,k), evaluated using Carlson's functions Rf and Rd [1]. The argument ranges are 0 ≤ φ ≤ π/2 and 0 ≤ k⋅sin(φ) ≤ 1 

Computes:
                   φ
                  ⌠     _______________
      E(φ, k)  =  │   \╱ 1 - k² sin²(t)  dt
                  ⌡
                 0
where:
         0 ≤ φ ≤ π/2
         0 ≤ k·sin(φ) ≤ 1

References:
[1] Press WH, Teukolsky SA, Vetterling WT, Fnannery BP (2007) Numerical Recipes: The Art of
    Scientific Computing. Third Edition. Cambridge University Press. 1235p.

func Elliptic3 

func Elliptic3(n, φ, k float64) float64

Elliptic3 computes Legendre elliptic integral of the third kind Π(n,φ,k), evaluated using Carlson's functions Rf and Rj. NOTE that the sign convention on n corresponds to that of Abramowitz and Stegun [2] and not to [1]. The argument ranges are 0 ≤ φ ≤ π/2 and 0 ≤ k⋅sin(φ) ≤ 1 

Computes:
                      φ
                     ⌠                  dt
      Π(n, φ, k)  =  │  ___________________________________
                     │                     _______________
                     ⌡   (1 - n sin²(t)) \╱ 1 - k² sin²(t)
                    0
where:
         0 ≤ φ ≤ π/2
         0 ≤ k·sin(φ) ≤ 1

References:
[1] Press WH, Teukolsky SA, Vetterling WT, Fnannery BP (2007) Numerical Recipes: The Art of
    Scientific Computing. Third Edition. Cambridge University Press. 1235p.
[2] Abramowitz M, Stegun IA (1972) Handbook of Mathematical Functions with Formulas, Graphs,
    and Mathematical Tables. U.S. Department of Commerce, NIST

func ExpMix 

func ExpMix(x float64) complex128

ExpMix uses Euler's formula to compute exp(-i⋅x) = cos(x) - i⋅sin(x)

func ExpPix 

func ExpPix(x float64) complex128

ExpPix uses Euler's formula to compute exp(+i⋅x) = cos(x) + i⋅sin(x)

func Factorial100 

func Factorial100(n int) big.Float

Factorial100 returns the factorial n! up to 100! using the math/big package

func Factorial22 

func Factorial22(n int) float64

Factorial22 implements the factorial function; i.e. computes n! up to 22! According to [1], factorials up to 22! have exact double precision representations (52 bits of mantissa, not counting powers of two that are absorbed into the exponent) 

References
[1] Press WH, Teukolsky SA, Vetterling WT, Fnannery BP (2007) Numerical Recipes: The Art of
     Scientific Computing. Third Edition. Cambridge University Press. 1235p.

func Hat added in v1.0.1 

func Hat(x, xc, y0, h, l float64) float64

Hat implements the hat function 

   --———--   o (xc,y0+h)
      |     / \
      h    /   \    m = h/l
      |   /m    \
y0 ——————o       o—————————

         |<  2l >|

func HatD1 added in v1.0.1 

func HatD1(x, xc, y0, h, l float64) float64

HatD1 returns the first derivative of the hat function NOTE: the discontinuity is ignored ⇒ D1(xc-l)=D1(xc+l)=D1(xc)=0

func Heav 

func Heav(x float64) float64

Heav computes the Heaviside step function (== derivative of Ramp(x)) 

          │ 0    if x < 0
Heav(x) = ┤ 1/2  if x = 0
          │ 1    if x > 0

func ImagPowN added in v1.0.1 

func ImagPowN(n int) complex128

ImagPowN computes iⁿ = (√-1)ⁿ 

i¹ = i      i²  = -1      i³  = -i      i⁴  = 1
i⁵ = i      i⁶  = -1      i⁷  = -i      i⁸  = 1
i⁹ = i      i¹⁰ = -1      i¹¹ = -i      i¹² = 1

func ImagXpowN added in v1.0.1 

func ImagXpowN(x float64, n int) complex128

ImagXpowN computes (x⋅i)ⁿ 

(x⋅i)¹ = x¹⋅i      (x⋅i)²  = -x²       (x⋅i)³  = -x³ ⋅i      (x⋅i)⁴  = x⁴
(x⋅i)⁵ = x⁵⋅i      (x⋅i)⁶  = -x⁶       (x⋅i)⁷  = -x⁷ ⋅i      (x⋅i)⁸  = x⁸
(x⋅i)⁹ = x⁹⋅i      (x⋅i)¹⁰ = -x¹⁰      (x⋅i)¹¹ = -x¹¹⋅i      (x⋅i)¹² = x¹²

func Logistic added in v1.1.0 

func Logistic(z float64) float64

Logistic implements the sigmoid/logistic function

func LogisticD1 added in v1.1.0 

func LogisticD1(z float64) float64

LogisticD1 implements the first derivative of the sigmoid/logistic function

func ModBesselI0 

func ModBesselI0(x float64) (ans float64)

ModBesselI0 returns the modified Bessel function I0(x) for any real x.

func ModBesselI1 

func ModBesselI1(x float64) (ans float64)

ModBesselI1 returns the modified Bessel function I1(x) for any real x.

func ModBesselIn 

func ModBesselIn(n int, x float64) (ans float64)

ModBesselIn returns the modified Bessel function In(x) for any real x and n ≥ 0

func ModBesselK0 

func ModBesselK0(x float64) float64

ModBesselK0 returns the modified Bessel function K0(x) for positive real x. 

Special cases:
  K0(x=0) = +Inf
  K0(x<0) = NaN

func ModBesselK1 

func ModBesselK1(x float64) float64

ModBesselK1 returns the modified Bessel function K1(x) for positive real x. 

Special cases:
  K0(x=0) = +Inf
  K0(x<0) = NaN

func ModBesselKn 

func ModBesselKn(n int, x float64) float64

ModBesselKn returns the modified Bessel function Kn(x) for positive x and n ≥ 0

func NegOnePowN added in v1.0.1 

func NegOnePowN(n int) float64

NegOnePowN computes (-1)ⁿ

func Pow2 added in v1.1.0 

func Pow2(x float64) float64

Pow2 computes x²

func Pow3 added in v1.1.0 

func Pow3(x float64) float64

Pow3 computes x³

func PowP added in v1.0.1 

func PowP(x float64, n uint32) (r float64)

PowP computes real raised to positive integer xⁿ

func Ramp 

func Ramp(x float64) float64

Ramp function => MacAulay brackets

func Rbinomial 

func Rbinomial(x, y float64) float64

Rbinomial computes the binomial coefficient with real (non-negative) arguments by calling the Gamma function

func Rect 

func Rect(x float64) float64

Rect implements the rectangular function 

Rect(x) = Boxcar(x;-0.5,0.5)

          │ 0    if |x| > 1/2
Rect(x) = ┤ 1/2  if |x| = 1/2
          │ 1    if |x| < 1/2

func Sabs 

func Sabs(x, eps float64) float64

Sabs implements a smooth abs function: Sabs(x) = x*x / (sign(x)*x + eps)

func SabsD1 

func SabsD1(x, eps float64) float64

SabsD1 returns the first derivative of Sabs

func SabsD2 

func SabsD2(x, eps float64) float64

SabsD2 returns the first derivative of Sabs

func Sign 

func Sign(x float64) float64

Sign implements the sign function 

          │ -1   if x < 0
Sign(x) = ┤  0   if x = 0
          │  1   if x > 0

func Sinc 

func Sinc(x float64) float64

Sinc computes the sine cardinal (sinc) function 

Sinc(x) = |     1      if x = 0
          | sin(x)/x   otherwise

func Sramp 

func Sramp(x, β float64) float64

Sramp implements a smooth ramp function. Ramp

func SrampD1 

func SrampD1(x, β float64) float64

SrampD1 returns the first derivative of Sramp

func SrampD2 

func SrampD2(x, β float64) float64

SrampD2 returns the second derivative of Sramp

func SuqCos 

func SuqCos(angle, expon float64) float64

SuqCos implements the superquadric auxiliary function that uses cos(x)

func SuqSin 

func SuqSin(angle, expon float64) float64

SuqSin implements the superquadric auxiliary function that uses sin(x)

func UintBinomial 

func UintBinomial(n, k uint64) uint64

UintBinomial implements the Binomial coefficient using uint64. Panic happens on overflow Also, this function uses a loop so it may not be very efficient for large k The code below comes from https://en.wikipedia.org/wiki/Binomial_coefficient [cannot find a reference to cite]

Types

type Axis added in v1.1.0 

type Axis struct {

	// configuration data
	DisableHunt bool // do not use hunt code at all
	// contains filtered or unexported fields
}

Axis implements a type to hold an arbitrarily spaced discrete data

func NewAxis added in v1.1.0 

func NewAxis(data []float64, interpType InterpType) (o *Axis)

NewAxis builds a new Axis type from a data slice for an InterpType

func (*Axis) Get added in v1.1.0 

func (o *Axis) Get(i int) float64

Get returns the value at data[i]

type BiLinear added in v1.1.0 

type BiLinear struct {
	// contains filtered or unexported fields
}

BiLinear implements a two dimensional interpolant

func NewBiLinear added in v1.1.0 

func NewBiLinear(f, xx, yy []float64) (o *BiLinear)

NewBiLinear builds a two dimensional bi-linear interpolant Input: 

xx -- function sample points abscissas
yy -- function sample points ordinates
f  -- function values
f(i,j) is stored at f[len(xx)*j + i]

Ref: 

f(x,y) = x^2 + 2y^2         2 h    i    j
xx = [0.00,0.50,1.00]         |
yy = [0.00,1.00,2.00]       1 d    e    f
f  = [a:0.00,b:0.25,c:1.00,   |
      d:2.00,e:2,25,f:3.00,   a____b____c
      h:8.00,i:8.25,j:9.00]  0   0.5   1.0

func (*BiLinear) P added in v1.1.0 

func (o *BiLinear) P(x, y float64) float64

P is the interpolation polynomial

func (*BiLinear) Reset added in v1.1.0 

func (o *BiLinear) Reset(f, xx, yy []float64)

Reset (Re)Set the axis and matrix for the interpolant

func (*BiLinear) SetDisableHunt added in v1.1.0 

func (o *BiLinear) SetDisableHunt(disable bool)

SetDisableHunt disables the hunt function for both axis

type ChebyInterp added in v1.0.1 

type ChebyInterp struct {

	// input
	N     int  // degree of polynomial
	Gauss bool // use roots (Gauss) or points (Lobatto)?

	// options
	Bary  bool // [default=true] use barycentric formulae in for ℓ_i and I{f} [default=true]
	Nst   bool // [default=true] use the "negative sum trick" to compute the diagonal components according to:
	Trig  bool // [default=false] use trigonometric identities (to reduce round-off errors)
	Flip  bool // [default=false] compute lower-diagonal part from upper diagonal part with D_{N-j,N-l} = -D_{j,l}
	StdD2 bool // [default=false] compute D2 using standard formula, otherwise use D1

	// constants
	EstimationN int // N to use when estimating CoefP [default=128]

	// derived
	X     []float64 // points. NOTE: mirrowed version of Chebyshev X; i.e. from +1 to -1
	Wb    []float64 // weights for Gaussian quadrature
	Gamma []float64 // denominator of coefficients equation ~ ‖p[i]‖²
	Lam   []float64 // λ_i barycentric weights (also w_i in some papers)

	// computed by auxiliary methods
	CoefI  []float64 // coefficients of interpolant
	CoefP  []float64 // coefficients of projection (estimated)
	CoefIs []float64 // coefficients of interpolation using Lagrange cardinal functions

	// computed
	C  *la.Matrix // physical to transform space conversion matrix
	Ci *la.Matrix // transform to physical space conversion matrix
	D1 *la.Matrix // (dℓj/dx)(xi)
	D2 *la.Matrix // (d²ℓj/dx²)(xi)
}

ChebyInterp defines a structure for efficient computations with Chebyshev polynomials such as projection or interpolation 

Some equations are based on [1,2,3]

References:
  [1] Canuto C, Hussaini MY, Quarteroni A, Zang TA (2006) Spectral Methods: Fundamentals in
      Single Domains. Springer. 563p
  [2] Webb M, Trefethen LN, Gonnet P (2012) Stability of barycentric interpolation formulas for
      extrapolation, SIAM J. Sci. Comput. Vol. 34, No. 6, pp. A3009-A3015
  [3] Baltensperger R, Trummer M (2003) Spectral Differencing with a twist,
      SIAM J. Sci. Comput., Vol. 24, No. 5, pp. 1465-1487

func NewChebyInterp added in v1.0.1 

func NewChebyInterp(N int, gaussChebyshev bool) (o *ChebyInterp)

NewChebyInterp returns a new ChebyInterp structure 

gaussChebyshev == true:

                  /  (2⋅j+1)⋅π  \
        X_j = cos | ——————————— |       j = 0 ... N
                  \   2⋅N + 2   /

gaussChebyshev == false: (Gauss-Lobatto-Chebyshev)

                  /  j⋅π  \
        X_j = cos | ————— |       j = 0 ... N
                  \   N   /

NOTE: X here is the mirrowed version of Chebyshev X; i.e. from +1 to -1

func (*ChebyInterp) CalcCoefI added in v1.0.1 

func (o *ChebyInterp) CalcCoefI(f Ss)

CalcCoefI computes the coefficients of the interpolant by (slow) direct formula 

           1    N
CoefI_k = ——— ⋅ Σ  f(x_j) ⋅ T_k(x_j) ⋅ wb_j
          γ_k  j=0

Thus (for Gauss-Lobatto):

            2       N    1                   /  k⋅j⋅π  \
CoefI_k = —————— ⋅  Σ  —————— ⋅ f(x_j) ⋅ cos | ——————— |
          N⋅cb_k   j=0  cb_j                 \    N    /

where:
        cb_k = 2 if j=0,N   or   1 if j=1...N-1

NOTE: the results will be stored in o.CoefI

func (*ChebyInterp) CalcCoefIs added in v1.0.1 

func (o *ChebyInterp) CalcCoefIs(f Ss)

CalcCoefIs computes the coefficients for interpolation with Lagrange cardinal functions ℓ_l(x)

func (*ChebyInterp) CalcCoefP added in v1.0.1 

func (o *ChebyInterp) CalcCoefP(f Ss)

CalcCoefP computes the coefficients of the projection (slow) using o.EstimationN + 1 points 

            ∫ f(x)⋅T_k(x)⋅w(x) dx      (f, T_k)_w
CoefP_k = ————————————————————————— = ————————————
           ∫ T_k(x)⋅T_k(x)⋅w(x) dx      ‖ T_k ‖²

NOTE: the results will be stored in o.CoefP

func (*ChebyInterp) CalcConvMats added in v1.0.1 

func (o *ChebyInterp) CalcConvMats()

CalcConvMats computes conversion matrices C and Ci 

              N
trans(u)_k =  Σ  C_{kj} ⋅ u(x_j)      e.g. trans(u)_k == coefI_k
             j=0

       N
u_j =  Σ   C⁻¹_{jk} ⋅ trans(u)_k
      k=0

               2             /  j⋅k⋅π  \
C_{kj} = ————————————— ⋅ cos | ——————— |
          cb_k⋅cb_j⋅N        \    N    /

               /  j⋅k⋅π  \
C⁻¹_{jk} = cos | ——————— |
               \    N    /

func (*ChebyInterp) CalcD1 added in v1.0.1 

func (o *ChebyInterp) CalcD1()

CalcD1 computes the differentiation matrix D1 of the function L_i 

 d I{f}(x)     N            d ℓ_l(x)
——————————— =  Σ   f(x_l) ⋅ ————————
     dx       l=0              dx

 d I{f}(x)  |         N
——————————— |      =  Σ   D1_jl ⋅ f(x_l)
     dx     |x=x_j   l=0

where:

         dℓ_l  |
 D1_jl = ————— |
          dx   |x=x_j

Equations (2.4.31) and (2.4.33), page 89 of [1]

If Nst==true (negative-sum-trick):

                                  N
                       D_{jj} = - Σ  D_{jl}
                                 l=0
                                 l≠j

NOTE: (1) the signs are swapped (compared to [1]) because X are reversed here (from -1 to +1)
      (2) this method is only available for Gauss-Lobatto points

func (*ChebyInterp) CalcD2 added in v1.0.1 

func (o *ChebyInterp) CalcD2()

CalcD2 calculates the second derivative 

          d²ℓ_l  |
  D2_jl = —————— |
           dx²   |x=x_j

NOTE: this function will call CalcD1() because the D1 values required to compute D2,
      unless StdD2=true where the "standard" formula (Eq. 2.4.32) is used instead => less accurate

func (*ChebyInterp) CalcErrorD1 added in v1.0.1 

func (o *ChebyInterp) CalcErrorD1(dfdxAna Ss) (maxDiff float64)

CalcErrorD1 computes the maximum error due to differentiation (@ X[i]) using the D1 matrix 

NOTE: CoefIs and D1 matrix must be computed previously

func (*ChebyInterp) CalcErrorD2 added in v1.0.1 

func (o *ChebyInterp) CalcErrorD2(d2fdx2Ana Ss) (maxDiff float64)

CalcErrorD2 computes the maximum error due to differentiation (@ X[i]) using the D2 matrix 

NOTE: CoefIs and D2 matrix must be computed previously

func (*ChebyInterp) EstimateMaxErr added in v1.0.1 

func (o *ChebyInterp) EstimateMaxErr(f Ss, projection bool) (maxerr, xloc float64)

EstimateMaxErr estimates the maximum error using 10000 stations along [-1,1] This function also returns the location (xloc) of the estimated max error 

maxerr = max(|f - I{f}|)  or  maxerr = max(|f - P{f}|)

NOTE: CoefI or CoefP must be computed first

func (*ChebyInterp) HierarchicalT added in v1.0.1 

func (o *ChebyInterp) HierarchicalT(i int, x float64) float64

HierarchicalT computes Tn(x) using hierarchical definition (but NOT recursive) 

NOTE: this function is not as efficient as ChebyshevT and should be used for testing only

func (*ChebyInterp) I added in v1.0.1 

func (o *ChebyInterp) I(x float64) (res float64)

I computes the interpolation 

           N
I{f}(x) =  Σ  CoefI_k ⋅ T_k(x)
          k=0

Thus:

          N
f(x_j) =  Σ   CoefI_k ⋅ T_k(x_j)
         k=0

NOTE: CoefI coefficients must be computed first

func (*ChebyInterp) Il added in v1.0.1 

func (o *ChebyInterp) Il(x float64) (res float64)

Il computes the interpolation using the Lagrange cardinal functions ℓ_i(x) 

           N                                         N
I{f}(x) =  Σ   f(x_i) ⋅ ℓ_i(x)    or      I{f}(x) =  Σ  CoefIs_i ⋅ ℓ_i(x)
          l=0                                       i=0

NOTE: (1) CoefIs == f(x_i) coefficients must be computed (or set) first
      (2) ℓ is symbolized by ℓ in [1]

func (*ChebyInterp) L added in v1.0.1 

func (o *ChebyInterp) L(i int, x float64) float64

L evaluates the Lagrange cardinal function ℓ_i(x) of degree N with Gauss-Lobatto points 

           N
I{f}(x) =  Σ   f(x_i) ⋅ ℓ_i(x)
          i=0

Equation (2.4.30), page 88 of [1]

NOTE: must not use with Gauss (roots) points

func (*ChebyInterp) P added in v1.0.1 

func (o *ChebyInterp) P(x float64) (res float64)

P computes the (approximated) projection 

           ∞
S{f}(x) =  Σ  CoefP_k ⋅ T_k(x)   (series representation)
          k=0

Thus:

           N
P{f}(x) =  Σ  CoefP_k ⋅ T_k(x)   (truncated series)
          k=0

NOTE: CoefP coefficients must be computed first

type DataInterp added in v1.0.1 

type DataInterp struct {

	// configuration data
	DisableHunt bool // do not use hunt code at all

	// output data
	Dy float64 // error estimate
	// contains filtered or unexported fields
}

DataInterp implements numeric interpolators to be used with discrete data

func NewDataInterp added in v1.0.1 

func NewDataInterp(Type string, p int, xx, yy []float64) (o *DataInterp)

NewDataInterp creates new interpolator for data point sets xx and yy (with same lengths) 

Type -- type of interpolator
   "lin"  : linear
   "poly" : polynomial

p  -- order of interpolator
xx -- x-data
yy -- y-data

func (*DataInterp) P added in v1.0.1 

func (o *DataInterp) P(x float64) float64

P computes P(x); i.e. performs the interpolation

func (*DataInterp) Reset added in v1.0.1 

func (o *DataInterp) Reset(xx, yy []float64)

Reset re-assigns xx and yy data sets

type FourierInterp added in v1.0.1 

type FourierInterp struct {

	// main
	N int        // number of terms. must be power of 2; i.e. N = 2ⁿ
	X la.Vector  // point coordinates == 2⋅π.j/N
	K la.Vector  // k values computed from j such that j = 0...N-1 ⇒ k = -N/2...N/2-1
	A la.VectorC // coefficients for interpolation. from FFT
	S la.VectorC // smoothing coefficients

	// computed (U may be set externally)
	U      la.Vector  // values of f(x) at grid points (nodes) X[j]
	Du     la.Vector  // p-order derivative of u
	Du1    la.Vector  // 1st derivative of f(x) at grid points (nodes) X[j]
	Du2    la.Vector  // 2nd derivative of f(x) at grid points (nodes) X[j]
	DuHat  la.VectorC // spectral coefficient corresponding to p-derivative
	Du1Hat la.VectorC // spectral coefficient corresponding to 1st derivative
	Du2Hat la.VectorC // spectral coefficient corresponding to 1st derivative
	// contains filtered or unexported fields
}

FourierInterp performs interpolation using truncated Fourier series 

            N/2 - 1
             ————          +i k X[j]
  f(x[j]) =  \     A[k] ⋅ e                   with    X[j] = 2 π j / N
             /
             ————
            k = -N/2                 Eq (2.1.27) of [1]    x ϵ [0, 2π]

  where:

              N - 1
          1   ————             -i k X[j]
  A[k] = ———  \     f(x[j]) ⋅ e              with    X[j] = 2 π j / N
          N   /
              ————
             j = 0                                  Eq (2.1.25) of [1]

NOTE: (1) f=u in [1] and A[k] is the tilde(u[k]) of [1]
      (2) FFTW says "so you should initialize your input data after creating the plan."
          Therefore, the plan can be created and reused several times.
          [http://www.fftw.org/fftw3_doc/Planner-Flags.html]
          Also: "The plan can be reused as many times as needed. In typical high-performance
          applications, many transforms of the same size are computed"
          [http://www.fftw.org/fftw3_doc/Introduction.html]

Create a new object with NewFourierInterp(...) AND deallocate memory with Free()

Reference:
  [1] Canuto C, Hussaini MY, Quarteroni A, Zang TA (2006) Spectral Methods: Fundamentals in
      Single Domains. Springer. 563p

func NewFourierInterp added in v1.0.1 

func NewFourierInterp(N int, smoothing string) (o *FourierInterp)

NewFourierInterp allocates a new FourierInterp object 

N -- number of terms. must be even; ideally power of 2, e.g. N = 2ⁿ

smoothing -- type of smoothing: use SmoNoneKind for no smoothing
  "" or "none" : no smoothing
  "lanc"       : Lanczos (sinc)
  "rcos"       : Raised Cosine
  "ces"        : Cesaro

NOTE: remember to call Free in the end to release memory allocated by FFTW; e.g.
      defer o.Free()

func (*FourierInterp) CalcA added in v1.0.1 

func (o *FourierInterp) CalcA()

CalcA calculates the coefficients A of the interpolation using (fwd) FFT 

               N - 1
           1   ————             -i k X[j]
   A[k] = ———  \     f(x[j]) ⋅ e              with    X[j] = 2 π j / N
           N   /
               ————
              j = 0                                  Eq (2.1.25) of [1]

NOTE: remember to set U (or call CalcU) first

func (*FourierInterp) CalcAwithAliasRemoval added in v1.0.1 

func (o *FourierInterp) CalcAwithAliasRemoval(f Ss)

CalcAwithAliasRemoval calculates the coefficients A by using the 3/2-rule to remove alias error via the padding method 

NOTE: with the 3/2-rule, the intepolatory property is not exact; i.e. I(xi)≈f(xi) only

func (*FourierInterp) CalcD added in v1.0.1 

func (o *FourierInterp) CalcD(p int)

CalcD calculates the p-derivative of the interpolated function @ grid points using the FFT (with smoothing or not) 

                p      |
               d(I{f}) |
       dfdx =  ——————— |             len(res) must be equal to N
                   p   |
                 dx    |x=x[j]

 INPUT:
    p -- derivative order

 OUTPUT:
    Du and DuHat will contain the results

NOTE: remember to call CalcA first

func (*FourierInterp) CalcD1 added in v1.0.1 

func (o *FourierInterp) CalcD1()

CalcD1 calculates the 1st derivative using function CalcD and internal arrays. See function CalcD for further details 

OUTPUT: the results will be stored in Du1 and Du1Hat

func (*FourierInterp) CalcD2 added in v1.0.1 

func (o *FourierInterp) CalcD2()

CalcD2 calculates the 2nd derivative using function CalcD and internal arrays. See function CalcD for further details 

OUTPUT: the results will be stored in Du2 and Du2Hat

func (*FourierInterp) CalcJ added in v1.0.1 

func (o *FourierInterp) CalcJ(k float64) int

CalcJ computes j-index from k-index where j corresponds to the FFT index 

k ϵ [-N/2, N/2-1]
j ϵ [0, N-1]

Example with N = 8:

     k=0 ⇒ j=0      k=-4 ⇒ j=4
     k=1 ⇒ j=1      k=-3 ⇒ j=5      j = { N + k  if  k < 0
     k=2 ⇒ j=2      k=-2 ⇒ j=6          {     k  otherwise
     k=3 ⇒ j=3      k=-1 ⇒ j=7

func (*FourierInterp) CalcK added in v1.0.1 

func (o *FourierInterp) CalcK(j int) float64

CalcK computes k-index from j-index where j corresponds to the FFT index 

FFT returns the A coefficients as:

   {A[0], A[1], ..., A[N/2-1], A[-N/2], A[-N/2+1], ... A[-1]}

k ϵ [-N/2, N/2-1]
j ϵ [0, N-1]

Example with N = 8:

     j=0 ⇒ k=0      j=4 ⇒ k=-4
     j=1 ⇒ k=1      j=5 ⇒ k=-3
     j=2 ⇒ k=2      j=6 ⇒ k=-2
     j=3 ⇒ k=3      j=7 ⇒ k=-1

func (*FourierInterp) CalcU added in v1.0.1 

func (o *FourierInterp) CalcU(f Ss)

CalcU calculates f(x) at grid points (to be used later with CalcA and/or CalcD)

func (*FourierInterp) Free added in v1.0.1 

func (o *FourierInterp) Free()

Free releases resources allocated for FFTW

func (*FourierInterp) I added in v1.0.1 

func (o *FourierInterp) I(x float64) float64

I computes the interpolation (with smoothing or not) 

             N/2 - 1
               ————          +i k x
   I {f}(x) =  \     A[k] ⋅ e                 x ϵ [0, 2π]
    N          /
               ————
              k = -N/2                 Eq (2.1.28) of [1]

NOTE: remember to call CalcA first

func (*FourierInterp) Idiff added in v1.0.1 

func (o *FourierInterp) Idiff(p int, x float64) float64

Idiff performs the differentiation of the interpolation; i.e. computes the p-derivative of the interpolation (with smoothing or not) 

                 p       N/2 - 1
      p         d(I{f})    ————       p           +i k x
res: DI{f}(x) = ——————— =  \     (i⋅k)  ⋅ A[k] ⋅ e
      N             p      /
                  dx       ————
                          k = -N/2                   x ϵ [0, 2π]

NOTE: remember to call CalcA first

type GeneralOrthoPoly 

type GeneralOrthoPoly struct {

	// input
	Kind string // type of orthogonal polynomial
	N    int    // (max) degree of polynomial. Lower order can be quickly obtained after this polynomial with max(N) is generated
	// contains filtered or unexported fields
}

GeneralOrthoPoly implements general orthogonal polynomials. It uses a general format and is NOT very efficient for large degrees. For efficiency, use the OrthoPoly structure instead. 

Reference:
[1] Abramowitz M, Stegun IA (1972) Handbook of Mathematical Functions with Formulas, Graphs,
    and Mathematical Tables. U.S. Department of Commerce, NIST

NOTE: this structure should be not be used for high-performance computing;
      it's probably useful for verifications or learning purposes only

func NewGeneralOrthoPoly 

func NewGeneralOrthoPoly(kind string, N int, alpha, beta float64) (o *GeneralOrthoPoly)

NewGeneralOrthoPoly creates a new orthogonal polynomial 

kind -- is the type of orthognal polynomial:
  "J" or "jac"    : Jacobi
  "L" or "leg"    : Legendre
  "H" or "her"    : Hermite
  "T" or "cheby1" : Chebyshev first kind
  "U" or "cheby2" : Chebyshev second kind

N -- is the (max) degree of the polynomial.
     Lower order can later be quickly obtained after this
     polynomial with max(N) is created

alpha -- Jacobi only: α coefficient

beta -- Jacobi only: β coefficient

NOTE: all coefficients for the 0...N polynomials will be generated

NOTE: this structure should be not be used for high-performance computing;
      it's probably useful for verifications or learning purposes only

func (*GeneralOrthoPoly) F 

func (o *GeneralOrthoPoly) F(x float64) (res float64)

F computes P(n,x) with n=N (max) 

Since GeneralOrthoPoly is a general form, the summations are directly implement; i.e. no
advantages are taken w.r.t the structure of the polynomial. Thus, these functions are not
highly efficient for large degrees N

func (*GeneralOrthoPoly) P 

func (o *GeneralOrthoPoly) P(n int, x float64) (res float64)

P computes P(n,x) where n must be ≤ N 

Since GeneralOrthoPoly is a general form, the summations are directly implement; i.e. no
advantages are taken w.r.t the structure of the polynomial. Thus, these functions are not
highly efficient for large degrees N

type InterpCubic added in v1.1.0 

type InterpCubic struct {
	A, B, C, D float64 // coefficients of polynomial
	TolDen     float64 // tolerance to avoid zero denominator
}

InterpCubic computes a cubic polynomial to perform interpolation either using 4 points or 3 points and a known derivative

func NewInterpCubic added in v1.1.0 

func NewInterpCubic() (o *InterpCubic)

NewInterpCubic returns a new object

func (*InterpCubic) Critical added in v1.1.0 

func (o *InterpCubic) Critical() (xmin, xmax, xifl float64, hasMin, hasMax, hasIfl bool)

Critical returns the critical points 

xmin -- x @ min and y(xmin)
xmax -- x @ max and y(xmax)
xifl -- x @ inflection point and y(ifl)
hasMin, hasMax, hasIfl -- flags telling what is available

func (*InterpCubic) F added in v1.1.0 

func (o *InterpCubic) F(x float64) float64

F computes y = f(x) curve

func (*InterpCubic) Fit3pointsD added in v1.1.0 

func (o *InterpCubic) Fit3pointsD(x0, y0, x1, y1, x2, y2, x3, d3 float64) (err error)

Fit3pointsD fits polynomial to 3 points and known derivative 

(x0, y0) -- first point
(x1, y1) -- second point
(x2, y2) -- third point
(x3, d3) -- derivative @ x3

func (*InterpCubic) Fit4points added in v1.1.0 

func (o *InterpCubic) Fit4points(x0, y0, x1, y1, x2, y2, x3, y3 float64) (err error)

Fit4points fits polynomial to 3 points 

(x0, y0) -- first point
(x1, y1) -- second point
(x2, y2) -- third point
(x3, y3) -- fourth point

func (*InterpCubic) G added in v1.1.0 

func (o *InterpCubic) G(x float64) float64

G computes y' = df/x|(x) curve

type InterpQuad added in v1.1.0 

type InterpQuad struct {
	A, B, C float64 // coefficients of polynomial
	TolDen  float64 // tolerance to avoid zero denominator
}

InterpQuad computes a quadratic polynomial to perform interpolation either using 3 points or 2 points and a known derivative

func NewInterpQuad added in v1.1.0 

func NewInterpQuad() (o *InterpQuad)

NewInterpQuad returns a new object

func (*InterpQuad) F added in v1.1.0 

func (o *InterpQuad) F(x float64) float64

F computes y = f(x) curve

func (*InterpQuad) Fit2pointsD added in v1.1.0 

func (o *InterpQuad) Fit2pointsD(x0, y0, x1, y1, x2, d2 float64) (err error)

Fit2pointsD fits polynomial to 2 points and known derivative 

(x0, y0) -- first point
(x1, y1) -- second point
(x2, d2) -- derivative @ x2

func (*InterpQuad) Fit3points added in v1.1.0 

func (o *InterpQuad) Fit3points(x0, y0, x1, y1, x2, y2 float64) (err error)

Fit3points fits polynomial to 3 points 

(x0, y0) -- first point
(x1, y1) -- second point
(x2, y2) -- third point

func (*InterpQuad) G added in v1.1.0 

func (o *InterpQuad) G(x float64) float64

G computes y' = df/x|(x) curve

func (*InterpQuad) Optimum added in v1.1.0 

func (o *InterpQuad) Optimum() (xopt, fopt float64)

Optimum returns the minimum or maximum point; i.e. the point with zero derivative 

xopt -- x @ optimum
fopt -- f(xopt) = y @ optimum

type InterpType added in v1.1.0 

type InterpType int

InterpType specifies the type of interpolant 

const (

	// BiLinearType defines the bi-linear type
	BiLinearType InterpType = 2

	// BiCubicType defines the bi-cubic type
	BiCubicType InterpType = 3
)

type LagIntSet added in v1.1.0 

type LagIntSet []*LagrangeInterp

LagIntSet is groups interpolators together; e.g. 2D, 3D

func NewLagIntSet added in v1.1.0 

func NewLagIntSet(ndim int, degrees []int, gridTypes []string) (lis LagIntSet)

NewLagIntSet returns a set of LagrangeInterp

type LagrangeInterp 

type LagrangeInterp struct {

	// general
	N int       // degree: N = len(X)-1
	X la.Vector // grid points: len(X) = P+1; generated in [-1, 1]
	U la.Vector // function evaluated @ nodes: f(x_i)

	// barycentric
	Bary   bool      // [default=true] use barycentric weights
	UseEta bool      // [default=true] use ηk when computing D1
	Eta    la.Vector // sum of log of differences: ηk = Σ ln(|xk-xl|) (k≠l)
	Lam    la.Vector // normalized barycentric weights λk = pow(-1, k+N) ⋅ ηk / (2ⁿ⁻¹/n)

	// computed
	D1 *la.Matrix // (dℓj/dx)(xi)
	D2 *la.Matrix // (d²ℓj/dx²)(xi)
}

LagrangeInterp implements Lagrange interpolators associated with a grid X 

An interpolant I^X_N{f} (associated with a grid X; of degree N; with N+1 points)
is expressed in the Lagrange form as follows:

                  N
      X          ————             X
     I {f}(x) =  \     f(x[i]) ⋅ ℓ (x)
      N          /                i
                 ————
                 i = 0

where ℓ^X_i(x) is the i-th Lagrange cardinal polynomial associated with grid X and given by:

              N
      N      ━━━━    x  -  X[j]
     ℓ (x) = ┃  ┃  —————————————           0 ≤ i ≤ N
      i      ┃  ┃   X[i] - X[j]
            j = 0
            j ≠ i

or, barycentric form:

                  N   λ[i] ⋅ f[i]
                  Σ   ———————————
      X          i=0   x - x[i]
     I {f}(x) = ——————————————————
      N            N     λ[i]
                   Σ   ————————
                  i=0  x - x[i]

with:

                  λ[i]
                ————————
      N         x - x[i]
     ℓ (x) = ———————————————
      i        N     λ[k]
               Σ   ————————
              k=0  x - x[k]

The barycentric weights λk are normalized and computed from ηk as follows:

   ηk = Σ ln(|xk-xl|) (k≠l)

         a ⋅ b             k+N
   λk =  —————     a = (-1)        b = exp(m)    m = -ηk
          lf0

   lf0 = 2ⁿ⁻¹/n

 or, if N > 700:

         / a ⋅ b \   /  b  \   /  b  \
   λk =  | ————— | ⋅ | ——— | ⋅ | ——— |      b = exp(m/3)
         \  lf0  /   \ lf1 /   \ lf2 /

   lf0⋅lf1⋅lf2 = 2ⁿ⁻¹/n

References:
  [1] Canuto C, Hussaini MY, Quarteroni A, Zang TA (2006) Spectral Methods: Fundamentals in
      Single Domains. Springer. 563p
  [2] Berrut JP, Trefethen LN (2004) Barycentric Lagrange Interpolation,
      SIAM Review Vol. 46, No. 3, pp. 501-517
  [3] Costa B, Don WS (2000) On the computation of high order pseudospectral derivatives,
      Applied Numerical Mathematics, 33:151-159.

func NewLagrangeInterp 

func NewLagrangeInterp(N int, gridType string) (o *LagrangeInterp)

NewLagrangeInterp allocates a new LagrangeInterp 

N -- degree

gridType -- type of grid:
   "uni" : uniform 1D grid kind
   "cg"  : Chebyshev-Gauss 1D grid kind
   "cgl" : Chebyshev-Gauss-Lobatto 1D grid kind

NOTE: the grid will be generated in [-1, 1]

func (*LagrangeInterp) CalcD1 added in v1.0.1 

func (o *LagrangeInterp) CalcD1()

CalcD1 computes the differentiation matrix D1 of the function L_i 

 d I{f}(x)  |         N
——————————— |      =  Σ   D1_kj ⋅ f(x_j)
     dx     |x=x_k   j=0

see [2]

func (*LagrangeInterp) CalcD2 added in v1.0.1 

func (o *LagrangeInterp) CalcD2()

CalcD2 calculates the second derivative 

          d²ℓ_l  |
  D2_jl = —————— |
           dx²   |x=x_j

NOTE: this function will call CalcD1() because the D1 values required to compute D2

func (*LagrangeInterp) CalcErrorD1 added in v1.0.1 

func (o *LagrangeInterp) CalcErrorD1(dfdxAna Ss) (maxDiff float64)

CalcErrorD1 computes the maximum error due to differentiation (@ X[i]) using the D1 matrix 

NOTE: U and D1 matrix must be computed previously

func (*LagrangeInterp) CalcErrorD2 added in v1.0.1 

func (o *LagrangeInterp) CalcErrorD2(d2fdx2Ana Ss) (maxDiff float64)

CalcErrorD2 computes the maximum error due to differentiation (@ X[i]) using the D2 matrix 

NOTE: U and D2 matrix must be computed previously

func (*LagrangeInterp) CalcU added in v1.0.1 

func (o *LagrangeInterp) CalcU(f Ss)

CalcU computes f(x_i); i.e. function f(x) @ all nodes

func (*LagrangeInterp) EstimateLebesgue 

func (o *LagrangeInterp) EstimateLebesgue() (ΛN float64)

EstimateLebesgue estimates the Lebesgue constant by using 10000 stations along [-1,1]

func (*LagrangeInterp) EstimateMaxErr 

func (o *LagrangeInterp) EstimateMaxErr(nStations int, f Ss) (maxerr, xloc float64)

EstimateMaxErr estimates the maximum error using 10000 stations along [-1,1] This function also returns the location (xloc) of the estimated max error 

Computes:
          maxerr = max(|f(x) - I{f}(x)|)

e.g. nStations := 10000 (≥2) will generate several points along [-1,1]

func (*LagrangeInterp) I 

func (o *LagrangeInterp) I(x float64) (res float64)

I computes the interpolation I^X_N{f}(x) @ x 

                  N
      X          ————          X
     I {f}(x) =  \     U[i] ⋅ ℓ (x)       with   U[i] = f(x[i])
      N          /             i
                 ————
                 i = 0

or (barycentric):

                 N   λ[i] ⋅ f[i]
                 Σ   ———————————
     X          i=0   x - x[i]
    I {f}(x) = ——————————————————
     N            N     λ[i]
                  Σ   ————————
                 i=0  x - x[i]

NOTE: U[i] = f(x[i]) must be calculated with o.CalcU or set first

func (*LagrangeInterp) L 

func (o *LagrangeInterp) L(i int, x float64) (lix float64)

L computes the i-th Lagrange cardinal polynomial ℓ^X_i(x) associated with grid X 

              N
      X      ━━━━    x  -  X[j]
     ℓ (x) = ┃  ┃  —————————————           0 ≤ i ≤ N
      i      ┃  ┃   X[i] - X[j]
            j = 0
            j ≠ i

or (barycentric):

                 λ[i]
               ————————
     X         x - x[i]
    ℓ (x) = ———————————————
     i        N     λ[k]
              Σ   ————————
             k=0  x - x[k]

Input:
   i -- index of X[i] point
   x -- where to evaluate the polynomial
Output:
   lix -- ℓ^X_i(x)

func (*LagrangeInterp) Om added in v1.0.1 

func (o *LagrangeInterp) Om(x float64) (ω float64)

Om computes the generating (nodal) polynomial associated with grid X. The nodal polynomial is the unique polynomial of degree N+1 and leading coefficient whose zeros are the N+1 nodes of X. 

         N
 X      ━━━━
ω (x) = ┃  ┃ (x - X[i])
N+1     ┃  ┃
       i = 0

type Mm 

type Mm func(f, m *la.Matrix)

Mm defines a matrix function f(m) of a matrix argument m (matrix matrix)) 

Input:
  m -- input matrix
Output:
  M -- output matrix

type Mss added in v1.1.0 

type Mss func(m *la.Matrix, a, b float64)

Mss defines a matrix function f(a,b) of two scalar arguments (matrix scalar scalar) 

Input:
  a -- first input scalar
  b -- second input scalar
Output:
  m -- output matrix

type Mv 

type Mv func(f *la.Matrix, v la.Vector)

Mv defines a matrix function f(v) of a vector argument v (matrix vector)) 

Input:
  v -- input vector
Output:
  f -- output matrix

type Sinusoid 

type Sinusoid struct {

	// input
	Period     float64 // T: period; e.g. [s]
	MeanValue  float64 // A0: mean value; e.g. [m]
	Amplitude  float64 // C1: amplitude; e.g. [m]
	PhaseShift float64 // θ: phase shift; e.g. [rad]

	// derived
	Frequency   float64 // f: frequency; e.g. [Hz] or [1 cycle/s]
	AngularFreq float64 // ω0 = 2⋅π⋅f: angular frequency; e.g. [rad⋅s⁻¹]
	TimeShift   float64 // ts = θ / ω0: time shift; e.g. [s]

	// derived: coefficients
	A []float64 // A0, A1, A2, ... (if series mode)
	B []float64 // B0, B1, B2, ... (if series mode)
}

Sinusoid implements the sinusoid equation: 

y(t) = A0 + C1⋅cos(ω0⋅t + θ)             [essential-form]

y(t) = A0 + A1⋅cos(ω0⋅t) + B1⋅sin(ω0⋅t)  [basis-form]

A1 =  C1⋅cos(θ)
B1 = -C1⋅sin(θ)
θ  = arctan(-B1 / A1)   if A1<0, θ += π
C1 = sqrt(A1² + B1²)

func NewSinusoidBasis 

func NewSinusoidBasis(T, A0, A1, B1 float64) (o *Sinusoid)

NewSinusoidBasis creates a new Sinusoid object with the "basis" parameters set 

T  -- period; e.g. [s]
A0 -- mean value; e.g. [m]
A1 -- coefficient of the cos term
B1 -- coefficient of the sin term

func NewSinusoidEssential 

func NewSinusoidEssential(T, A0, C1, θ float64) (o *Sinusoid)

NewSinusoidEssential creates a new Sinusoid object with the "essential" parameters set 

T  -- period; e.g. [s]
A0 -- mean value; e.g. [m]
C1 -- amplitude; e.g. [m]
θ  -- phase shift; e.g. [rad]

func (*Sinusoid) ApproxSquareFourier 

func (o *Sinusoid) ApproxSquareFourier(N int)

ApproxSquareFourier approximates sinusoid using Fourier series with N terms

func (*Sinusoid) TestPeriodicity 

func (o *Sinusoid) TestPeriodicity(tmin, tmax float64, npts int) bool

TestPeriodicity tests that f(t) = f(T + t)

func (*Sinusoid) Ybasis 

func (o *Sinusoid) Ybasis(t float64) (res float64)

Ybasis computes y(t) = A0 + A1⋅cos(ω0⋅t) + B1⋅sin(ω0⋅t) [basis-form]

func (*Sinusoid) Yessen 

func (o *Sinusoid) Yessen(t float64) float64

Yessen computes y(t) = A0 + C1⋅cos(ω0⋅t + θ [essential-form]

type Ss 

type Ss func(s float64) float64

Ss defines a scalar function f(s) of a scalar argument s (scalar scalar) 

Input:
  s -- input scalar
Returns:
  scalar

type Sss added in v1.1.0 

type Sss func(s1, s2 float64) float64

Sss defines a scalar function f(r,s) of two scalar arguments (scalar scalar scalar) 

Input:
  s1, s2 -- input scalar
Returns:
  scalar

type Sv 

type Sv func(v la.Vector) float64

Sv defines a scalar function f(v) of a vector argument v (scalar vector) 

Input:
  v -- input vector
Returns:
  scalar

type Svs added in v1.0.1 

type Svs func(v la.Vector, s float64) float64

Svs defines a scalar function f(v,s) of a vector and a scalar 

Input:
  s -- the scalar
  v -- the vector
Returns:
  scalar

type Tt 

type Tt func(f, t *la.Triplet)

Tt defines a triplet (matrix) function f(t) of a triplet (matrix) argument t (triplet triplet) 

Input:
  t -- input triplet
Output:
  f -- output triplet

type Tv 

type Tv func(f *la.Triplet, v la.Vector)

Tv defines a triplet (matrix) function f(v) of a vector argument v (triplet vector) 

Input:
  v -- input vector
Output:
  f -- output triplet

type Vs 

type Vs func(f la.Vector, s float64)

Vs defines a vector function f(s) of a scalar argument s (vector scalar) 

Input:
  s -- input scalar
Output:
  f -- output vector

type Vss added in v1.1.0 

type Vss func(f la.Vector, a, b float64)

Vss defines a vector function f(a,b) of two scalar arguments (vector scalar scalar) 

Input:
  a -- first input scalar
  b -- second input scalar
Output:
  f -- output vector

type Vv 

type Vv func(f, v la.Vector)

Vv defines a vector function f(v) of a vector argument v (vector vector) 

Input:
  v -- input vector
Output:
  f -- output vector

type Vvss added in v1.1.0 

type Vvss func(u, v la.Vector, a, b float64)

Vvss defines two vector functions u(a,b) and v(a,b) of two scalar arguments (vector vector scalar scalar) 

Input:
  a -- first input scalar
  b -- second input scalar
Output:
  u -- first output vector
  v -- second output vector

type Vvvss added in v1.1.0 

type Vvvss func(u, v, w la.Vector, a, b float64)

Vvvss defines three vector functions u(a,b), v(a,b) and w(a,b) of two scalar arguments (vector vector vector scalar scalar) 

Input:
  a -- first input scalar
  b -- second input scalar
Output:
  u -- first output vector
  v -- second output vector
  w -- second output vector

Source Files

View all Source files

Directories

Path Synopsis
Package fftw wraps the FFTW library to perform Fourier Transforms using the "fast" method by Cooley and Tukey
 Package fftw wraps the FFTW library to perform Fourier Transforms using the "fast" method by Cooley and Tukey

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL
go.dev uses cookies from Google to deliver and enhance the quality of its services and to analyze traffic. Learn more.