fun #

fn atan2p(y f64, x f64) f64

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

fn atan2pdeg(y f64, x f64) f64

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

fn bessel_j0(x_ f64) f64

bessel_j0 returns the order-zero Bessel function of the first kind.

special cases are: bessel_j0(±inf) = 0 bessel_j0(0) = 1 bessel_j0(nan) = nan

fn bessel_j1(x_ f64) f64

bessel_j1 returns the order-one Bessel function of the first kind.

special cases are: bessel_j1(±inf) = 0 bessel_j1(nan) = nan

fn bessel_jn(n_ int, x_ f64) f64

bessel_jn returns the order-n Bessel function of the first kind.

special cases are: bessel_jn(n, ±inf) = 0 bessel_jn(n, nan) = nan

fn bessel_y0(x f64) f64

bessel_y0 returns the order-zero Bessel function of the second kind.

special cases are: bessel_y0(+inf) = 0 bessel_y0(0) = -inf bessel_y0(x < 0) = nan bessel_y0(nan) = nan

fn bessel_y1(x f64) f64

bessel_y1 returns the order-one Bessel function of the second kind.

special cases are: bessel_y1(+inf) = 0 bessel_y1(0) = -inf bessel_y1(x < 0) = nan bessel_y1(nan) = nan

fn bessel_yn(n_ int, x f64) f64

bessel_yn returns the order-n Bessel function of the second kind.

special cases are: bessel_yn(n, +inf) = 0 bessel_yn(n ≥ 0, 0) = -inf bessel_yn(n < 0, 0) = +inf if n is odd, -inf if n is even bessel_yn(n, x < 0) = nan bessel_yn(n, nan) = nan

fn beta(a f64, b f64) f64

beta computes the beta function by calling the log_gamma_sign function

fn binomial(n int, k_ int) f64

binomial comptues the binomial coefficient (n k)^t

fn boxcar(x f64, a f64, b f64) f64

boxcar implements the boxcar function

boxcar(x;a,b) = heav(x-a) - heav(x-b)

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

fn cgamma(z cmplx.Complex) cmplx.Complex

gamma computes the gamma function value

fn choose(n int, p int) f64

Compute the binomial coefficient

fn clog_gamma(z cmplx.Complex) cmplx.Complex

log_gamma computes the log-gamma function value

fn cos(x f64) f64

fn cos_e(x f64) (f64, f64)

fn digamma(x f64) f64

fn erf(x_ f64) f64

erf returns the error function of x.

special cases are: erf(+inf) = 1 erf(-inf) = -1 erf(nan) = nan

fn erfc(x_ f64) f64

erfc returns the complementary error function of x.

special cases are: erfc(+inf) = 0 erfc(-inf) = 2 erfc(nan) = nan

fn exp_mix(x f64) cmplx.Complex

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

fn exp_pix(x f64) cmplx.Complex

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

fn fib(n int) int

fib returns the nth number in the Fibonacci sequence using O(1) space in O(logn) time complexity

fn gamma(x_ f64) f64

gamma returns the gamma function of x.

special ifs are: gamma(+inf) = +inf gamma(+0) = +inf gamma(-0) = -inf gamma(x) = nan for integer x < 0 gamma(-inf) = nan gamma(nan) = nan

fn hat(x f64, xc f64, y0 f64, h f64, l f64) f64

hat implements the hat function

fn hatd1(x f64, xc f64, y0 f64, h f64, l f64) f64

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

fn heav(x f64) f64

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

fn hypot(x f64, y f64) f64

fn hypot_e(x f64, y f64) (f64, f64)

fn imag_pow_n(n int) cmplx.Complex

imag_pow_n computes iⁿ = (√-1)ⁿ

fn imag_x_pow_n(x f64, n int) cmplx.Complex

imag_x_pow_n computes (x⋅i)ⁿ

fn lin_interp(mut o DataInterp, j int, x f64) f64

lin_interp implements linear interpolator

fn log_gamma(x f64) f64

log_gamma returns the natural logarithm and sign (-1 or +1) of Gamma(x).

special ifs are: log_gamma(+inf) = +inf log_gamma(0) = +inf log_gamma(-integer) = +inf log_gamma(-inf) = -inf log_gamma(nan) = nan

fn log_gamma_sign(x_ f64) (f64, int)

fn logistic(z f64) f64

logistic implements the sigmoid/logistic function

fn logistic_d1(z f64) f64

logistic_d1 implements the first derivative of the sigmoid/logistic function

fn n_combos_w_replacement(n f64, r f64) u64

(n+r-1)! / r! / (n-1)! when n > 0.

fn neg_one_pow_n(n int) f64

neg_one_pow_n computes (-1)ⁿ

fn poly_interp(mut o DataInterp, jl int, x f64) f64

poly_interp performs a polynomial interpolation. this routine returns an interpolated value y, and stores an error estimate dy. the returned value is obtained by m-point polynomial interpolation on the subrange xx[jl..jl+m-1].

fn pone(x f64) f64

