diff --git a/.gitignore b/.gitignore
index 6c4a347..ba1871c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,4 +6,5 @@ __pycache__/
*.egg-info/
# Sphinx documentation
+docs/build/.doctrees/
docs/build/doctrees/
diff --git a/docs/build/html/.buildinfo b/docs/build/html/.buildinfo
index 8106eaa..2c89628 100644
--- a/docs/build/html/.buildinfo
+++ b/docs/build/html/.buildinfo
@@ -1,4 +1,4 @@
# Sphinx build info version 1
# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
-config: 94375f4299332632f508e2242b7c30d8
+config: aaf67f6f94ce2e6ce3750a4b226f6461
tags: 645f666f9bcd5a90fca523b33c5a78b7
diff --git a/docs/build/html/_modules/data.html b/docs/build/html/_modules/data.html
new file mode 100644
index 0000000..679abdd
--- /dev/null
+++ b/docs/build/html/_modules/data.html
@@ -0,0 +1,179 @@
+
+
+
+
+
+
+ data — pylib 2019.5.19 documentation
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for data
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Read and write data to or from file.
+
+.. module:: data
+ :platform: *nix, Windows
+ :synopsis: Handle data files.
+
+.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
+"""
+from__future__importprint_function
+importpickle
+
+
[docs]defdata_read(file_name,x_column,y_column):
+ """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
+
+ :returns: x and y
+ :rtype: tuple(list, list)
+ """
+ importre
+ file=open(file_name)
+ x=[]
+ y=[]
+ forrowinfile:
+ fields=re.split(r'\s+',row.strip())
+ #print(filds)
+ x.append(float(fields[x_column]))
+ y.append(float(fields[y_column]))
+ file.close()
+ returnx,y
+
+
[docs]defdata_load(file_name,verbose=False):
+ """Load stored program objects from binary file.
+
+ :param file_name: file to load
+ :type file_name: str
+ :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:
+ object_data=pickle.load(input)# one load for every dump is needed to load all the data
+ ifverbose:
+ print('found:')
+ print(object_data)
+ exceptIOError:
+ object_data=None
+ ifverbose:
+ print('no saved datas found')
+ returnobject_data
+
+
[docs]defdata_store(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:
+ pickle.dump(object_data,output,pickle.HIGHEST_PROTOCOL)# every dump needs a load
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/docs/build/html/_modules/date.html b/docs/build/html/_modules/date.html
new file mode 100644
index 0000000..3ab7901
--- /dev/null
+++ b/docs/build/html/_modules/date.html
@@ -0,0 +1,222 @@
+
+
+
+
+
+
+ date — pylib 2019.5.19 documentation
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for date
+#!/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>
+"""
+from__future__importdivision,print_function,unicode_literals
+
+
[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]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
+ """
+ vx,vy=vec
+ returntuple([(x+vx,y+vy)for(x,y)inpts])
+
+
+
[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
+ """
+ ox,oy=origin
+
+ # add first point to the end
+ ifkwargsisnotNoneand"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])
+
+ 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]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)
+
+
+#
+# matplotlib format, return lists for x and y
+#
+
+
[docs]defline(point1,point2,samples=2):
+ """\
+ .. 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
+ :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
+ """
+ p1x,p1y=point1
+ p2x,p2y=point2
+
+ denominator=(p1x-p2x)
+
+ ifsamples>2anddenominator>0:
+ x=np.linspace(p1x,p2x)
+ 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=50):
+ """\
+ :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
+ :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):
+ """\
+ :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)
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Function and approximation.
+
+.. 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))
+ x[i+1,:]=x[i,:]+dxdt*Dt# Approximate solution at next value of x
+ 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))
+ x[i+1,:]=x[i,:]+k_2*Dt# Approximate solution at next value of x
+ 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))
+ x[i+1,:]=x[i,:]+1./6*(k_1+2*k_2+2*k_3+k_4)*Dt# Approximate solution at next value of x
+ ifverbose:
+ print('Numerical integration of ODE using explicit 4th-order method (Runge-Kutta) was successful.')
+ returnx
[docs]deffixed_point_iteration(f,xi,t,max_iterations=1000,tol=1e-9,verbose=False):
+ r"""
+ :param f: the function to iterate :math:`f = \Delta{x}(t)`
+ :type f: function
+ :param xi: initial condition :math:`x_i`
+ :type xi: list
+ :param t: time :math:`t`
+ :type t: 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 (default = 1e-9)
+ :type tol: float
+ :param verbose: print information (default = False)
+ :type verbose: bool
+
+ :returns: :math:`x_{i+1}`
+
+ .. math ::
+ x_{i+1} = x_i + \Delta x
+
+ .. seealso::
+ :meth:`dxdt_Dt` for :math:`\Delta x`
+ """
+ xi=x0
+ forjinrange(max_iterations):# Fixed-point iteration
+ Dx=array(f(xi,t,*p))
+ xi1=x0+Dx# Approximate solution at next value of x
+ residuum=norm(xi1-xi)/norm(xi1)
+ xi=xi1
+ ifresiduum<tol:
+ break
+ iterations=j+1# number beginning with 1 therefore + 1
+ returnxi,iterations
+
+
[docs]defi1n(f,x0,t,*p,max_iterations=1000,tol=1e-9,verbose=False):
+ iterations=zeros((len(t),1))
+ x=zeros((len(t),len(x0)))# Preallocate array
+ x[0,:]=x0# Initial condition gives solution at first t
+ foriinrange(len(t)-1):
+ Dt=t[i+1]-t[i]
+ xi=x[i,:]
+ Dx=dxdt_Dt(f,xi,t[i+1],Dt,*p)
+ x[i+1,:],iterations[i]=fixed_point_iteration(Dx,xi,t,max_iterations,tol,verbose)
+ ifverbose:
+ print('Numerical integration of ODE using implicite first-order method (Euler) was successful.')
+ returnx,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
+ foriinrange(len(t)-1):
+ Dt=t[i+1]-t[i]
+ xi=x[i,:]
+ # 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 therefor using Newton-Raphson method
+ forjinrange(max_iterations):# Fixed-point iteration
+ dxdt=array(f(xi,t[i+1],*p))
+ xi1=x[i,:]+dxdt*Dt# Approximate solution at next value of x
+ residuum=norm(xi1-xi)/norm(xi1)
+ xi=xi1
+ ifresiduum<tol:
+ break
+ iterations[i]=j+1
+ x[i+1,:]=xi
+ 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))
+ #x11 = x[i,:] + dxdt*Dt # Approximate solution at next value of x
+
+ 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,maxIterations=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(maxIterations):# Fixed-point iteration
+ #dxdt = array(f(t[i+1], x1, p))
+ #x11 = x[i,:] + dxdt*Dt # Approximate solution at next value of x
+
+ 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
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)¶
-
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)
-
-
-
-
-
-
-
\ No newline at end of file
diff --git a/docs/build/html/genindex.html b/docs/build/html/genindex.html
index 423ac4b..38c3bb8 100644
--- a/docs/build/html/genindex.html
+++ b/docs/build/html/genindex.html
@@ -47,6 +47,7 @@
| L
| M
| N
+ | O
| P
| R
| S
@@ -88,11 +89,13 @@
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.
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/docs/build/html/objects.inv b/docs/build/html/objects.inv
index ab7137e..710740c 100644
--- a/docs/build/html/objects.inv
+++ b/docs/build/html/objects.inv
@@ -2,10 +2,6 @@
# Project: pylib
# Version:
# The remainder of this file is compressed using zlib.
-xڥMr s
-RVS++ Tʋ,)Zb~%$AxPA^CZyFӝbP@m
-5H;Uȷ^ZCn2E:F[YͤF|73`0Pү.tCN
-vƵ5[p
ڨkӠz0 xu ϋ2XX wa +[0 X
ޝM.|8H>-UK
-%J{+ztjC
-g*lm
m 3jV}Wr
--x 3sPʼ~h h8 њR' fL8)$gū*j? K;?A/ɾIJ_ {SMM'&#gNCjn&c ṯ㗧x3/x*J ~?>! ؠC#gBovfem&էt*Sj̳;tMܧϬ9wsL?
\ No newline at end of file
+xڭ0<h! 40CAQcGD# A&< ddH4}fL;mŬ&"e2uHZ*Dp/!LQQmV3inQpY
d |
f |
g |
- s
+ i |
+ n |
+ o
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.
