[docs]defto_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 (fisrt 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:`pylib.function.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=""
+ ifchar_setin["line","braille"]:
+ 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)
+ forx_iinseq(x_0*window_factor*pixels_horizontal,x_1*window_factor*pixels_horizontal,1/density):
+ canvas.set(x_i,f(x_i,t))
+ #frame = canvas.frame(min_x=a*pixel_per_char, min_y=1, max_x=b*pixel_per_char, max_y=21)
+ frame=canvas.frame()
+ elifchar_setin["histogram","block"]:
+ 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)forxiinx]
+ frame=drawblock.histogram(f,x)
+
+ returnframe
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Read and write data to or from file and manipulate data structures.
+
+:Date: 2019-10-11
+
+.. module:: data
+ :platform: *nix, Windows
+ :synopsis: Handle data files and structures.
+
+.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
+"""
+importmath
+importpickle
+
+
[docs]defread(file_name,x_column,y_column,default=None,verbose=False):
+ """Read ascii data file.
+
+ :param filename: file to read
+ :type filename: str
+ :param x_column: column index for the x data (first column is 0)
+ :type x_column: int
+ :param y_column: column index for the y data (first column is 0)
+ :type y_column: int
+ :param default: return object if data loading fails
+ :type default: object
+ :param verbose: verbose information (default = False)
+ :type verbose: bool
+
+ :returns: x and y
+ :rtype: tuple(list, list)
+ """
+ importre
+
+ x=default
+ y=default
+
+ ifverbose:
+ print('check if data is available')
+ try:
+ file=open(file_name)
+ x=[]
+ y=[]
+ forrowinfile:
+ fields=re.split(r'\s+',row.strip())
+ x.append(float(fields[x_column]))
+ y.append(float(fields[y_column]))
+ file.close()
+ exceptIOError:
+ ifverbose:
+ print('data file not found')
+ returnx,y
[docs]defload(file_name,default=None,verbose=False):
+ """Load stored program objects from binary file.
+
+ :param file_name: file to load
+ :type file_name: str
+ :param default: return object if data loading fails
+ :type default: object
+ :param verbose: verbose information (default = False)
+ :type verbose: bool
+
+ :returns: loaded data
+ :rtype: object
+ """
+ ifverbose:
+ print('check if data is available')
+ try:
+ withopen(file_name,'rb')asinput:
+ # one load for every dump is needed to load all the data
+ object_data=pickle.load(input)
+ ifverbose:
+ print('found:')
+ print(object_data)
+ exceptIOError:
+ object_data=default
+ ifverbose:
+ print('no saved datas found')
+ returnobject_data
+
+
[docs]defstore(file_name,object_data):
+ """Store program objects to binary file.
+
+ :param file_name: file to store
+ :type file_name: str
+ :param object_data: data to store
+ :type object_data: object
+ """
+ withopen(file_name,'wb')asoutput:
+ # every dump needs a load
+ pickle.dump(object_data,output,pickle.HIGHEST_PROTOCOL)
+
+
[docs]deffold_list(lst,n):
+ """Convert one-dimensional kx1 array (list) to two-dimensional mxn
+ array. m = k / n
+
+ :param lst: list to convert
+ :type lst: list
+ :param n: length of the second dimenson
+ :type n: int
+
+ :returns: two-dimensional array (list of lists)
+ :rtype: list
+ """
+ k=len(lst)
+ ifk%n==0:
+ length=int(k/n)
+ return[lst[i*n:i*n+n]foriinrange(length)]
+
+
[docs]defseq(start,stop=None,step=1):
+ r"""Create an arithmetic bounded sequence.
+
+ The sequence is one of the following;
+
+ - empty :math:`\{\}=\emptyset`, if start and stop are the same
+ - degenerate :math:`\{a\}`, if the sequence has only one element.
+ - left-close and right-open :math:`[a, b)`
+
+ :param start: start of the sequence, the lower bound. If only start
+ is given than it is interpreted as stop and start will be 0.
+ :type start: int or float
+ :param stop: stop of sequence, the upper bound.
+ :type stop: int or float
+ :param step: step size, the common difference (constant difference
+ between consecutive terms).
+ :type step: int or float
+ :returns: arithmetic bounded sequence
+ :rtype: list
+ """
+ # example of seq(4, 0, -0.4)
+ # without round:
+ # [4.0, 3.6, 3.2, 2.8, 2.4, 2.0, 1.5999999999999996,
+ # 1.1999999999999997, 0.7999999999999998, 0.3999999999999999]
+ # with round:
+ # [4.0, 3.6, 3.2, 2.8, 2.4, 2.0, 1.6, 1.2, 0.8, 0.4]
+ # example of seq(4, 0, -0.41)
+ # without round:
+ # [4.0, 3.59, 3.18, 2.77, 2.3600000000000003,
+ # 1.9500000000000002, 1.54, 1.1300000000000003,
+ # 0.7200000000000002, 0.31000000000000005]
+ # with round:
+ # [4.0, 3.59, 3.18, 2.77, 2.36, 1.95, 1.54, 1.13, 0.72, 0.31]
+ ifstopisNone:
+ returnseq(0,start,step)
+
+ start_str=str(start)
+ start_exp=0
+ if'.'instart_str:
+ start_exp=len(start_str.split('.')[1])
+
+ step_str=str(step)
+ step_exp=0
+ if'.'instep_str:
+ step_exp=len(step_str.split('.')[1])
+
+ exponent=max(start_exp,step_exp)# no stop because it is an open bound
+
+ n=int(math.ceil((stop-start)/float(step)))
+ lst=[]
+ ifn>0:
+ lst=[round(start+step*i,exponent)foriinrange(n)]
+ returnlst
+
+
[docs]defunique_ending(ids,n=1):
+ """From id list get list with unique ending.
+
+ :param ids: ids
+ :type ids: list
+ :param n: minumum chars or ints
+ :type n: int
+
+ :returns: unique ending of ids
+ :rtype: list
+ """
+ ifidsisnotNone:
+ x=[idi[-n:]foridiinids]
+ iflen(x)>len(set(x)):
+ returnunique_ending(ids,n+1)
+ else:
+ returnx
+
+
[docs]defget_id(ids,uide):
+ """Get full id from unique id ending.
+
+ :param ids: ids
+ :type ids: list
+ :param uide: unique id ending
+ :type uide: str
+
+ :returns: full id
+ :rtype: str or int
+ """
+ # take first element, because we know it is a unique ending
+ return[idiforidiinidsifidi.endswith(uide)][0]
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Calculate spacial dates.
+
+:Date: 2018-01-15
+
+.. module:: date
+ :platform: *nix, Windows
+ :synopsis: Special dates.
+
+.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
+"""
+
+
[docs]defgaußsche_osterformel(year):
+ """Gaußsche Osterformel.
+
+ :param year: the year to calculate the Easter Sunday
+ :type year: int
+
+ :returns: the day of Easter Sunday as a day in march.
+ :rtype: int
+
+ :ivar X: Das Jahr / year
+ :vartype X: int
+ :ivar K(X): Die Säkularzahl
+ :vartype K(X): int
+ :ivar M(X): Die säkulare Mondschaltung
+ :vartype M(X): int
+ :ivar S(K): Die säkulare Sonnenschaltung
+ :vartype S(K): int
+ :ivar A(X): Den Mondparameter
+ :vartype A(X): int
+ :ivar D(A,M): Den Keim für den ersten Vollmond im Frühling
+ :vartype D(A,M): int
+ :ivar R(D,A): Die kalendarische Korrekturgröße
+ :vartype R(D,A): int
+ :ivar OG(D,R): Die Ostergrenze
+ :vartype OG(D,R): int
+ :ivar SZ(X,S): Den ersten Sonntag im März
+ :vartype SZ(X,S): int
+ :ivar OE(OG,SZ): Die Entfernung des Ostersonntags von der Ostergrenze (Osterentfernung in Tagen)
+ :vartype OE(OG,SZ): int
+ :ivar OS(OG,OE): Das Datum des Ostersonntags als Märzdatum (32. März = 1. April usw.)
+ :vartype OS(OG,OE): int
+
+ Algorithmus gilt für den Gregorianischen Kalender.
+
+ source: https://de.wikipedia.org/wiki/Gau%C3%9Fsche_Osterformel
+ """
+ x=year
+ k=x//100
+ m=15+(3*k+3)//4-(8*k+13)//25
+ s=2-(3*k+3)//4
+ a=x%19
+ d=(19*a+m)%30
+ r=(d+a//11)//29
+ og=21+d-r
+ sz=7-(x+x//4+s)%7
+ oe=7-(og-sz)%7
+ os=og+oe
+ returnos
+
+
[docs]defeaster_sunday(year):
+ """Easter Sunday.
+
+ :param year: the year to calculate the Easter Sunday
+ :type year: int
+
+ :returns: the day of Easter Sunday
+ :rtype: datetime.date"""
+ importdatetime
+ march=datetime.date(year,3,1)
+ day=march+datetime.timedelta(days=gaußsche_osterformel(year))
+ returnday
+
+
[docs]defeaster_friday(year):
+ """Easter Friday.
+
+ :param year: the year to calculate the Easter Friday
+ :type year: int
+
+ :returns: the day of Easter Friday
+ :rtype: datetime.date"""
+ importdatetime
+ day=easter_sunday(year)+datetime.timedelta(days=-2)
+ returnday
+
+
[docs]defeaster_monday(year):
+ """Easter Monday.
+
+ :param year: the year to calculate the Easter Monday
+ :type year: int
+
+ :returns: the day of Easter Monday
+ :rtype: datetime.date"""
+ importdatetime
+ day=easter_sunday(year)+datetime.timedelta(days=+1)
+ returnday
+
+
[docs]defascension_of_jesus(year):
+ """Ascension of Jesus.
+
+ :param year: the year to calculate the ascension of Jesus
+ :type year: int
+
+ :returns: the day of ascension of Jesus
+ :rtype: datetime.date"""
+ importdatetime
+ day=easter_sunday(year)+datetime.timedelta(days=+39)
+ returnday
+
+
[docs]defpentecost(year):
+ """Pentecost.
+
+ :param year: the year to calculate the Pentecost
+ :type year: int
+
+ :returns: the day of Pentecost
+ :rtype: datetime.date"""
+ importdatetime
+ day=easter_sunday(year)+datetime.timedelta(days=+49)
+ returnday
[docs]deftransformation(f,scale_vertical=1,scale_horizontal=1,
+ shift_horizontal=0,shift_vertical=0):
+ r"""Transform functions.
+
+ :param f: function or list of functions
+ :type f: function or list
+ :param scale_vertical: "a" scale factor in vertical direction (default
+ = 1)
+ :type height: float
+ :param scale_horizontal: "b" scale factor in horizontal direction
+ (default = 1)
+ :type height: float
+ :param shift_horizontal: "c" shift factor in horizontal direction
+ (default = 0)
+ :type height: float
+ :param shift_vertical: "d" shift factor in vertical direction (default
+ = 0)
+ :type height: float
+
+ :returns: transformed function or list of transformed functions
+ :rtype: function or list
+
+ .. math::
+ y = a \, f(b\,(x-c)) + d
+ """
+ # shorter variables
+ a=scale_vertical
+ b=scale_horizontal
+ c=shift_horizontal
+ d=shift_vertical
+
+ # check if f is a function than put it in a list and return only
+ # the function, not the one element list
+ ifcallable(f):
+ returntransformation(
+ [f],scale_vertical=a,scale_horizontal=b,shift_horizontal=c,shift_vertical=d)[0]
+
+ # otherwise assume list of functions
+ ifnotf:# if f is empty. End of the recursive fucntion
+ return[]
+ return[lambdax,t:a*f[0](b*(x-c),t)+d]+\
+ transformation(
+ f[1:],scale_vertical=a,scale_horizontal=b,shift_horizontal=c,shift_vertical=d)
+
+
[docs]defsine_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
+
+ In general, the function is:
+
+ .. math::
+ y(x,t) = A\sin(kx + 2\pi f t + \varphi) + D \\
+ y(x,t) = A\sin(kx + \omega 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`
+ """
+ ifdegree:
+ phi=math.radians(phi)
+ returnlambdax,t=0:A*math.sin(k*x+2*math.pi*f*t+phi)+D
+
+
[docs]defcosine_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
+
+ .. seealso::
+ :meth:`sine_wave`
+ """
+ ifdegree:
+ phi=phi+90
+ else:
+ phi=phi+math.pi/2
+ returnsine_wave(A=A,k=k,f=f,phi=phi,D=D,degree=degree)
+
+
+#
+# Parametric equations
+# roulette
+#
+
+
[docs]defhypotrochoid(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
+ :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)[:1]
+ >>> x, y, theta_end = hyotrochoid(20, 6, 6)
+
+ .. seealso::
+ :meth:`mathematics.lcm`
+ """
+ x=lambdatheta:(R-r)*math.cos(theta)+d*math.cos((R-r)/r*theta)
+ y=lambdatheta:(R-r)*math.sin(theta)-d*math.sin((R-r)/r*theta)
+ theta_end=2*math.pi*lcm(r,R)/R
+ returnx,y,theta_end
+
+
[docs]defepitrochoid(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
+ :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)[:1]
+ >>> x, y, theta_end = epitrochoid(3, 1, 0.5)
+ """
+ x=lambdatheta:(R+r)*math.cos(theta)-d*math.cos((R+r)/r*theta)
+ y=lambdatheta:(R+r)*math.sin(theta)-d*math.sin((R+r)/r*theta)
+ theta_end=2*math.pi
+ returnx,y,theta_end
+
+
[docs]defto_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 (fisrt 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:`pylib.function.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=""
+ ifchar_setin["line","braille"]:
+ 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)
+ forx_iinseq(x_0*window_factor*pixels_horizontal,x_1*window_factor*pixels_horizontal,1/density):
+ canvas.set(x_i,f(x_i,t))
+ #frame = canvas.frame(min_x=a*pixel_per_char, min_y=1, max_x=b*pixel_per_char, max_y=21)
+ frame=canvas.frame()
+ elifchar_setin["histogram","block"]:
+ 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)forxiinx]
+ frame=drawblock.histogram(f,x)
+
+ returnframe
[docs]defdistance(point1,point2):
+ """Distance between two points (or length of a straight line).
+
+ :param point1: first point (first end point of straight line)
+ :type point1: tuple
+ :param point2: second point (second end point of straight line)
+ :type point2: tuple
+
+ :returns: distance between the two points
+ :rtype: float
+ """
+ returnmath.sqrt((point2[0]-point1[0])**2+(point2[1]-point1[1])**2)
+
+
+
[docs]defangle(point1,point2=None):
+ """Angle of point or between two points.
+
+ :param point1: (first) point
+ :type point1: tuple
+ :param point2: second point (default = None)
+ :type point2: tuple
+
+ :returns: angle of point or between two points
+ :rtype: float
+ """
+ ifpoint2isNone:
+ returnmath.atan2(point1[1],point1[0])
+ returnmath.atan2(point2[1]-point1[1],point2[0]-point1[0])
+
+
+
[docs]deftranslate(vec,*pts):
+ """Translate a point or polygon by a given vector.
+
+ :param vec: translation vector
+ :type vec: tuple
+ :param `*pts`: points to translate
+
+ :returns: (point_x, point_y) or (point1, point2, ...)
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`translate_xy`
+ """
+ vx,vy=vec
+ returntuple([(x+vx,y+vy)for(x,y)inpts])
+
+
+
[docs]deftranslate_xy(vec,x,y):
+ """Translate a point or polygon by a given vector.
+
+ :param vec: translation vector
+ :type vec: tuple
+ :param x: points to translate
+ :type x: int or float or list
+ :param y: points to translate
+ :type y: int or float or list
+
+ :returns: (x', y')
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`translate`
+ """
+ vx,vy=vec
+
+ ifnothasattr(x,"__len__"):
+ x=[x]
+
+ ifnothasattr(y,"__len__"):
+ y=[y]
+
+ xp=[xi+vxforxiinx]
+ yp=[yi+vyforyiiny]
+
+ # no list if it is only one value
+ iflen(xp)==1:
+ returnxp[0],yp[0]
+
+ returnxp,yp
+
+
+
[docs]defrotate(origin,angle,*pts,**kwargs):
+ """Rotate a point or polygon counterclockwise by a given angle
+ around a given origin. The angle should be given in radians.
+
+ :param origin: the center of rotation
+ :type origin: tuple
+ :param angle: the rotation angle
+ :type angle: int or float
+ :param `*pts`: points to rotate
+ :param `**kwargs`: options
+
+ :returns: (point_x, point_y) or (point1, point2, ...)
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`rotate_xy`
+ """
+ ox,oy=origin
+
+ # add first point to the end
+ if"closed"inkwargsandkwargs["closed"]isTrue:
+ pts+=(pts[0],)
+
+ result=tuple([(ox+math.cos(angle)*(px-ox)-math.sin(angle)*(py-oy),
+ oy+math.sin(angle)*(px-ox)+math.cos(angle)*(py-oy))
+ for(px,py)inpts])
+
+ # no tuple in tuple if it is only one point
+ iflen(pts)==1:
+ returnresult[0][0],result[0][1]
+
+ returnresult
+
+
+
[docs]defrotate_deg(origin,angle,*pts,**kwargs):
+ """Rotate a point or polygon counterclockwise by a given angle
+ around a given origin. The angle should be given in degrees.
+
+ :param origin: the center of rotation
+ :type origin: tuple
+ :param angle: the rotation angle
+ :type angle: int or float
+ :param `*pts`: points to rotate
+ :param `**kwargs`: options
+
+ :returns: (point_x, point_y) or (point1, point2, ...)
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`rotate`
+ """
+ returnrotate(origin,angle*math.pi/180,*pts,**kwargs)
+
+
+
[docs]defrotate_xy(origin,angle,x,y,**kwargs):
+ """Rotate x and y coordinates counterclockwise by a given angle
+ around a given origin. The angle should be given in radians.
+
+ :param origin: the center of rotation
+ :type origin: tuple
+ :param angle: the rotation angle
+ :type angle: int or float
+ :param x: x coordinates
+ :type x: int or float or list
+ :param y: y coordinates
+ :type y: int or float or list
+ :param `**kwargs`: options
+
+ .. seealso::
+ :meth:`rotate`
+ """
+ ox,oy=origin
+
+ ifnothasattr(x,"__len__"):
+ x=[x]
+
+ ifnothasattr(y,"__len__"):
+ y=[y]
+
+ # add first point to the end
+ if"closed"inkwargsandkwargs["closed"]isTrue:
+ x.append(x[0])
+ y.append(y[0])
+
+ x_result=[ox+math.cos(angle)*(xi-ox)-math.sin(angle)*(yi-oy)
+ forxi,yiinzip(x,y)]
+
+ y_result=[oy+math.sin(angle)*(xi-ox)+math.cos(angle)*(yi-oy)
+ forxi,yiinzip(x,y)]
+
+ # no list if it is only one value
+ iflen(x_result)==1:
+ returnx_result[0],y_result[0]
+
+ returnx_result,y_result
+
+
+
[docs]defrectangle(width,height):
+ """\
+ :param width: the width of the rectangle
+ :type width: int or float
+ :param height: the height of the rectangle
+ :type height: int or float
+
+ :returns: (point1, point2, point3, point4)
+ :rtype: tuple
+ """
+ pt1=(-width/2,-height/2)
+ pt2=(width/2,-height/2)
+ pt3=(width/2,height/2)
+ pt4=(-width/2,height/2)
+ returnpt1,pt2,pt3,pt4,pt1
+
+
+
[docs]defsquare(width):
+ """\
+ :param width: the edge size of the square
+ :type width: int or float
+
+ :returns: (point1, point2, point3, point4)
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`rectangle`
+ """
+ returnrectangle(width,width)
+
+
+
[docs]deflines(pts,**kwargs):
+ """Lines defined by a list of end points.
+
+ :param pts: list of points in absolute global coordinate system. If
+ keyword inc is given than the inc decides what the left and the
+ right end point of the line is, otherwise it is assumed that the
+ points build a solid line, that is lines between the given points
+ in given order.
+ :type pts: list
+
+ :param `**kwargs`: options:
+
+ * deformation -- list of points. Additional deformation
+ (translation) at point.
+ * factor -- factor of the deformation (default = 1).
+ * inc -- the incidence table, a list of 2 element lists. The inc
+ decides what the left and the right end point of the line is.
+ * index_offset -- starting index of lists (default = 0).
+
+ :returns: list of endpoints for each line;
+ [((point1_x, point1_y), (point2_x, point2_y)),
+ (p1, p2),
+ ...]
+ :rtype: list
+
+ .. seealso::
+ :meth:`~geometry_plot.plot_lines` of the :mod:`geometry_plot`
+ module to plot the lines
+ """
+ if'index_offset'notinkwargs:
+ kwargs['index_offset']=0
+ if'deformation'inkwargs:
+ if'factor'notinkwargs:
+ kwargs['factor']=1
+ pts=[(p[0]+d[0]*kwargs['factor'],p[1]+d[1]*kwargs['factor'])for
+ p,dinzip(pts,kwargs['deformation'])]
+ if'inc'inkwargs:
+ return[(pts[l-kwargs['index_offset']],
+ pts[r-kwargs['index_offset']])forl,rinkwargs['inc']]
+ returnlist(zip(pts[:-1],pts[1:]))
+
+
+
[docs]defcubics(pts,**kwargs):
+ """Cubic lines defined by a list of two end points. The deformation
+ as displacement and rotation (radians) is defined element wise as
+ keyword argument deformation or global node wise as
+ global_deformation. The global coordinate system is xy. x in the
+ right direction and y in the top direction.
+
+ :param pts: list of points in absolute global coordinate system. If
+ keyword inc is given than the inc decides what the left and the
+ right end point of the line is, otherwise it is assumed that the
+ points build a solid line, that is lines between the given points
+ in given order.
+ :type pts_rot: list
+
+ :param `**kwargs`: options:
+
+ * deformation -- list of deformation element wise. Additional
+ deformation (translation and rotation in radians) at element
+ left and right node.
+ * rotation_plane -- rotation plane of the element wise
+ deformation defined by a string; either 'xy' or 'xz' (default
+ = 'xy'). x in the right direction and y in the top direction
+ or z in the bottom direction.
+ * global_deformation -- list of deformation global node wise.
+ Additional deformation (horizontal translation, vertical
+ translation and rotation in radians) at node.
+ * factor -- factor of the derformation (default = 1).
+ * inc -- the incidence table, a list of 2 element lists. The inc
+ decides what the left and the right end point of the line is.
+ * index_offset -- starting index of lists (default = 0).
+
+ :returns: list of endpoints for each line;
+ [(((point1_x, point1_y) angle1), ((point2_x, point2_y), angle2),
+ (p1, angle1, p2, angle2),
+ ...]
+ :rtype: list
+ """
+ if'index_offset'notinkwargs:
+ kwargs['index_offset']=0
+ if'deformation'inkwargsor'global_deformation'inkwargs:
+ if'factor'notinkwargs:
+ kwargs['factor']=1
+ if'inc'inkwargs:
+ if'global_deformation'inkwargs:
+ lr=[(pts[l-kwargs['index_offset']],
+ pts[r-kwargs['index_offset']])forl,rinkwargs['inc']]
+ ang=[angle(l,r)forl,rinlr]
+ # system deformation
+ U=kwargs['global_deformation']
+ if'rotation_plane'inkwargsand'xz'==kwargs['rotation_plane']:
+ # system deformation left X Z RY right X Z RY element wise
+ # Z downwards -> convert to X Y RZ
+ Ue=[[U[(l-kwargs['index_offset'])][0],
+ -U[(l-kwargs['index_offset'])][1],
+ U[(l-kwargs['index_offset'])][2],
+ U[(r-kwargs['index_offset'])][0],
+ -U[(r-kwargs['index_offset'])][1],
+ U[(r-kwargs['index_offset'])][2]]forl,rinkwargs['inc']]
+ else:
+ # system deformation left X Y RZ right X Y RZ element wise
+ Ue=[[U[(l-kwargs['index_offset'])][0],
+ U[(l-kwargs['index_offset'])][1],
+ U[(l-kwargs['index_offset'])][2],
+ U[(r-kwargs['index_offset'])][0],
+ U[(r-kwargs['index_offset'])][1],
+ U[(r-kwargs['index_offset'])][2]]forl,rinkwargs['inc']]
+ # element deformation, X Y RZ to x r rz
+ # back transformation T^T = [[c, s, 0], [-s, c, 0], [0, 0, 1]]
+ u=[[(math.cos(angi)*Uei[0]+math.sin(angi)*Uei[1])*kwargs['factor'],
+ (-math.sin(angi)*Uei[0]+math.cos(angi)*Uei[1])*kwargs['factor'],
+ Uei[2]*kwargs['factor'],
+ (math.cos(angi)*Uei[3]+math.sin(angi)*Uei[4])*kwargs['factor'],
+ (-math.sin(angi)*Uei[3]+math.cos(angi)*Uei[4])*kwargs['factor'],
+ Uei[5]*kwargs['factor']]
+ forUei,angiinzip(Ue,ang)]
+ else:# deformation
+ # the deformation is in element coordinate system, therefore the angle is needed
+ if'rotation_plane'inkwargsand'xz'==kwargs['rotation_plane']:
+ u=[[ue[0]*kwargs['factor'],
+ -ue[1]*kwargs['factor'],
+ ue[2]*kwargs['factor'],
+ ue[3]*kwargs['factor'],
+ -ue[4]*kwargs['factor'],
+ ue[5]*kwargs['factor']]forueinkwargs['deformation']]
+ else:
+ u=[[ui*kwargs['factor']foruiinue]forueinkwargs['deformation']]
+ return[(pts[l-kwargs['index_offset']],
+ pts[r-kwargs['index_offset']],
+ d)for(l,r),dinzip(kwargs['inc'],u)]
+ returnlist(zip(pts[:-1],pts[1:]))
+
+
+
[docs]definterpolate_hermite(lvd,lr,rvd,rr,lhd=0,rhd=0,scale_x=1,scale_y=1,samples=10):
+ r"""Interpolate cubic line with hermite boundary conditions.
+
+ :param lvd: left vertcal deflection
+ :type lvd: int or float
+ :param lr: left rotation
+ :type lr: int or float
+ :param rvd: right vertical deflection
+ :type rvd: int or float
+ :param rr: right rotation
+ :type rr: int or float
+ :param lhd: left horizontal deformation (default = 0)
+ :type lhd: int or float
+ :param rhd: right horizontal deformation (default = 0)
+ :type rhd: int or float
+ :param scale_x: length of element (default = 1)
+ :type scale_x: int or float
+ :param scale_y: factor of the deformation (default = 1).
+ This does not change the length.
+ :type scale_y: int or float
+ :param samples: number of sampling points (default = 10)
+ :type samples: int
+
+ .. math::
+ s = \frac{x - x_1}{L} \\
+ x = s\,L + x_1
+
+ """
+ L=scale_x
+ Lp=L+rhd-lhd
+
+ # x=[0,1] in non-dimensional coordinates
+ x=np.linspace(0,1,num=samples)
+ N1=1-3*x**2+2*x**3
+ N2=(x-2*x**2+x**3)*Lp
+ N3=3*x**2-2*x**3
+ N4=(-x**2+x**3)*Lp
+ x=x*Lp
+
+
+ # x=[0,L] in global coordinates
+ x=np.linspace(0,Lp,num=samples)
+ N1=1-3*x**2/Lp**2+2*x**3/Lp**3
+ N2=x-2*x**2/Lp+x**3/Lp**2
+ N3=3*x**2/Lp**2-2*x**3/Lp**3
+ N4=-x**2/Lp+x**3/Lp**2
+
+
+ v=N1*lvd+N2*lr+N3*rvd+N4*rr
+ x=x+lhd
+ y=v*scale_y
+
+ returnx,y
+
+
+#
+# matplotlib format, return lists for x and y
+#
+
+
[docs]defline(point1,point2,samples=2):
+ """Line defined by two end points.
+
+ .. math::
+ y = \\frac{y_2-y_1}{x_2-x_1}(x-x_1) + y_1
+
+ :param point1: one end point
+ :type point1: tuple
+ :param point2: other end point
+ :type point2: tuple
+ :param samples: number of sampling points (default = 2)
+ :type samples: int
+
+ :returns: ((point1_x, point2_x), (points1_y, point2_y)) or
+ ([sample_point1_x, sample_point2_x, ...],
+ [sample_points1_y, sample_point2_y, ...])
+ :rtype: tuple
+
+ :Example:
+
+ >>> x, y = line((0, 0), (1, 0))
+ >>> print(x, y)
+ ((0, 1), (0, 0))
+ """
+ p1x,p1y=point1
+ p2x,p2y=point2
+
+ denominator=(p1x-p2x)
+
+ ifsamples>2anddenominator>0:
+ x=np.linspace(p1x,p2x,samples)
+ a=(p1y-p2y)/denominator
+ b=(p1x*p2y-p2x*p1y)/denominator
+ y=a*x+b
+ returnx,y
+ return(p1x,p2x),(p1y,p2y)# matplotlib format
+
+
+
[docs]defcubic(point1,angle1,point2,angle2,samples=10):
+ """Cubic line defined by two end points and the rotation in radians
+ at the points.
+
+ :param point1: one end point
+ :type point1: tuple
+ :param angle1: the slope at the one end point
+ :type angle1: int or float
+ :param point2: other end point
+ :type point2: tuple
+ :param angle2: the slope at the other end point
+ :type angle2: int or float
+ :param samples: number of sampling points (default = 10)
+ :type samples: int
+
+ :returns: ([sample_point1_x, sample_point2_x, ...],
+ [sample_points1_y, sample_point2_y, ...])
+ :rtype: tuple
+ """
+ p1x,p1y=point1
+ p2x,p2y=point2
+
+ x=np.linspace(p1x,p2x,num=samples)
+
+ p1ys=math.tan(angle1)
+ p2ys=math.tan(angle2)
+ a=(p1x*p1ys+p1x*p2ys-p2x*p1ys-p2x*p2ys-2*p1y+2*p2y)/(p1x**3-3*p1x**2*p2x+3*p1x*p2x**2-p2x**3)
+ b=(-p1x**2*p1ys-2*p1x**2*p2ys-p1x*p2x*p1ys+p1x*p2x*p2ys+3*p1x*p1y-3*p1x*p2y+2*p2x**2*p1ys+p2x**2*p2ys+3*p2x*p1y-3*p2x*p2y)/(p1x**3-3*p1x**2*p2x+3*p1x*p2x**2-p2x**3)
+ c=(p1x**3*p2ys+2*p1x**2*p2x*p1ys+p1x**2*p2x*p2ys-p1x*p2x**2*p1ys-2*p1x*p2x**2*p2ys-6*p1x*p2x*p1y+6*p1x*p2x*p2y-p2x**3*p1ys)/(p1x**3-3*p1x**2*p2x+3*p1x*p2x**2-p2x**3)
+ d=(-p1x**3*p2x*p2ys+p1x**3*p2y-p1x**2*p2x**2*p1ys+p1x**2*p2x**2*p2ys-3*p1x**2*p2x*p2y+p1x*p2x**3*p1ys+3*p1x*p2x**2*p1y-p2x**3*p1y)/(p1x**3-3*p1x**2*p2x+3*p1x*p2x**2-p2x**3)
+ y=a*x**3+b*x**2+c*x+d
+ returnx,y
+
+
+
[docs]defcubic_deg(point1,angle1,point2,angle2):
+ """Cubic line defined by two end points and the roation in degree
+ at the points.
+
+ :param point1: one end point
+ :type point1: tuple
+ :param angle1: the slope at the one end point
+ :type angle1: int or float
+ :param point2: other end point
+ :type point2: tuple
+ :param angle2: the slope at the other end point
+ :type angle2: int or float
+
+ :returns: ([sample_point1_x, sample_point2_x, ...],
+ [sample_points1_y, sample_point2_y, ...])
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`cubic`
+ """
+ returncubic(point1,angle1*math.pi/180,point2,angle2*math.pi/180)
[docs]defplot_lines(lns,**kwargs):
+ if'color'notinkwargs:
+ # Colors from "tab10" colormap
+ cmap=pylab.get_cmap("tab10")
+ kwargs['color']=cmap(0)
+
+ # Colors from color cycle
+ cycle=pylab.rcParams['axes.prop_cycle'].by_key()['color']
+ kwargs['color']=cycle[0]
+
+ # Colors from CN notation (same as the first 10 colors from the color cycle)
+ kwargs['color']="C0"
+
+ forlninlns:
+ x,y=line(*ln)
+ pylab.plot(x,y,**kwargs)
+
+
+
[docs]defplot_cubic_lines(lns,**kwargs):
+ if'color'notinkwargs:
+ # Colors from "tab10" colormap
+ cmap=pylab.get_cmap("tab10")
+ kwargs['color']=cmap(0)
+
+ # Colors from color cycle
+ cycle=pylab.rcParams['axes.prop_cycle'].by_key()['color']
+ kwargs['color']=cycle[0]
+
+ # Colors from CN notation (same as the first 10 colors from the color cycle)
+ kwargs['color']="C0"
+
+ factor=1
+ if'factor'inkwargs:
+ factor=kwargs['factor']
+
+ samples=10
+ if'samples'inkwargs:
+ samples=kwargs['samples']
+
+ kwargs={k:vfork,vinkwargs.items()ifkinpylab_Line2D_properties}
+
+ forlninlns:
+ l,r,d=ln
+ L=distance(l,r)
+ ang=angle(l,r)
+ x,y=interpolate_hermite(d[1],d[2],d[4],d[5],d[0],d[3],
+ scale_x=L,scale_y=factor,samples=samples)
+ x,y=rotate_xy((0,0),ang,x,y)
+ x,y=translate_xy(l,x,y)
+
+ pylab.plot(x,y,**kwargs)
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Function and approximation.
+
+:Date: 2019-10-15
+
+.. module:: fit
+ :platform: *nix, Windows
+ :synopsis: Function and approximation.
+
+.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
+"""
+from__future__importprint_function
+frompylabimportarray,argmax,gradient,exp,sqrt,log,linspace
+fromscipy.optimizeimportcurve_fit
+
+
[docs]defgauss(x,*p):
+ """Gauss distribution function.
+
+ .. math::
+ f(x)=ae^{-(x-b)^{2}/(2c^{2})}
+
+ :param x: positions where the gauss function will be calculated
+ :type x: int or float or list or numpy.ndarray
+ :param p: gauss parameters [a, b, c, d]:
+
+ * a -- amplitude (:math:`\\int y \\,\\mathrm{d}x=1 \\Leftrightarrow a=1/(c\\sqrt{2\\pi})` )
+ * b -- expected value :math:`\\mu` (position of maximum, default = 0)
+ * c -- standard deviation :math:`\\sigma` (variance :math:`\\sigma^2=c^2`)
+ * d -- vertical offset (default = 0)
+ :type p: list
+
+ :returns: gauss values at given positions x
+ :rtype: numpy.ndarray
+ """
+ x=array(x)# cast e. g. list to numpy array
+ a,b,c,d=p
+ returna*exp(-(x-b)**2./(2.*c**2.))+d
+
+
[docs]defgauss_fit(x,y,e=None,x_fit=None,verbose=False):
+ """Fit Gauss distribution function to data.
+
+ :param x: positions
+ :type x: int or float or list or numpy.ndarray
+ :param y: values
+ :type y: int or float or list or numpy.ndarray
+ :param e: error values (default = None)
+ :type e: int or float or list or numpy.ndarray
+ :param x_fit: positions of fitted function (default = None, if None then x
+ is used)
+ :type x_fit: int or float or list or numpy.ndarray
+ :param verbose: verbose information (default = False)
+ :type verbose: bool
+
+ :returns:
+ * numpy.ndarray -- fitted values (y_fit)
+ * numpy.ndarray -- parameters of gauss distribution function (popt:
+ amplitude a, expected value :math:`\\mu`, standard deviation
+ :math:`\\sigma`, vertical offset d)
+ * numpy.float64 -- full width at half maximum (FWHM)
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`gauss`
+ """
+ x=array(x)# cast e. g. list to numpy array
+ y=array(y)# cast e. g. list to numpy array
+ y_max=max(y)
+ y_max_pos=argmax(y)
+ x_y_max=x[y_max_pos]
+
+ # starting parameter
+ p0=[y_max,x_y_max,.1,y[0]]
+ ifverbose:
+ print('p0:',end=' ')
+ print(p0)
+ ifeisnotNone:
+ popt,pcov=curve_fit(gauss,x,y,p0=p0,sigma=e)
+ else:
+ popt,pcov=curve_fit(gauss,x,y,p0=p0)
+ ifverbose:
+ print('popt:',end=' ')
+ print(popt)
+ #print(pcov)
+
+ FWHM=2*sqrt(2*log(2))*popt[2]
+ ifverbose:
+ print('FWHM',FWHM)
+
+ ifx_fitisNone:
+ x_fit=x
+ y_fit=gauss(x_fit,*popt)
+
+ returny_fit,popt,FWHM
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Numerical solver of ordinary differential equations.
+
+Solves the initial value problem for systems of first order
+ordinary differential equations.
+
+:Date: 2015-09-21
+
+.. module:: ode
+ :platform: *nix, Windows
+ :synopsis: Numerical solver.
+
+.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
+"""
+
+from__future__importdivision,print_function
+fromnumpyimportarray,isnan,sum,zeros,dot
+fromnumpy.linalgimportnorm,inv
+
+
[docs]defe1(f,x0,t,*p,verbose=False):
+ r"""Explicit first-order method /
+ (standard, or forward) Euler method /
+ Runge-Kutta 1st order method.
+
+ de:
+ Euler'sche Polygonzugverfahren / explizite Euler-Verfahren /
+ Euler-Cauchy-Verfahren / Euler-vorwärts-Verfahren
+
+ :param f: the function to solve
+ :type f: function
+ :param x0: initial condition
+ :type x0: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function (thickness, diameter,
+ ...)
+ :param verbose: print information (default = False)
+ :type verbose: bool
+
+ Approximate the solution of the initial value problem
+
+ .. math ::
+ \dot{x} &= f(t,x) \\
+ x(t_0) &= x_0
+
+ Choose a value h for the size of every step and set
+
+ .. math ::
+ t_i = t_0 + i h ~,\quad i=1,2,\ldots,n
+
+ The derivative of the solution is approximated as the forward
+ difference equation
+
+ .. math ::
+ \dot{x}_i = f(t_i, x_i) = \frac{x_{i+1} - x_i}{t_{i+1}-t_i}
+
+ Therefore one step :math:`h` of the Euler method from
+ :math:`t_i` to :math:`t_{i+1}` is
+
+ .. math ::
+ x_{i+1} &= x_i + (t_{i+1}-t_i) f(t_i, x_i) \\
+ x_{i+1} &= x_i + h f(t_i, x_i) \\
+
+ Example 1:
+
+ .. math ::
+ m\ddot{u} + d\dot{u} + ku = f(t) \\
+ \ddot{u} = m^{-1}(f(t) - d\dot{u} - ku) \\
+
+ with
+
+ .. math ::
+ x_1 &= u &\quad \dot{x}_1 = \dot{u} = x_2 \\
+ x_2 &= \dot{u} &\quad \dot{x}_2 = \ddot{u} \\
+
+ becomes
+
+ .. math ::
+ \dot{x}_1 &= x_2 \\
+ \dot{x}_2 &= m^{-1}(f(t) - d x_2 - k x_1) \\
+
+ or
+
+ .. math ::
+ \dot{x} &= f(t,x) \\
+ \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} &=
+ \begin{bmatrix} x_2 \\ m^{-1}(f(t) - d x_2 - k x_1)
+ \end{bmatrix} \\
+ &=
+ \begin{bmatrix} 0 \\ m^{-1} f(t) \end{bmatrix} +
+ \begin{bmatrix} 0 & 1 \\ -m^{-1} k & -m^{-1} d \end{bmatrix}
+ \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}
+
+ Example 2:
+
+ .. math ::
+ m(u)\ddot{u} + d(u,\dot{u})\dot{u} + k(u)u = f(t) \\
+ \ddot{u} = m^{-1}(u)(f(t) - d(u,\dot{u})\dot{u} - k(u)u) \\
+
+ with
+
+ .. math ::
+ x_1 &= u &\quad \dot{x}_1 = \dot{u} = x_2 \\
+ x_2 &= \dot{u} &\quad \dot{x}_2 = \ddot{u} \\
+
+ becomes
+
+ .. math ::
+ \dot{x}_1 &= x_2 \\
+ \dot{x}_2 &=
+ m^{-1}(x_1)(f(t) - d(x_1,x_2) x_2 - k(x_1) x_1) \\
+
+ or
+
+ .. math ::
+ \dot{x} &= f(t,x) \\
+ \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} &=
+ \begin{bmatrix}
+ x_2 \\ m^{-1}(x_1)(f(t) - d(x_1,x_2) x_2 - k(x_1) x_1)
+ \end{bmatrix} \\
+ &=
+ \begin{bmatrix} 0 \\ m^{-1}(x_1) f(t) \end{bmatrix} +
+ \begin{bmatrix}
+ 0 & 1 \\ -m^{-1}(x_1) k(x_1) & -m^{-1} d(x_1,x_2)
+ \end{bmatrix}
+ \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}
+
+ The Euler method is a first-order method, which means that the
+ local error (error per step) is proportional to the square of
+ the step size, and the global error (error at a given time) is
+ proportional to the step size.
+ """
+ x=zeros((len(t),len(x0)))# Preallocate array
+ x[0,:]=x0# Initial condition gives solution at first t
+ foriinrange(len(t)-1):# Calculation loop
+ Dt=t[i+1]-t[i]
+ dxdt=array(f(x[i,:],t[i],*p))
+ # Approximate solution at next value of x
+ x[i+1,:]=x[i,:]+dxdt*Dt
+ ifverbose:
+ print('Numerical integration of ODE using explicit '+
+ 'first-order method (Euler / Runge-Kutta) was successful.')
+ returnx
+
+
[docs]defe2(f,x0,t,*p,verbose=False):
+ r"""Explicit second-order method / Runge-Kutta 2nd order method.
+
+ :param f: the function to solve
+ :type f: function
+ :param x0: initial condition
+ :type x0: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function (thickness, diameter,
+ ...)
+ :param verbose: print information (default = False)
+ :type verbose: bool
+ """
+ x=zeros((len(t),len(x0)))# Preallocate array
+ x[0,:]=x0# Initial condition gives solution at first t
+ foriinrange(len(t)-1):# Calculation loop
+ Dt=t[i+1]-t[i]
+ k_1=array(f(x[i,:],t[i],*p))
+ k_2=array(f(x[i,:]+0.5*Dt*k_1,t[i]+0.5*Dt,*p))
+ # Approximate solution at next value of x
+ x[i+1,:]=x[i,:]+k_2*Dt
+ ifverbose:
+ print('Numerical integration of ODE using explicit '+
+ '2th-order method (Runge-Kutta) was successful.')
+ returnx
+
+
[docs]defe4(f,x0,t,*p,verbose=False):
+ r"""Explicit fourth-order method / Runge-Kutta 4th order method.
+
+ :param f: the function to solve
+ :type f: function
+ :param x0: initial condition
+ :type x0: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function (thickness, diameter,
+ ...)
+ :param verbose: print information (default = False)
+ :type verbose: bool
+ """
+ x=zeros((len(t),len(x0)))# Preallocate array
+ x[0,:]=x0# Initial condition
+ foriinrange(len(t)-1):# Calculation loop
+ Dt=t[i+1]-t[i]
+ k_1=array(f(x[i,:],t[i],*p))
+ k_2=array(f(x[i,:]+0.5*Dt*k_1,t[i]+0.5*Dt,*p))
+ k_3=array(f(x[i,:]+0.5*Dt*k_2,t[i]+0.5*Dt,*p))
+ k_4=array(f(x[i,:]+k_3*Dt,t[i]+Dt,*p))
+ # Approximate solution at next value of x
+ x[i+1,:]=x[i,:]+1./6*(k_1+2*k_2+2*k_3+k_4)*Dt
+ ifverbose:
+ print('Numerical integration of ODE using explicit '+
+ '4th-order method (Runge-Kutta) was successful.')
+ returnx
+
+
[docs]deffpi(f,xi,ti,ti1,*p,max_iterations=1000,tol=1e-9,
+ verbose=False):
+ r"""Fixed-point iteration.
+
+ :param f: the function to iterate :math:`f = \dot{x}(x,t)`
+ :type f: function
+ :param xi: initial condition :math:`x_i`
+ :type xi: list
+ :param ti: time :math:`t_i`
+ :type ti: float
+ :param ti1: time :math:`t_{i+1}`
+ :type ti1: float
+ :param `*p`: parameters of the function (thickness, diameter,
+ ...)
+ :param max_iterations: maximum number of iterations
+ :type max_iterations: int
+ :param tol: tolerance against residuum :math:`\varepsilon`
+ (default = 1e-9)
+ :type tol: float
+ :param verbose: print information (default = False)
+ :type verbose: bool
+
+ :returns: :math:`x_{i}`
+
+ .. math ::
+ x_{i,j=0} = x_{i}
+
+ .. math ::
+ x_{i,j+1} = x_i + \dot{x}(x_{i,j}, t_{i+1})\cdot(t_{i+1}-t_i)
+
+ .. math ::
+ \text{residuum} = \frac{\lVert x_{i,j+1}-x_{i,j}\rVert}
+ {\lVert x_{i,j+1} \rVert} < \varepsilon
+
+ .. math ::
+ x_{i} = x_{i,j=\text{end}}
+ """
+ xij=xi
+ forjinrange(max_iterations):
+ dxdt=array(f(xij,ti1,*p))
+ # Approximate solution at next value of x
+ xij1=xi+dxdt*(ti1-ti)
+ residuum=norm(xij1-xij)/norm(xij1)
+ xij=xij1
+ ifresiduum<tol:
+ break
+ iterations=j+1# number beginning with 1 therefore + 1
+ returnxij,iterations
+
+
[docs]defi1(f,x0,t,*p,max_iterations=1000,tol=1e-9,
+ verbose=False):
+ r"""Implicite first-order method / backward Euler method.
+
+ :param f: the function to solve
+ :type f: function
+ :param x0: initial condition
+ :type x0: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function (thickness, diameter,
+ ...)
+ :param max_iterations: maximum number of iterations
+ :type max_iterations: int
+ :param tol: tolerance against residuum (default = 1e-9)
+ :type tol: float
+ :param verbose: print information (default = False)
+ :type verbose: bool
+
+ The backward Euler method has order one and is A-stable.
+ """
+ iterations=zeros((len(t),1))
+ x=zeros((len(t),len(x0)))# Preallocate array
+ x[0,:]=x0# Initial condition gives solution at first t
+ # x(i+1) = x(i) + f(x(i+1), t(i+1)), exact value of
+ # f(x(i+1), t(i+1)) is not available therefore using
+ # Newton-Raphson method
+ foriinrange(len(t)-1):
+ Dt=t[i+1]-t[i]
+ xi=x[i,:]
+ xi,iteration=fpi(f,xi,t[i],t[i+1],*p,max_iterations,
+ tol,verbose)
+ x[i+1,:]=xi
+ iterations[i]=iteration
+ ifverbose:
+ print('Numerical integration of ODE using implicite '+
+ 'first-order method (Euler) was successful.')
+ returnx,iterations
+
+
[docs]defnewmark_newtonraphson(f,x0,xp0,xpp0,t,*p,gamma=.5,
+ beta=.25,max_iterations=1000,tol=1e-9,verbose=False):
+ r"""Newmark method.
+
+ :param f: the function to solve
+ :type f: function
+ :param x0: initial condition
+ :type x0: list
+ :param xp0: initial condition
+ :type xp0: list
+ :param xpp0: initial condition
+ :type xpp0: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function (thickness, diameter,
+ ...)
+ :param gamma: newmark parameter for velocity (default = 0.5)
+ :type gamma: float
+ :param beta: newmark parameter for displacement (default = 0.25)
+ :type beta: float
+ :param max_iterations: maximum number of iterations
+ :type max_iterations: int
+ :param tol: tolerance against residuum (default = 1e-9)
+ :type tol: float
+ :param verbose: print information (default = False)
+ :type verbose: bool
+ """
+ iterations=zeros((len(t),1))
+ x=zeros((len(t),len(x0)))# Preallocate array
+ xp=zeros((len(t),len(xp0)))# Preallocate array
+ xpp=zeros((len(t),len(xpp0)))# Preallocate array
+ x[0,:]=x0# Initial condition gives solution at first t
+ xp[0,:]=xp0# Initial condition gives solution at first t
+ xpp[0,:]=xpp0# Initial condition gives solution at first t
+ foriinrange(len(t)-1):
+ Dt=t[i+1]-t[i]
+
+ xi=x[i,:].reshape(3,1)
+ xpi=xp[i,:].reshape(3,1)
+ xppi=xpp[i,:].reshape(3,1)
+ x1=xi
+ xp1=xpi
+ xpp1=xppi
+ j=0
+ forjinrange(max_iterations):# Fixed-point iteration
+ #dxdt = array(f(t[i+1], x1, p))
+ # Approximate solution at next value of x
+ #x11 = x[i,:] + dxdt*Dt
+
+ N,dN,dNp,dNpp=f(x1.reshape(-1,).tolist(),
+ xp1.reshape(-1,).tolist(),xpp1.reshape(-1,).tolist(),
+ t[i],*p)
+ ifisnan(sum(dN))orisnan(sum(dNp))orisnan(sum(dNpp)):
+ print('divergiert')
+ break
+
+ xpp11=xpp1-dot(inv(dNpp),(N+dot(dN,(x1-xi))+ \
+ dot(dNp,(xp1-xpi))))
+ xp1=xpi+Dt*((1-gamma)*xppi+gamma*xpp11)
+ x1=xi+Dt*xpi+Dt**2*((.5-beta)*xppi+beta*xpp11)
+
+ residuum=norm(xpp11-xpp1)/norm(xpp11)
+ xpp1=xpp11
+ ifresiduum<tol:
+ break
+ iterations[i]=j+1
+
+ xpp[i+1,:]=xpp1.reshape(-1,).tolist()
+ xp[i+1,:]=xp1.reshape(-1,).tolist()
+ x[i+1,:]=x1.reshape(-1,).tolist()
+ ifverbose:
+ print('Numerical integration of ODE using explicite '+
+ 'newmark method was successful.')
+ returnx,xp,xpp,iterations
+ # x = concatenate((x, xp, xpp), axis=1)
+
+
[docs]defnewmark_newtonraphson_rdk(fnm,x0,xp0,xpp0,t,*p,gamma=.5,
+ beta=.25,max_iterations=1000,tol=1e-9,verbose=False):
+ r"""Newmark method.
+
+ :param f: the function to solve
+ :type f: function
+ :param x0: initial condition
+ :type x0: list
+ :param xp0: initial condition
+ :type xp0: list
+ :param xpp0: initial condition
+ :type xpp0: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function (thickness, diameter,
+ ...)
+ :param gamma: newmark parameter for velocity (default = 0.5)
+ :type gamma: float
+ :param beta: newmark parameter for displacement (default = 0.25)
+ :type beta: float
+ :param max_iterations: maximum number of iterations
+ :type max_iterations: int
+ :param tol: tolerance against residuum (default = 1e-9)
+ :type tol: float
+ :param verbose: print information (default = False)
+ :type verbose: bool
+ """
+ iterations=zeros((len(t),1))
+ x=zeros((len(t),len(x0)))# Preallocate array
+ xp=zeros((len(t),len(xp0)))# Preallocate array
+ xpp=zeros((len(t),len(xpp0)))# Preallocate array
+ x[0,:]=x0# Initial condition gives solution at first t
+ xp[0,:]=xp0# Initial condition gives solution at first t
+ xpp[0,:]=xpp0# Initial condition gives solution at first t
+ foriinrange(len(t)-1):
+ Dt=t[i+1]-t[i]
+
+ rm,rmx,rmxpp,rd,rdx,rdxp,rk,rkx,f=fnm(x[i,:],
+ xp[i,:],xpp[i,:],t[i],*p)
+
+ xi=x[i,:].reshape(3,1)
+ xpi=xp[i,:].reshape(3,1)
+ xppi=xpp[i,:].reshape(3,1)
+ x1=xi
+ xp1=xpi
+ xpp1=xppi
+ j=0
+ forjinrange(max_iterations):# Fixed-point iteration
+ #dxdt = array(f(t[i+1], x1, p))
+ # Approximate solution at next value of x
+ #x11 = x[i,:] + dxdt*Dt
+
+ r=(rmx+rdx+rkx)*Dt**2./4+rdxp*Dt/2+rmxpp
+ rp=f-(rm+dot(rmx,(Dt*xpi+Dt**2./4*xppi))- \
+ dot(rmxpp,xppi)+ \
+ rd+dot(rdx,(Dt*xpi+Dt**2./4*xppi))+ \
+ dot(rdxp,Dt/2*xppi)+ \
+ rk+dot(rkx,(Dt*xpi+Dt**2./4*xppi)))
+ xpp11=dot(inv(r),rp)
+ xp1=xpi+Dt*((1-gamma)*xppi+gamma*xpp11)
+ x1=xi+Dt*xpi+Dt**2*((.5-beta)*xppi+beta*xpp11)
+
+ residuum=norm(xpp11-xpp1)/norm(xpp11)
+ xpp1=xpp11
+ ifresiduum<tol:
+ break
+ iterations[i]=j+1
+
+ xpp[i+1,:]=xpp1.reshape(-1,).tolist()
+ xp[i+1,:]=xp1.reshape(-1,).tolist()
+ x[i+1,:]=x1.reshape(-1,).tolist()
+ ifverbose:
+ print('Numerical integration of ODE using explicite '+
+ 'newmark method was successful.')
+ returnx,xp,xpp,iterations
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Mathmatical models governed by ordinary differential equations.
+
+Describes initial value problems as systems of first order ordinary differential
+equations.
+
+:Date: 2019-05-25
+
+.. module:: ode_model
+ :platform: *nix, Windows
+ :synopsis: Models of ordinary differential equations.
+
+.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
+"""
+from__future__importdivision,print_function
+fromnumpyimportarray,cos,sin,dot,square
+fromnumpy.linalgimportinv
+
+
[docs]defdisk(x,t,*p):
+ """Rotation of an eccentric disk.
+
+ :param x: values of the function
+ :type x: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function
+
+ * diameter
+ * eccentricity
+ * torque
+ """
+ qp1=x[3]
+ qp2=x[4]
+ qp3=x[5]
+ M=array([[1,0,cos(x[2])],[0,1,-sin(x[2])],[0,0,1]])
+ y=array([[-2*p[0]*x[3]+sin(x[2])*x[5]**2-2*p[0]*cos(x[2])*x[5]-x[0]], \
+ [-2*p[0]*x[4]+cos(x[2])*x[5]**2-2*p[0]*sin(x[2])*x[5]-x[1]], \
+ [p[2]-p[1]*x[1]*sin(x[2])+p[1]*x[0]*cos(x[2])]])
+ qp46=dot(inv(M),y)
+ qp4,qp5,qp6=qp46.reshape(-1,).tolist()# 2d array to 1d array to list
+ returnqp1,qp2,qp3,qp4,qp5,qp6
+
+
[docs]defdisk_nm(xn,xpn,xppn,t,*p):
+ """Rotation of an eccentric disk.
+
+ :param xn: values of the function
+ :type xn: list
+ :param xpn: first derivative values of the function
+ :type xpn: list
+ :param xppn: second derivative values of the function
+ :type xppn: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function
+
+ * diameter
+ * eccentricity
+ * torque
+ """
+ N=array([[xppn[0]+cos(xn[2])*xppn[2]+2*p[0]*xpn[0]+2*p[0]*cos(xn[2])*xpn[2]-sin(xn[2])*square(xpn[2])+xn[0]],
+ [xppn[1]-sin(xn[2])*xppn[2]+2*p[0]*xpn[1]-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])+xn[1]],
+ [xppn[2]+p[1]*(-cos(xn[2])*xn[0]+sin(xn[2])*xn[1])-p[2]]])
+ dN=array([[1,0,-sin(xn[2]*xppn[2])-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])],
+ [0,1,-cos(xn[2]*xppn[2])-2*p[0]*cos(xn[2])*xpn[2]+sin(xn[2])*square(xpn[2])],
+ [-p[1]*cos(xn[2]),p[1]*cos(xn[2]),p[1]*(sin(xn[2])*xn[0]+cos(xn[2])*xn[1])]])
+ dNp=array([[2*p[0],0,2*p[0]*cos(xn[2])-2*sin(xn[2])*xpn[2]],
+ [0,2*p[0],-2*p[0]*sin(xn[2])-2*cos(xn[2])*xpn[2]],
+ [0,0,0]])
+ dNpp=array([[1,0,cos(xn[2])],
+ [0,1,-sin(xn[2])],
+ [0,0,1]])
+ returnN,dN,dNp,dNpp
+
+
[docs]defdisk_nmmdk(xn,xpn,xppn,t,*p):
+ """Rotation of an eccentric disk.
+
+ :param xn: values of the function
+ :type xn: list
+ :param xpn: derivative values of the function
+ :type xpn: list
+ :param xppn: second derivative values of the function
+ :type xppn: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function
+
+ * diameter
+ * eccentricity
+ * torque
+ """
+ rm=array([[xppn[0]+cos(xn[2])*xppn[2]],
+ [xppn[1]-sin(xn[2])*xppn[2]],
+ [xppn[2]]])
+ rmx=array([[0,0,-sin(xn[2]*xppn[2])],
+ [0,0,-cos(xn[2]*xppn[2])],
+ [0,0,0]])
+ rmxpp=array([[1,0,cos(xn[2])],
+ [0,1,-sin(xn[2])],
+ [0,0,1]])
+ rd=array([[2*p[0]*xpn[0]+2*p[0]*cos(xn[2])*xpn[2]-sin(xn[2])*square(xpn[2])],
+ [2*p[0]*xpn[1]-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])],
+ [0]])
+ rdx=array([[0,0,-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])],
+ [0,0,-2*p[0]*cos(xn[2])*xpn[2]+sin(xn[2])*square(xpn[2])],
+ [0,0,0]])
+ rdxp=array([[2*p[0],0,2*p[0]*cos(xn[2])-2*sin(xn[2])*xpn[2]],
+ [0,2*p[0],-2*p[0]*sin(xn[2])-2*cos(xn[2])*xpn[2]],
+ [0,0,0]])
+ rk=array([[xn[0]],
+ [xn[1]],
+ [p[1]*(-cos(xn[2])*xn[0]+sin(xn[2])*xn[1])]])
+ rkx=array([[1,0,0],
+ [0,1,0],
+ [-p[1]*cos(xn[2]),p[1]*cos(xn[2]),p[1]*(sin(xn[2])*xn[0]+cos(xn[2])*xn[1])]])
+ f=array([[0],[0],[p[2]]])
+ returnrm,rmx,rmxpp,rd,rdx,rdxp,rk,rkx,f
[docs]defin_seconds(time):
+ """If time is `time.struct_time` convert to float seconds.
+
+ :param time: the time in seconds
+ :type time: float or `time.struct_time`
+
+ :returns: the time in seconds
+ :rtype: float
+ """
+ ifisinstance(time,struct_time):
+ time=mktime(time)
+ returntime
+
+
[docs]defseconds(time):
+ """The seconds of the time.
+
+ :param time: the time in seconds
+ :type time: float or `time.struct_time`
+
+ :returns: seconds, range [0, 60]
+ :rtype: float
+ """
+ returnin_seconds(time)%60
+
+
[docs]defseconds_norm(time):
+ """The seconds normalized to 60 seconds.
+
+ :param time: the time in seconds
+ :type time: float or `time.struct_time`
+
+ :returns: the normalized seconds, range [0, 1]
+ :rtype: float
+ """
+ returnseconds(time)/60
+
+
[docs]defminutes(time):
+ """The minutes of the time.
+
+ :param time: the time in seconds
+ :type time: float or `time.struct_time`
+
+ :returns: minutes, range [0, 60]
+ :rtype: float
+ """
+ returnin_seconds(time)/60%60
+
+
[docs]defminutes_norm(time):
+ """The minutes normalized to 60 minutes.
+
+ :param time: the time in seconds
+ :type time: float or `time.struct_time`
+
+ :returns: the normalized minutes, range [0, 1]
+ :rtype: float
+ """
+ returnminutes(time)/60
+
+
[docs]defhours(time):
+ """The hours of the time.
+
+ :param time: the time in seconds
+ :type time: float or `time.struct_time`
+
+ :returns: hours, range [0, 24]
+ :rtype: float
+ """
+ returnin_seconds(time)/60/60%24
+
+
[docs]defhours_norm(time):
+ """The hours normalized to 24 hours.
+
+ :param time: the time in seconds
+ :type time: float or `time.struct_time`
+
+ :returns: the normalized hours, range [0, 1]
+ :rtype: float
+ """
+ returnhours(time)/24
+
+
[docs]defdays(time):
+ """The days of the time (year).
+
+ :param time: the time in seconds
+ :type time: float or `time.struct_time`
+
+ :returns: hours, range [0, 365.2425]
+ :rtype: float
+ """
+ returnin_seconds(time)/60/60/24%365.2425
+
+
[docs]defdays_norm(time):
+ """The days normalized to 365.2425 (Gregorian, on average) days.
+
+ :param time: the time in seconds
+ :type time: float or `time.struct_time`
+
+ :returns: the normalized days, range [0, 1]
+ :rtype: float
+ """
+ returndays(time)/365.2425
+
+
[docs]deftransform(time_norm,length,offset=0):
+ """Transform normalized time value to new length.
+
+ :param position_norm: the normalized time value to transform
+ :type position_norm: float
+ :param length: the transformation
+ :type length: float
+ :param offset: the offset (default = 0)
+ :type offset: float
+
+ :returns: the transformation value
+ :rtype: float
+ """
+ returntime_norm*length+offset
Represent functions as string frame with a specific character set.
+which are normed to the range of [0, 1] to
+
+
Parameters
+
+
f (function or list) – function or list of functions normed to the range of [0, 1]
+
h (int) – number of chars in vertical direction
+
w (int) – number of chars in horizontal direction
+
char_set (str) – either “braille” or “block”. “braille” uses Unicode
+Characters in the Braille Patterns Block (fisrt 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
+
+
+
+
+
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”
+
+
+
+
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)
A cosine wave is said to be sinusoidal, because,
+\(\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.
+
+
Parameters
+
+
A (float or int) – amplitude
+
k (float or int) – (angular) wave number
+
f (float or int) – ordinary frequency
+
phi (float or int) – phase
+
D (float or int) – non-zero center amplitude
+
degree (bool) – boolean to switch between radians and degree. If
+False phi is interpreted in radians and if True then phi is
+interpreted in degrees.
+
+
+
Results
+
sine wave function of spatial variable x and optional
+time t
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.
+
+
Parameters
+
+
R (float) – radius of the fixed interior circle
+
r – radius of the rolling circle outside of the fixed circle
+
d (float) – distance from the center of the exterior circle
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.
+
+
Parameters
+
+
R (float) – radius of the fixed exterior circle
+
r – radius of the rolling circle inside of the fixed circle
+
d (float) – distance from the center of the interior circle
A sine wave or sinusoid is a mathematical curve that describes a
+smooth periodic oscillation.
+
+
Parameters
+
+
A (float or int) – amplitude
+
k (float or int) – (angular) wave number
+
f (float or int) – ordinary frequency
+
phi (float or int) – phase
+
D (float or int) – non-zero center amplitude
+
degree (bool) – boolean to switch between radians and degree. If
+False phi is interpreted in radians and if True then phi is
+interpreted in degrees.
+
+
+
Results
+
sine wave function of spatial variable x and optional
+time t
+
+
+
In general, the function is:
+
+\[\begin{split}y(x,t) = A\sin(kx + 2\pi f t + \varphi) + D \\
+y(x,t) = A\sin(kx + \omega t + \varphi) + D\end{split}\]
+
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:
+
+\[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.
Represent functions as string frame with a specific character set.
+which are normed to the range of [0, 1] to
+
+
Parameters
+
+
f (function or list) – function or list of functions normed to the range of [0, 1]
+
h (int) – number of chars in vertical direction
+
w (int) – number of chars in horizontal direction
+
char_set (str) – either “braille” or “block”. “braille” uses Unicode
+Characters in the Braille Patterns Block (fisrt 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
+
+
+
+
+
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”
+
+
+
+
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)
Cubic lines defined by a list of two end points. The deformation
+as displacement and rotation (radians) is defined element wise as
+keyword argument deformation or global node wise as
+global_deformation. The global coordinate system is xy. x in the
+right direction and y in the top direction.
+
+
Parameters
+
+
pts – list of points in absolute global coordinate system. If
+keyword inc is given than the inc decides what the left and the
+right end point of the line is, otherwise it is assumed that the
+points build a solid line, that is lines between the given points
+in given order.
+
**kwargs –
options:
+
+
deformation – list of deformation element wise. Additional
+deformation (translation and rotation in radians) at element
+left and right node.
+
rotation_plane – rotation plane of the element wise
+deformation defined by a string; either ‘xy’ or ‘xz’ (default
+= ‘xy’). x in the right direction and y in the top direction
+or z in the bottom direction.
+
global_deformation – list of deformation global node wise.
+Additional deformation (horizontal translation, vertical
+translation and rotation in radians) at node.
+
factor – factor of the derformation (default = 1).
+
inc – the incidence table, a list of 2 element lists. The inc
+decides what the left and the right end point of the line is.
+
index_offset – starting index of lists (default = 0).
+
+
+
+
+
Returns
+
list of endpoints for each line;
+[(((point1_x, point1_y) angle1), ((point2_x, point2_y), angle2),
+(p1, angle1, p2, angle2),
+…]
pts (list) – list of points in absolute global coordinate system. If
+keyword inc is given than the inc decides what the left and the
+right end point of the line is, otherwise it is assumed that the
+points build a solid line, that is lines between the given points
+in given order.
+
**kwargs –
options:
+
+
deformation – list of points. Additional deformation
+(translation) at point.
+
factor – factor of the deformation (default = 1).
+
inc – the incidence table, a list of 2 element lists. The inc
+decides what the left and the right end point of the line is.
+
index_offset – starting index of lists (default = 0).
+
+
+
+
+
Returns
+
list of endpoints for each line;
+[((point1_x, point1_y), (point2_x, point2_y)),
+(p1, p2),
+…]
x (int or float or list or numpy.ndarray) – positions where the gauss function will be calculated
+
p (list) –
gauss parameters [a, b, c, d]:
+
+
a – amplitude (\(\int y \,\mathrm{d}x=1 \Leftrightarrow a=1/(c\sqrt{2\pi})\) )
+
b – expected value \(\mu\) (position of maximum, default = 0)
+
c – standard deviation \(\sigma\) (variance \(\sigma^2=c^2\))
+
d – vertical offset (default = 0)
+
+
+
+
+
Returns
+
gauss values at given positions x
+
+
Return type
+
numpy.ndarray
+
+
+
+
+
+
+gauss_fit(x, y, e=None, x_fit=None, verbose=False)[source]¶
+
Fit Gauss distribution function to data.
+
+
Parameters
+
+
x (int or float or list or numpy.ndarray) – positions
+
y (int or float or list or numpy.ndarray) – values
+
e (int or float or list or numpy.ndarray) – error values (default = None)
+
x_fit (int or float or list or numpy.ndarray) – positions of fitted function (default = None, if None then x
+is used)
+
verbose (bool) – verbose information (default = False)
+
+
+
Returns
+
+
numpy.ndarray – fitted values (y_fit)
+
numpy.ndarray – parameters of gauss distribution function (popt:
+amplitude a, expected value \(\mu\), standard deviation
+\(\sigma\), vertical offset d)
Integration of \(f(x)\) using the trapezoidal rule
+(Simpson’s rule, Kepler’s rule).
+
de: Trapezregel, Simpsonregel (Thomas Simpson), Keplersche
+Fassregel (Johannes Kepler)
+
+
Parameters
+
+
f (function or list) – function to integrate.
+
a (float) – lower limit of integration (default = 0).
+
b (float) – upper limit of integration (default = 1).
+
N (int) – specify the number of subintervals.
+
x (list) – variable of integration, necessary if f is a list
+(default = None).
+
verbose (bool) – print information (default = False)
+
+
+
Returns
+
the definite integral as approximated by trapezoidal
+rule.
+
+
Return type
+
float
+
+
+
The trapezoidal rule approximates the integral by the area of a
+trapezoid with base h=b-a and sides equal to the values of the
+integrand at the two end points.
The Euler method is a first-order method, which means that the
+local error (error per step) is proportional to the square of
+the step size, and the global error (error at a given time) is
+proportional to the step size.