For x >= 8, the asymptotic expansions of pone is 1 + 15/128 s2 - 4725/215 s4 - ..., where s = 1/x. We approximate pone by pone(x) = 1 + (R/S) where R = pr0 + pr1*s2 + pr2s**4 + ... + pr5s10 S = 1 + ps0*s2 + ... + ps4*s10 and | pone(x)-1-R/S | <= 2(-60.06)

fn pow2(x f64) f64

pow2 computes x²

fn pow3(x f64) f64

pow3 computes x³

fn powp(x_ f64, n u32) f64

powp computes real raised to positive integer xⁿ

fn psi(x_ f64) f64

fn pzero(x f64) f64

The asymptotic expansions of pzero is 1 - 9/128 s2 + 11025/98304 s4 - ..., where s = 1/x. For x >= 2, We approximate pzero by pzero(x) = 1 + (R/S) where R = pj0r0 + pR1s**2 + pR2s4 + ... + pR5*s10 S = 1 + pj0s0s**2 + ... + pS4s**10 and | pzero(x)-1-R/S | <= 2 ** ( -60.26)

fn qone(x f64) f64

For x >= 8, the asymptotic expansions of qone is 3/8 s - 105/1024 s3 - ..., where s = 1/x. We approximate qone by qone(x) = s*(0.375 + (R/S)) where R = qr1*s2 + qr2s**4 + ... + qr5s10 S = 1 + qs1*s2 + ... + qs6*s12 and | qone(x)/s -0.375-R/S | <= 2(-61.13)

fn qzero(x f64) f64

For x >= 8, the asymptotic expansions of qzero is -1/8 s + 75/1024 s3 - ..., where s = 1/x. We approximate pzero by qzero(x) = s*(-1.25 + (R/S)) where R = qj0r0 + qR1*s2 + qR2s**4 + ... + qR5s10 S = 1 + qj0s0*s2 + ... + qS5*s12 and | qzero(x)/s +1.25-R/S | <= 2(-61.22)

fn ramp(x f64) f64

ramp function => macaulay brackets

fn rbinomial(x f64, y f64) f64

rbinomial computes the binomial coefficient with real (non-negative) arguments by calling the gamma function

fn rect(x f64) f64

rect implements the rectangular function

rect(x) = boxcar(x;-0.5,0.5)

fn sabs(x f64, eps f64) f64

sabs implements a smooth abs f: sabs(x) = x*x / (sign(x)*x + eps)

fn sabs_d1(x f64, eps f64) f64

sabs_d1 returns the first derivative of sabs

fn sabs_d2(x f64, eps f64) f64

sabs_d2 returns the first derivative of sabs

fn sign(x f64) f64

sign implements the sign function

fn sin(x f64) f64

fn sin_e(x f64) (f64, f64)

fn sinc(x f64) f64

sinc computes the sine cardinal (sinc) function

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

fn sramp(x f64, beta f64) f64

sramp implements a smooth ramp function. ramp

fn srampd1(x f64, beta f64) f64

srampd1 returns the first derivative of sramp

fn srampd2(x f64, beta f64) f64

srampd2 returns the second derivative of sramp

fn suqcos(angle f64, expon f64) f64

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

fn suqsin(angle f64, expon f64) f64

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

fn uint_binomial(n_ u64, k_ u64) u64

uint_binomial implements the binomial coefficient using u64. 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]

fn DataInterp.new(itype string, p int, xx []f64, yy []f64) &DataInterp

DataInterp.new 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

fn InterpCubic.new() &InterpCubic

InterpCubic.new returns a new InterpCubic instance

fn InterpQuad.new() &InterpQuad

InterpQuad.new returns a new InterpQuad instance

fn Sinusoid.basis(t f64, a0 f64, a1 f64, b1 f64) &Sinusoid

Sinusoid.basis 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

fn Sinusoid.essential(t f64, a0 f64, c1 f64, theta f64) &Sinusoid

Sinusoid.essential 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] theta -- phase shift; e.g. [rad]

type InterpFn = fn (mut o DataInterp, j int, x f64) f64

InterpFn defines the type of the implementation of the data interpolation

struct ChebSeries {
pub:
	c     []f64 // coefficients
	order int   // order of expansion
	a     f64   // lower interval point
	b     f64   // upper interval point
}

data for a Chebyshev series over a given interval

fn (cs ChebSeries) eval_e(x f64) (f64, f64)

struct DataInterp {
mut:
	// input data
	itype string // type of interpolator
	xx    []f64  // x-data values
	yy    []f64  // y-data values
	// derived data
	m        int  // number of points of interpolating formula; e.g. 2 for segments, 3 for 2nd order polynomials
	n        int  // length of xx
	j_hunt   int  // temporary j to decide on using hunt
	dj_hunt  int  // increent of j to decide on using hunt function or locate
	use_hunt bool // use hunt code instead of locate
	ascnd    bool // ascending order of x-values
	// implementation
	interp InterpFn = unsafe { nil }
pub mut:
	// configuration data
	disable_hunt bool // do not use hunt code at all
	// output data
	dy f64 // error estimate
}

