#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Mathematical equations.
:Date: 2019-11-15
.. module:: function
:platform: *nix, Windows
:synopsis: Mathematical equations.
.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
Functions returns function to apply conditions.
"""
import math
from .data import seq
from .mathematics import lcm
[docs]def sine_wave(A=1, k=1, f=1, phi=0, D=0, degree=False):
r"""A sine wave or sinusoid is a mathematical curve that describes a
smooth periodic oscillation.
:param A: amplitude
:type A: float or int
:param k: (angular) wave number
:type k: float or int
:param f: ordinary frequency
:type f: float or int
:param phi: phase
:type phi: float or int
:param D: non-zero center amplitude
:type D: float or int
:param degree: boolean to switch between radians and degree. If
False phi is interpreted in radians and if True then phi is
interpreted in degrees.
:type degree: bool
:results: sine wave function of spatial variable x and optional
time t
:rtype: function
In general, the function is:
.. math::
y(x,t) = A\sin(kx + \omega t + \varphi) + D \\
y(x,t) = A\sin(kx + 2\pi f t + \varphi) + D
where:
* A, amplitude, the peak deviation of the function from zero.
* f, ordinary frequency, the number of oscillations (cycles) that
occur each second of time.
* ω = 2πf, angular frequency, the rate of change of the function
argument in units of radians per second. If ω < 0 the wave is
moving to the right, if ω > 0 the wave is moving to the left.
* φ, phase, specifies (in radians) where in its cycle the
oscillation is at t = 0.
* x, spatial variable that represents the position on the
dimension on which the wave propagates.
* k, characteristic parameter called wave number (or angular wave
number), which represents the proportionality between the
angular frequency ω and the linear speed (speed of propagation)
ν.
* D, non-zero center amplitude.
The wavenumber is related to the angular frequency by:
.. math::
k={\omega \over v}={2\pi f \over v}={2\pi \over \lambda }
where λ (lambda) is the wavelength, f is the frequency, and v is the
linear speed.
.. seealso::
:meth:`cosine_wave`
"""
if degree:
phi = math.radians(phi)
return lambda x, t=0: A*math.sin(k*x + 2*math.pi*f*t + phi) + D
[docs]def cosine_wave(A=1, k=1, f=1, phi=0, D=0, degree=False):
r"""A cosine wave is said to be sinusoidal, because,
:math:`\cos(x) = \sin(x + \pi/2)`, which is also a sine wave with a
phase-shift of π/2 radians. Because of this head start, it is often
said that the cosine function leads the sine function or the sine
lags the cosine.
:param A: amplitude
:type A: float or int
:param k: (angular) wave number
:type k: float or int
:param f: ordinary frequency
:type f: float or int
:param phi: phase
:type phi: float or int
:param D: non-zero center amplitude
:type D: float or int
:param degree: boolean to switch between radians and degree. If
False phi is interpreted in radians and if True then phi is
interpreted in degrees.
:type degree: bool
:results: sine wave function of spatial variable x and optional
time t
:rtype: function
.. seealso::
:meth:`sine_wave`
"""
if degree:
phi = phi + 90
else:
phi = phi + math.pi/2
return sine_wave(A=A, k=k, f=f, phi=phi, D=D, degree=degree)
[docs]def b_spline_basis(knots, knot_span, degree):
r"""Cox-de Boor algorithm / recursion formula.
Calculate the i-th B-spline basis function of degree p: N_{i,p}(u)
:param knots: Knot vector U. m + 1 non-decreasing numbers / knots,
:math:`u_0 <= u_1 <= u_2 <= ... <= u_m`
:type knots: list
:param knot_span: i-th knot span
:type knot_span: int
:param degree: degree of B-spline basis function
:type degree: int
:returns: B-spline basis function using variable, u \in [u_0, u_m]
:rtype: function
.. math::
N_{i,0}(u) &= \begin{cases} 1 & \text{if } u_i \le u \lt u_{i+1} \\
0 & \text{otherwise}\end{cases} \\
N_{i,p}(u) &= \frac{u - u_i}{u_{i+p} - u_i} N_{i,p-1}(u) +
\frac{u_{i+p+1} - u}{u_{i+p+1} - u_{i+1}} N_{i+1,p-1}(u)
"""
U = knots
N = b_spline_basis
i = knot_span
p = degree
if p == 0:
def Np(u):
return 1 if U[i] <= u < U[i+1] else 0
else:
def Np(u):
try:
term1 = (u - U[i])/(U[i+p] - U[i])*N(U, i, p-1)(u)
except:
term1 = 0
try:
term2 = (U[i+p+1] - u)/(U[i+p+1] - U[i+1])*N(U, i+1, p-1)(u)
except:
term2 = 0
return term1 + term2
return Np
# TODO: other types than clamped?
[docs]def b_spline_knots(control_point_spans, degree=3):
"""B-spline knots to generate a clamped uniform B-spline curve
of degree p (order + 1).
The internal knots are equally spaced (uniform B-spline curve)
:param control_point_spans: number of control points + 1
:type control_point_spans: int
:param degree: degree of B-spline basis functions (default = 3)
:type degree: int
:returns: knot vector
:rtype: tuple
.. seealso::
:meth:`b_spline_curve_with_knots`
"""
p = degree
n = control_point_spans
m = n + p + 1 # number of knot spans
# number of
U_outer = p + 1 # at each of the vector
U_inner = m + 1 - 2*(U_outer)
U = [0]*(U_outer)
U += [i for i in range(1, U_inner)]
U += [U_inner]*(U_outer)
return tuple(U) # tuples are hashable
[docs]def b_spline_curve_with_knots(degree, control_points, knots):
"""B-spline curve of degree p (order + 1) on a given set of knots.
n, m and p must satisfy m = n + p + 1.
:param degree: degree of B-spline basis functions
:type degree: int
:param control_points: control points P, n + 1 control points
:type control_points: list
:param knots: Knot vector U. m + 1 non-decreasing numbers / knots,
:math:`u_0 <= u_1 <= u_2 <= ... <= u_m`
:type knots: list
:returns: B-spline curve using variable, u \in [u_0, u_m]
:rtype: function
.. math::
\mathbf{C}_p(u) = \sum\limits_{i=0}^{n} N_{i,p}(u) \mathbf{P}_i
* open B-spline curves
* the curve will not touch the first and last legs of the
control polyline
* the knot vector does not have any particular structure
* for degree p, intervals [u_0, u_p) and [u_{n-p}, u_n) will not
have "full support" of basis functions and are ignored when a
B-spline curve is open. For open B-spline curves, the domain
is inteval [u_p, u_{n-p}]
* clamped B-spline curves, nonperiodic B-spline curves
* the curve is tangent to the first and the last legs just like
a Bézier curve
* the first knot and the last knot must be repeated p+1 times
(i. e., of multiplicity p+1)
* closed B-spline curves
* the start and the end of the generated curve join together
forming a closed loop
* repeating some knots and control points # TODO: which?
* uniform B-spline curves
* internal knots are equally spaced
* Bézier curves
* a B-spline with no internal knots.
.. seealso::
:meth:`b_spline_knots`
"""
dim = len(control_points[0])
def C(u):
NiPi = [0] * dim
for i in range(len(control_points)):
Ni = b_spline_basis(knots, i, degree)
for j in range(dim):
NiPi[j] += Ni(u)*control_points[i][j]
return NiPi
return C
[docs]def sample_half_open(f, a, b, n=50, endpoint_epsilon=1e-7):
# hack to sample close to the endpoint
x = seq(a, b, (b-a)/n) + [b-endpoint_epsilon]
# Sample the function
return [f(xi) for xi in x]
[docs]def sample_half_open_seq(f, x, endpoint_epsilon=1e-7):
# hack to sample close to the endpoint, x[-1] can be present
# multiple times.
xend = x[-1]
x = [xi for xi in x if xi != xend] + [xend-endpoint_epsilon]
# Sample the function
return [f(xi) for xi in x]
#
# Parametric equations
# roulette
#
[docs]def circle(r):
r"""Circle
:param r: radius of the circle
:type r: float
:results: functions for x of theta and y of theta and the interval
:rtype: tuple
.. math::
x(\theta) = r\cos\theta \\
y(\theta) = r\sin\theta \\
\theta = \left[0, 2\pi\right]
::
* *
* r *
* *
* *
* *
* *
>>> x, y = circle(20)[:2]
>>> x, y, _ = circle(20)
>>> x, y, interval = circle(20)
"""
return ellipse(r, r)
[docs]def ellipse(a, b):
r"""Ellipse
:param a: semi-major axis
:type a: float
:param b: semi-minor axis
:type b: float
:results: functions for x of theta and y of theta and the interval
:rtype: tuple
.. math::
x(\theta) = a\cos\theta \\
y(\theta) = b\sin\theta \\
\theta = \left[0, 2\pi\right]
::
* .*
* b : *
* :......*
* a *
* *
* *
>>> x, y = ellipse(10, 5)[:2]
>>> x, y, _ = ellipse(10, 5)
>>> x, y, interval = ellipse(10, 5)
"""
x = lambda theta: a*math.cos(theta)
y = lambda theta: b*math.sin(theta)
return x, y, [0, 2*math.pi]
[docs]def hypotrochoid(R, r, d):
r"""Hypotrochoid
A point is attached with a distance d from the center of a circle
of radius r. The circle is rolling around the inside of a fixed
circle of radius R.
:param R: radius of the fixed exterior circle
:type R: float
:param r: radius of the rolling circle inside of the fixed circle
:typre r: float
:param d: distance from the center of the interior circle
:type d: float
:results: functions for x of theta and y of theta and the interval
:rtype: tuple
.. math::
x(\theta) = (R - r)\cos\theta + d\cos\left(\frac{R-r}{r}\theta\right) \\
y(\theta) = (R - r)\sin\theta - d\sin\left(\frac{R-r}{r}\theta\right) \\
\theta = \left[0, 2\pi\frac{\mathrm{lcm}(r, R)}{R}\right]
::
* * *
* R *
* *
* * * *
* * r **
* * ,.. *
* * d *
* * **
* * * *
* *
* *
* * *
>>> x, y = hyotrochoid(20, 6, 6)[:2]
>>> x, y, _ = hyotrochoid(20, 6, 6)
>>> x, y, interval = hyotrochoid(20, 6, 6)
.. seealso::
:meth:`pylib.mathematics.lcm`
"""
x = lambda theta: (R-r)*math.cos(theta) + d*math.cos((R-r)/r * theta)
y = lambda theta: (R-r)*math.sin(theta) - d*math.sin((R-r)/r * theta)
return x, y, [0, 2*math.pi*lcm(r, R)/R]
[docs]def epitrochoid(R, r, d):
r"""Epitrochoid
A point is attached with a distance d from the center of a circle
of radius r. The circle is rolling around the outside of a fixed
circle of radius R.
:param R: radius of the fixed interior circle
:type R: float
:param r: radius of the rolling circle outside of the fixed circle
:typre r: float
:param d: distance from the center of the exterior circle
:type d: float
:results: functions for x of theta and y of theta and the interval
:rtype: tuple
.. math::
x(\theta) = (R + r)\cos\theta - d\cos\left(\frac{R+r}{r}\theta\right) \\
y(\theta) = (R + r)\sin\theta - d\sin\left(\frac{R+r}{r}\theta\right) \\
\theta = \left[0, 2\pi\right]
::
* * *
* R *
* *
* * * *
* * * r *
* ** .., *
* ** d *
* * * *
* * * *
* *
* *
* * *
>>> x, y = epitrochoid(3, 1, 0.5)[:2]
>>> x, y, _ = epitrochoid(3, 1, 0.5)
>>> x, y, interval = epitrochoid(3, 1, 0.5)
"""
x = lambda theta: (R+r)*math.cos(theta) - d*math.cos((R+r)/r * theta)
y = lambda theta: (R+r)*math.sin(theta) - d*math.sin((R+r)/r * theta)
return x, y, [0, 2*math.pi]
[docs]def to_str(f, x=None, x_0=0, x_1=1, t=None, h=10, w=80, density=1, char_set="line"):
"""Represent functions as string frame with a specific character set.
which are normed to the range of [0, 1] to
:param f: function or list of functions normed to the range of [0, 1]
:type f: function or list
:param h: number of chars in vertical direction
:type h: int
:param w: number of chars in horizontal direction
:type w: int
:param char_set: either "braille" or "block". "braille" uses Unicode
Characters in the Braille Patterns Block (first index U+2800, last
index U+28FF [CUDB]_) and the "block" uses part of the Unicode
Characters in the Block Elements Block (fisrt index U+2580, last
index U+259F [CUDB]_). Alias for braille is line and alias for
block is histogram
:type char_set: str
Usage:
* case dependent arguments
* 1 quasi line plot (braille characters)
f = lambda function (lambda x, t=0: ...)
x_1 = max x value
x_0 = min x value
t = time (animation)
* 2 histogram (block characters)
f = lambda function (lambda x, t=0: ...)
x_1 = max x value
x_0 = min x value
t = time (animation)
chart="histogram"
* general arguments
w = window width in number of chars
density = number of data point per pixel (for line chart higher density makes thicker lines)
chart = either "line" or "histogram"
.. rubric:: Braille Patterns Block
* Dots or pixels per character in vertical direction: 6
* Dots or pixels per character in horizontal direction: 2
First dot (bottom left) is the zero (0) position of the
function. So, a function value of zero does not mean to have zero
dots but one.
Example of 3 columns and 3 rows (upside down view of 'normal' braille/drawille position)
::
| ^ y axis
| |
| ,_____,,_____,,_____,
7 + | * *|| * *|| * *|
6 + | * *|| * *|| * *|
5 + | * *|| * *|| * *|
4 + | * *|| * *|| * *|
| ;=====;;=====;;=====;
3 + | * *|| * *|| * *|
2 + | * *|| * *|| * *|
1 + | * *|| * *|| * *|
0 + - | * *|| * *|| * *| ---> x axis
| ;=====;;=====;;=====;
-1 + | * *|| * *|| * *|
-2 + | * *|| * *|| * *|
-3 + | * *|| * *|| * *|
-4 + | * *|| * *|| * *|
| `````````````````````
| |
`-----+--+---+--+---+--+-------------
-2 -1 0 1 2 3
.. rubric:: Block Elements Block
* Dots or pixels per character in vertical direction: 8
* Dots or pixels per character in horizontal direction: 1
Example of 3 columns and 3 rows
::
| ^ y axis
| |
| ,_____,,_____,,_____,
15 + | --- || --- || --- |
14 + | --- || --- || --- |
13 + | --- || --- || --- |
12 + | --- || --- || --- |
11 + | --- || --- || --- |
10 + | --- || --- || --- |
9 + | --- || --- || --- |
8 + | --- || --- || --- |
| ;=====;;=====;;=====;
7 + | --- || --- || --- |
6 + | --- || --- || --- |
5 + | --- || --- || --- |
4 + | --- || --- || --- |
3 + | --- || --- || --- |
2 + | --- || --- || --- |
1 + | --- || --- || --- |
0 + - | --- || --- || --- | ---> x axis
| ;=====;;=====;;=====;
-1 + | --- || --- || --- |
-2 + | --- || --- || --- |
-3 + | --- || --- || --- |
-4 + | --- || --- || --- |
-5 + | --- || --- || --- |
-6 + | --- || --- || --- |
-7 + | --- || --- || --- |
-8 + | --- || --- || --- |
| `````````````````````
| |
`------+------+------+---------------
-1 0 1
.. rubric:: References
.. [CUDB] `Unicode Database - Blocks <ftp://ftp.unicode.org/Public/UNIDATA/Blocks.txt>`_
.. seealso::
:meth:`transformation`
"""
# scale function to used chars and dots/pixel in y direction (4 for braille characters): from [0, 1] to [0, chars*4-1]
# negate the function because the y axis is pointing downwards: from [0, 1] to [-1, 0] or from [0, chars*4-1] to [-(chars*4-1), 0]
# and becasue of the negation shift the function by one dot/pixel, otherwise the function would need another char, because the 0 is on an another char on the "positive side".
# devide the time by the number of dots/pixel in x direction (2 for braille characters)
window_factor = w/(x_1-x_0)
frame = ""
if char_set in ["line", "braille"]:
import drawille
pixels_horizontal = 2
pixels_vertical = 4
a = -(h*pixels_vertical-1)
b = 1/pixels_horizontal
f = transformation(f, a, b, 0, -1)
canvas = drawille.Canvas()
# divide step width of the sequence by 2 (double density, 2 dots/pixel per char)
# multiplicate x by 2 (2 dots/pixel per char)
for x_i in seq(x_0*window_factor*pixels_horizontal, x_1*window_factor*pixels_horizontal, 1/density):
y_i = f(x_i, t)
if 0 >= y_i >= -h*pixels_vertical:
canvas.set(x_i, y_i)
#frame = canvas.frame(min_x=a*pixel_per_char, min_y=1, max_x=b*pixel_per_char, max_y=21)
frame = canvas.frame()
elif char_set in ["histogram", "block"]:
from .drawblock import histogram
pixels_horizontal = 1
pixels_vertical = 8
a = h*pixels_vertical-1
f = transformation(f, a, 1, 0, 0)
density = min(density, 1) # density max 1!
x = seq(x_0*window_factor, x_1*window_factor, 1/pixels_horizontal/density)
f = [f(xi, t) for xi in x]
frame = histogram(f, x)
return frame