DataInterp implements numeric interpolators to be used with discrete data

fn (mut o DataInterp) reset(xx []f64, yy []f64)

reset re-assigns xx and yy data sets

fn (mut o DataInterp) p(x f64) f64

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

fn (mut o DataInterp) locate(x f64) int

locate returns a value j such that x is (insofar as possible) centered in the subrange xx[j..j+mm-1], where xx is the stored pointer. the values in xx must be monotonic, either increasing or decreasing. the returned value is not less than 0, nor greater than n-1.

fn (mut o DataInterp) hunt(x f64) int

hunt returns a value j such that x is (insofar as possible) centered in the subrange xx[j..j+mm-1], where xx is the stored pointer. the values in xx must be monotonic, either increasing or decreasing. the returned value is not less than 0, nor greater than n-1.

struct InterpCubic {
pub mut:
	// coefficients of polynomial a, b, c, d
	a f64
	b f64
	c f64
	d f64
	// tolerance to avoid zero denominator
	tol_den f64
}

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

fn (o InterpCubic) f(x f64) f64

f computes y = f(x) curve

fn (o InterpCubic) g(x f64) f64

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

fn (o InterpCubic) critical() (f64, f64, f64, bool, bool, 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) has_min, has_max, has_ifl -- flags telling what is available

fn (mut o InterpCubic) fit_4points(x0 f64, y0 f64, x1 f64, y1 f64, x2 f64, y2 f64, x3 f64, y3 f64) !

fit_4points fits polynomial to 3 points (x0, y0) -- first point (x1, y1) -- second point (x2, y2) -- third point (x3, y3) -- fourth point

fn (mut o InterpCubic) fit_3points_d(x0 f64, y0 f64, x1 f64, y1 f64, x2 f64, y2 f64, x3 f64, d3 f64) !

fit_3points_d fits polynomial to 3 points and known derivative (x0, y0) -- first point (x1, y1) -- second point (x2, y2) -- third point (x3, d3) -- derivative @ x3

struct InterpQuad {
pub mut:
	// coefficients of polynomial a, b, c
	a f64
	b f64
	c f64
	// tolerance to avoid zero denominator
	tol_den f64
}

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

fn (o InterpQuad) f(x f64) f64

f computes y = f(x) curve

fn (o InterpQuad) g(x f64) f64

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

fn (o InterpQuad) optimum() (f64, f64)

optimum returns the minimum or maximum point; i.e. the point with zero derivative xopt -- x @ optimum fopt -- f(xopt) = y @ optimum

fn (mut o InterpQuad) fit_3points(x0 f64, y0 f64, x1 f64, y1 f64, x2 f64, y2 f64) !

fit_3points fits polynomial to 3 points (x0, y0) -- first point (x1, y1) -- second point (x2, y2) -- third point

fn (mut o InterpQuad) fit_2points_d(x0 f64, y0 f64, x1 f64, y1 f64, x2 f64, d2 f64) !

fit_2points_d fits polynomial to 2 points and known derivative (x0, y0) -- first point (x1, y1) -- second point (x2, d2) -- derivate @ x2

struct Sinusoid {
pub mut:
	// input
	period      f64 // t: period; e.g. [s]
	mean_value  f64 // a0: mean value; e.g. [m]
	amplitude   f64 // c1: amplitude; e.g. [m]
	phase_shift f64 // theta: phase shift; e.g. [rad]
	// derived
	frequency    f64 // f: frequency; e.g. [Hz] or [1 cycle/s]
	angular_freq f64 // omega_0 = 2⋅π⋅f: angular frequency; e.g. [rad⋅s⁻¹]
	time_shift   f64 // ts = theta / omega_0: time shift; e.g. [s]
	// derived: coefficients
	a []f64 // a0, a1, A2, ... (if series mode)
	b []f64 // B0, b1, B2, ... (if series mode)
}

Sinusoid implements the sinusoid equation:

y(t) = a0 + c1⋅cos(omega_0⋅t + theta) [essential-form]

y(t) = a0 + a1⋅cos(omega_0⋅t) + b1⋅sin(omega_0⋅t) [basis-form]

a1 = c1⋅cos(theta) b1 = -c1⋅sin(theta) theta = arctan(-b1 / a1) if a1<0, theta += π c1 = sqrt(a1² + b1²)

fn (o Sinusoid) yessen(t f64) f64

yessen computes y(t) = a0 + c1⋅cos(omega_0⋅t + theta [essential-form]

fn (o Sinusoid) ybasis(t f64) f64

ybasis computes y(t) = a0 + a1⋅cos(omega_0⋅t) + b1⋅sin(omega_0⋅t) [basis-form]

fn (mut o Sinusoid) approx_square_fourier(n int)

approx_square_fourier approximates sinusoid using Fourier series with n terms

fn (o Sinusoid) test_periodicity(tmin f64, tmax f64, npts int) bool

test_periodicity tests that f(t) = f(t + t)