Benutzer:Alva2004/Programme/Numpyfunc
Zur Navigation springen
Zur Suche springen
import bisect, numpy, math, os, random, sys, time, platform
if sys.version_info[ 0 ] == 3: # python3
from collections import UserDict
else:
from UserDict import UserDict
#import scipy.fftpack, scipy.optimize
voigt = ( ( 0,0 ), ( 1,1 ), ( 2,2 ), ( 1,2 ), ( 0,2 ), ( 0,1 ) )
#==============================================================================
def test( text, soll, ist, tol=0, fout=1 ):
'''ausgabe:
fout=0: keine
fout=1: text und ok und bei fehler vergleich
fout=2: text und ok
fout=3: nur bei fehler text und ok
'''
if tol == 0:
ok = numpy.allclose( ist, soll )
else:
ok = numpy.allclose( ist, soll, atol=tol, rtol=tol )
if fout:
if not ok or fout < 3:
print( text.ljust( 8, ' ' ), ok )
if not ok and fout == 1:
print( 'ist '.ljust( 8, ' ' ), ist )
print( 'soll'.ljust( 8, ' ' ), soll )
return ok
#==============================================================================
def cumtest( cum, text, soll, ist, tol=0, fout=1, ferr=0 ):
'''cumulierter test, cum = ( ok, anz ), ok=anzahl oks, anz=anzahl tests
ausgabe:
fout=0: keine
fout=1: ok und bei fehler soll-ist-vergleich
fout=2: nur ok
fout=3: nur bei fehler ok
'''
#print( 'cumtest fout=', fout )
( ok0, anz ) = cum
if tol == 0:
ok = numpy.allclose( ist, soll )
else:
ok = numpy.allclose( ist, soll, atol=tol, rtol=tol )
if fout == 1: # ok und falls nicht ok soll-ist-vergleich
print( text.ljust( 8, ' ' ), ok )
if not ok: # vergleich
print( 'ist '.ljust( 8, ' ' ), ist )
print( 'soll'.ljust( 8, ' ' ), soll )
elif fout == 2: # ok jedenfalls
print( text.ljust( 8, ' ' ), ok )
elif fout == 3 and not ok: # ok nur falls nicht ok
print( text.ljust( 8, ' ' ), ok )
cum[ : ] = ( ok0 + ok, anz + 1 )
if not ok and ferr:
raise ValueError( 'cumtest %s falsch'%text )
return ok
#==============================================================================
def testn( text, *par, **dpar ):
dtol = { 'tol': 0, 'atol': 0, 'rtol': 0 }
dtol.update( dpar )
tol = dtol[ 'tol' ]
par = list( par )
par.reverse()
p0 = par.pop()
if tol == 0:
ok = 1
while par:
ok = ok and numpy.allclose( p0, par.pop() )
else:
atol = rtol = tol
if dtol[ 'atol' ] != 0:
atol = dtol[ 'atol' ]
if dtol[ 'rtol' ] != 0:
rtol = dtol[ 'rtol' ]
ok = 1
while par:
ok = ok and numpy.allclose( p0, par.pop(), atol=atol, rtol=rtol )
print( text.ljust( 8, ' ' ), ok )
if not ok:
print( par )
return ok
#==============================================================================
def testup( text, all, fout=0 ):
ok = 1
n = 0
for tup in all:
soll = tup[ 0 ]
m = 1
l = len( tup ) - 1
for ist in tup[ 1: ]:
okn = numpy.allclose( ist, soll )
if okn:
if fout > 1:
print( '%s %d %d/%d ok'%( text, n, m, l ) )
else:
print( '%s %d %d/%d ist '%( text, n, m, l ), ist )
print( '%s %d %d/%d soll'%( text, n, m, l ), soll )
ok = ok and okn
m += 1
n += 1
if fout:
print( text.ljust( 8, ' ' ), ok )
return ok
#==============================================================================
def performance( fkt, par, n=1000, name='' ):
if name == '':
name = str( fkt )
t = time.time()
for i in range( n ):
fkt( *par )
t = time.time() - t
print( 'Zeitverbrauch von %s fuer %d Aufrufe ist %s'%(
name, n, t ) )
#==============================================================================
def scale( a, b, n, meth=0 ):
'''n Punkte zwischen a und b einschliesslich
logarithmisch zunehmend
x = [ a, a + q, a + q**2, ..., a + q**n ]
xn = b = a + q**n
b - a = q**n
log( b - a ) = n*log( q )
logarithmisch abnehmend
x = [ b - q**n, b - q**2, b - q, b ]
a = b - q**n
b - a = q**n
log( b - a ) = n*log( q )
linear zunehmend
x = [ a, a+q, a+3q, a+6q, a+10q, a+15q, .., a+nq ]
= [ a + (i+1)*i/2*q, i=0..n ]
xn = b = a + (n+1)*n/2*q
q = 2*( b - a )/( n*( n + 1 ) )
linear abnehmend
x = [ b-ai*q, .., b-15q, b-10q, b-6q, b-3q, b-q, b ]
= [ b - 0.5*( i + 1 )*i*q, i=0..n ]
a = b - 0.5*( n + 1 )*n*q
0.5*( n + 1 )*n*q = b - a
q = 2.*( b - a )/( ( n + 1 )*n )
'''
if meth == 1: # abstand linear zunehmend
q = 2.*( b - a )/float( n*( n + 1 ) )
ans = [ a + 0.5*(i+1)*i*q for i in range( 1, n ) ]
elif meth == 2: # abstand linear abnehmend
q = 2.*( b - a )/float( n*( n + 1 ) )
ans = [ b - 0.5*(i+1)*i*q for i in range( n-1, 0, -1 ) ]
elif meth == 3: # abstand geometrisch zunehmend
q = math.exp( math.log( b - a )/float( n - 1 ) )
ans = [ a+q**i for i in range( 1, n-1 ) ]
elif meth == 4: # abstand geometrisch abnehmend
q = math.exp( math.log( b - a )/float( n - 1 ) )
ans = [ b-q**i for i in range( n-2, 0, -1 ) ]
else: # konstanter abstand
q = float( b - a )/float( n - 1 )
ans = [ a + q*i for i in range( 1, n-1 ) ]
return [ a ] + ans + [ b ]
#==============================================================================
def fourier( t, x ):
'''fft des Signals x zu Zeiten t.
Gibt Frequenzen ( nicht Kreis-Frequenzen ) und komplexe Werte'''
f = scipy.fftpack.fft( x )
anz = len( t )
dt = t[ -1 ] - t[ 0 ]
fak = 2./anz
freq = [ i/dt for i in range( anz ) ]
a = 2./anz * f
return [ freq, a ]
#==============================================================================
def bits( n, basis=2, dim=3 ):
b = []
if n == 0:
b = [ 0 ]
else:
while n > 0:
b.append( n%basis )
n /= basis
b.reverse()
while len( b ) < dim:
b.insert( 0, 0 )
return b
#==============================================================================
def stib( tup, basis=2 ):
'Umkehrfunktion von bits'
n = 0
z = 1
for i in range( len( tup ) - 1, -1, -1 ):
n += tup[ i ]*z
z *= basis
return n
#==============================================================================
def bits2zahl( tup ):
ok = 0
s = ''
for i in tup:
if ok:
s += str( i )
elif i != 0:
s += str( i )
ok = 1
if not ok: # 0
s = '0'
return int( s )
#==============================================================================
class LeviCivita:
'Levi-Civita Tensor'
#--------------------------------------------------------------------------
def dot( self, v ):
'''gibt E3*v = -X(v) = -v x I
(eijk ei o ej o ek).vl*el = eijk*vk*ei o ej = ei.(ej x ek)*vk*ei o ej
= (ej x v) o ej = -v x ej o ej = -v x 1
'''
return -v.cross( Eye() )
#--------------------------------------------------------------------------
def mmul( self, T ):
'''gibt E3*T
(eijk ei o (ej o ek)):Tmn*em o en = djm*dkn*eijk*Tmn*ei
= ei.(ejxek)*Tjk*ei = (ejxek)*Tjk = i(T)'''
return T.i()
#--------------------------------------------------------------------------
#==============================================================================
class Quaternion( numpy.ndarray ):
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0., 0. ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data, dtype=dtype, subok=1 )
).view( cls )
#--------------------------------------------------------------------------
def __array_finalize__(self, obj):
pass
#--------------------------------------------------------------------------
def real( x ):
'Produkt'
return x[ 0 ]
#--------------------------------------------------------------------------
def imag( x ):
'Produkt'
return Quaternion( [ 0., x[ 1 ], x[ 2 ], x[ 3 ] ] )
#--------------------------------------------------------------------------
def dot( x, y ):
'Skalarprodukt'
return x[ 0 ]*y[ 0 ] + x[ 1 ]*y[ 1 ] + x[ 2 ]*y[ 2 ] + x[ 3 ]*y[ 3 ]
#--------------------------------------------------------------------------
def __neg__( x ):
'Produkt'
return Quaternion( [ -x[ 0 ], -x[ 1 ], -x[ 2 ], -x[ 3 ] ] )
#--------------------------------------------------------------------------
def __abs__( x ):
'Produkt'
return math.sqrt( x.norm() )
#--------------------------------------------------------------------------
def norm( x ):
'Produkt'
return x.dot( x )
#--------------------------------------------------------------------------
def normalize( x ):
'Einheitsquaternion'
return ( 1./abs( x ) )*x
#--------------------------------------------------------------------------
def __rmul__( x, y ):
'Produkt mit einem Skalar'
return Quaternion( [ y*x[ 0 ], y*x[ 1 ], y*x[ 2 ], y*x[ 3 ] ] )
#--------------------------------------------------------------------------
def __rdiv__( x, y ):
'Division: ans = y/x'
return ( y/x.norm() )*x.conj()
#--------------------------------------------------------------------------
def conj( x ):
'Produkt'
return Quaternion( [ x[ 0 ], -x[ 1 ], -x[ 2 ], -x[ 3 ] ] )
#--------------------------------------------------------------------------
def __mul__( x, y ):
'Produkt'
return Quaternion( [
x[ 0 ]*y[ 0 ] - x[ 1 ]*y[ 1 ] - x[ 2 ]*y[ 2 ] - x[ 3 ]*y[ 3 ],
x[ 0 ]*y[ 1 ] + x[ 1 ]*y[ 0 ] - x[ 2 ]*y[ 3 ] + x[ 3 ]*y[ 2 ],
x[ 0 ]*y[ 2 ] + x[ 2 ]*y[ 0 ] - x[ 3 ]*y[ 1 ] + x[ 1 ]*y[ 3 ],
x[ 0 ]*y[ 3 ] + x[ 3 ]*y[ 0 ] - x[ 1 ]*y[ 2 ] + x[ 2 ]*y[ 1 ] ] )
#--------------------------------------------------------------------------
def __div__( x, y ):
'Division ans*y = x'
return ( 1./y.norm() )*x*y.conj()
#--------------------------------------------------------------------------
def cross( x, y ):
'Produkt'
return 0.5*( y*x - x*y )
#--------------------------------------------------------------------------
def drehmatrix( x ):
( a, b, c, d ) = x
return Tensor( [
[ a**2+b**2-c**2-d**2 , 2.*(b*c-a*d) , 2.*(b*d + a*c) ],
[ 2.*(b*c+a*d) , a**2+c**2-b**2-d**2 , 2.*(c*d - a*b) ],
[ 2.*(b*d-a*c) , 2.*(c*d+a*b) , a**2+d**2-b**2-c**2 ] ] )
#--------------------------------------------------------------------------
#==============================================================================
class Signal( numpy.ndarray ):
'Analyse-Werkzeuge fuer ein reelles Signal'
#--------------------------------------------------------------------------
def __new__( cls, data ):
return numpy.asarray( data ).view( cls )
#--------------------------------------------------------------------------
def wert( self, t=0, achse=[] ):
'gibt den Wert bei t auf der achse'
if achse == []: # t ist ein Index
achse = range( len( self ) )
i = bisect.bisect( achse, t )
if i == len( achse ):
ans = self[ -1 ]
else:
t1 = achse[ i ]
x1 = self[ i ]
if t == t1:
ans = x1
else:
t0 = achse[ i - 1 ]
x0 = self[ i - 1 ]
ans = x0 + ( t - t0 )*( x1 - x0 )/( t1 - t0 )
return ans
#--------------------------------------------------------------------------
def abl( self ):
'gibt die Differenzen self[ 2: ] - self[ :-2 ]'
ans = ( [ self[ 1 ] - self[ 0 ] ]
+ list( 0.5*( self[ 2: ] - self[ :-2 ] ) )
+ [ self[ -1 ] - self[ -2 ] ] )
return Signal( ans )
#--------------------------------------------------------------------------
def extrema( self ):
'''extrahiert Minima und Maxima des Real-Teils dieses Signals'''
a = self.real
( a0, a1 ) = a[ :2 ]
imin = []
imax = []
for i in range( 2, len( a ) ):
a2 = a[ i ]
if a1 >= a0 and a1 > a2:
imax.append( i - 1 )
elif a1 <= a0 and a1 < a2:
imin.append( i - 1 )
a0 = a1
a1 = a2
return ( imin, imax )
#--------------------------------------------------------------------------
def nullStellen( self ):
'gibt die Null-Durchgaenge der Real-Teile'
ans = []
( x0, x1 ) = [ x.real for x in self[ :2 ] ]
xamp = 0.01*max( self.real )
xmax = xmin = 0.
imax = imin = 0
i0 = 1
for i in range( 2, len( self )/2 ):
x2 = self[ i ].real
if abs( x1 ) > xamp:
if x1 > x0 and x1 > x2:
xmax = x1
imax = i0
if xmax*xmin < 0. and xmax - xmin > xamp:
y0 = xmin
for j in range( imin + 1, imax + 1 ):
y1 = self[ j ].real
if y0*y1 < 0.:
ans.append( j + y1/( y0 - y1 ) )
break
y0 = y1
if x1 < x0 and x1 < x2:
xmin = x1
imin = i0
if xmax*xmin < 0. and xmax - xmin > xamp:
y0 = xmax
for j in range( imax + 1, imin + 1 ):
y1 = self[ j ].real
if y0*y1 < 0.:
ans.append( j + y1/( y0 - y1 ) )
break
y0 = y1
x0 = x1
x1 = x2
i0 = i
return Vec( ans )
#--------------------------------------------------------------------------
def phase( self, other ):
'gibt die phase zwischen diesem Signal und other'
p = self
q = other
return ( p.conjugate()*q ).real**2/( p.conjugate()*p
*q.conjugate()*q ).real
#--------------------------------------------------------------------------
def kohaerenz( self, other ):
'gibt die kohaerenz. Ohne vorherige Mittelung ist sie konstant 1.'
pq = self.conjugate()*other
return ( ( pq.conjugate()*pq )/( self.conjugate()*self
*other.conjugate()*other ) ).real
#--------------------------------------------------------------------------
def kohaerenzMittel( self, other ):
'gibt die kohaerenz in 400 mal laengerer Zeit als korrelation'
l = len( self )
l2 = l/2
l4 = l/4
w1 = self.wavelet( 0, l4 )
w2 = self.wavelet( l4, l4 )
w3 = self.wavelet( l2, l4 )
w4 = self.wavelet( l2+l4, l4 )
w5 = self.wavelet( l-1, l4 )
p = zip( *[ ( self*w ).fft() for w in ( w1, w2, w3, w4, w5 ) ] )
q = zip( *[ ( other*w ).fft() for w in ( w1, w2, w3, w4, w5 ) ] )
Gpq = l*[ 0. ]
for i in range( l ):
pi = Vec( p[ i ] )
qi = Vec( q[ i ] )
pp = sum( pi.conjugate()*pi ).real
pq = sum( pi.conjugate()*qi )
qq = sum( qi.conjugate()*qi ).real
Gpq[ i ] = ( pq.conjugate()*pq ).real/( pp*qq )
return Signal( Gpq )
#--------------------------------------------------------------------------
def fft( self ):
return Signal( ( 2./len( self ) )*scipy.fftpack.fft( self ) )
#--------------------------------------------------------------------------
def frequenzen( self, dt ):
'''Wenn dt die Laenge des Messschriebes dieses Signals ist, dann
wird die zur Fourier-Transformierten gehoerenden Frequenz-Achse
gegeben.'''
return [ i/dt for i in range( len( self ) ) ]
#--------------------------------------------------------------------------
def autoLeistungsSpektrum( self ):
return self.kreuzLeistungsSpektrum( self ).real
#--------------------------------------------------------------------------
def kreuzLeistungsSpektrum( self, other ):
return self.conjugate()*other
#--------------------------------------------------------------------------
def wavelet( self, i0, l0 ):
'gibt eine Glocken-Kurve mit Maximum 1 bei i0 und Breite l0'
return Signal( [ math.exp( -0.04*( i - i0 )*( i - i0 ) )
for i in range( len( self ) ) ] )
#--------------------------------------------------------------------------
def verrauschen( self, ramp ):
return self + Signal( [ random.uniform( -ramp, ramp )
for i in range( len( self ) ) ] )
#--------------------------------------------------------------------------
def ausgabe( self, n=10 ):
'ausgabe von n Werten pro Zeile'
for i in range( 0, len( self ), n ):
print( '%d-%d:'%( i, i+n-1 ), self[ i:i+n ] )
#--------------------------------------------------------------------------
#==============================================================================
class SchwingSignal( Signal ):
'Analyse-Werkzeuge fuer ein Signal aus der Modal-Analyse'
#--------------------------------------------------------------------------
def __new__( cls, data, freq ):
return numpy.asarray( data ).view( cls )
#--------------------------------------------------------------------------
def __init__( self, data, freq ):
self.EMS=Vec( [(2.3196-1.0309j), (2.7737-2.6277j), -5.j,
(-2.2703-2.3784j), (-1.7516-0.9554j) ],
dtype=numpy.complex ).conjugate()
self.freq = freq
#--------------------------------------------------------------------------
def modaleParameter( self, name, ir=0, fout=0 ):
'''ermittelt best Fit Ein-Massen-Schwinger an Stelle ir mit
Gradienten-Verfahren.'''
if fout==2:
self.gnuplotABC( name, 0., 0.01, 1., 'Nur gucken' )
( c, x0, inx0 ) = self.realZero( ir ) # c = w0
f = abs( self[ ir + 1 ] )/abs( self[ ir ] )
h = self.freq[ ir + 1 ]/self.freq[ ir ] - 1.
if fout > 0:
print( 'Frequenz = %s, f=%9.2e, h=%9.2e, fout=%d'%(
self.freq[ ir ], f, h, fout ) )
D = f*h/math.sqrt( 1. - f*f )
b = D*c/math.sqrt( 1. - D*D )
if fout > 0:
print( 'f=%9.2e, h=%9.2e, D=%9.2e, b=%9.2e, c=%9.2e, x0=%s'%(
f, h, D, b, c, x0 ) )
# -j a = xk*( j wk - ( - b + j c ) ) = xk*( j wk + b - j c ) | *j
# a = xk*( c - wk + j b )
# xk.imag = -2 D c**2*a/( 4 D**2 c**4 ) = -a/( 2 D c**2 )
if 0:
a = -2.*D*c**2*x0.imag
else:
a = abs( self[ ir ]*( c - self.freq[ ir ] + 1.0j*b ) )
if x0.imag > 0.:
a = -a
da = db = dc = 0.
if 10:
( a, b, c ) = scipy.optimize.fmin( self.pot, [ a, b, c ] )
if fout > 0:
print( a, b, c, sum(
[ self.dk( a, b, c, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len(self ) ) ] ) )
if fout > 1:
self.gnuplotABC( name, a, b, c, 'Start' )
for iter in range( 10 ):
p0 = Vec( [ a, b, c ] )
( a, b ) = self.anpassungAB( a, b, c, fout )
if fout==2:
self.gnuplotABC( name, a, b, c, p0.norm() )
c = self.anpassungC( a, b, c, fout )
pn = ( Vec( [ a, b, c ] ) - p0 ).norm()
( w0, D ) = self.bc2wd( b, c )
if fout > 0:
print( 'global %d: a=%9.2e, D=%9.2e, w0=%9.2e, |dp|=%9.2e'%(
iter, a, D, w0, pn ) )
if pn < 1.e-3:
break
if fout:
self.gnuplotABC( name, a, b, c, 'Ende' )
return ( b, c )
#--------------------------------------------------------------------------
def bc2wd( self, b, c ):
w = math.sqrt( b**2 + c**2 )
D = b/w
return ( w, D )
#--------------------------------------------------------------------------
def dk( self, a, b, c, freq, xk, fabl=0 ):
'ein Summand der zu minimierenden Fehler-Summe'
j = complex( 0., 1. )
nk = complex( b, freq - c )
nk2 = nk*nk
grad = Vec( [ -j/nk, j*a/nk2, a/nk2 ] )
wert = ( -j*a/nk - xk ).conjugate()
if fabl == 0: # Summe
ans = 0.5*wert*wert.conjugate()
elif fabl == 1: # Vektor
ans = wert*grad
else: # Matrix
nk3 = nk2*nk
grada = Vec( [ 0., j/nk2, 1./nk2 ] )
gradb = Vec( [ j/nk2, -j*2.*a/nk3, -2.*a/nk3 ] )
gradc = Vec( [ 1./nk2, -2.*a/nk3, j*2.*a/nk3 ] )
ans = ( grad.outer( grad )
+ wert*Tensor( [ grada, gradb, gradc ] ) )
return ans.real
#--------------------------------------------------------------------------
def pot( self, param ):
( a, b, c ) = param
return sum( [ self.dk( a, b, c, self.freq[ i ], self[ i ] )
for i in range( len(self ) ) ] )
#--------------------------------------------------------------------------
def anpassungAB( self, a, b, c, fout=0 ):
if fout > 0:
print( 'Anpassung von a und b' )
for iter in range( 10 ):
rhs = sum( [ self.dk( a, b, c, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len( self ) ) ] )
mat = sum( [ self.dk( a, b, c, self.freq[ i ], self[ i ], fabl=2 )
for i in range( len( self ) ) ] )
dp = mat[ :2, :2 ].solve( rhs[ :2 ] )
rn = rhs.norm()
pn = dp.norm()
( w, D ) = self.bc2wd( b, c )
if fout > 0:
print( 'ab %d: a=%9.2e, D=%9.2e, |dp|=%9.2e, |rhs|=%9.2e'%(
iter, a, D, pn, rn ) )
a -= dp[ 0 ]
b -= dp[ 1 ]
if pn < 1.e-3:
break
return ( a, b )
#--------------------------------------------------------------------------
def anpassungC( self, a, b, c, fout=0 ):
if fout > 0:
print( 'Anpassung von c' )
for iter in range( 10 ):
rhs = sum( [ self.dk( a, b, c, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len( self ) ) ] )
mat = sum( [ self.dk( a, b, c, self.freq[ i ], self[ i ], fabl=2 )
for i in range( len( self ) ) ] )
dp = rhs[ 2 ]/mat[ 2, 2 ]
rn = abs( rhs[ 2 ] )
pn = abs( dp )
if fout > 0:
print( 'c %d: c=%9.2e, |dp|=%9.2e, |rhs|=%9.2e'%(
iter, c, pn, rn ) )
c -= dp
if pn < 1.e-3:
break
return c
#--------------------------------------------------------------------------
def gnuplotABC( self, name, a, b, c, pn ):
os.system( '''gnuplot -persist <<EOC
a = %s
b = %s
c = %s
set title "a=%s, b=%s, c=%s, |p|=%s"
set xzeroaxis
re( x ) = a*( c - x )/( b**2 + ( c - x )**2 )
im( x ) = -a*b/( b**2 + ( c - x )**2 )
plot [ %s:%s ] re( x ) w lines t "ReFit", im( x ) w lines t "ImFit", "%s" u 1:2 w lines t "Re", "%s" u 1:3 w lines t "Im"
q
EOC
'''%( a, b, c, a, b, c, pn, self.freq[ 0 ], self.freq[ -1 ], name, name ) )
#--------------------------------------------------------------------------
def realZero( self, ir ):
'''gibt exakten Null-Durchgang des Real-Teils
self = [ 2.39795795 -0.20787867j, 3.03729154 -0.33612649j,
4.12705650 -0.62526639j, 6.34354794 -1.51056344j,
12.17101567 -6.78776367j, -13.48430162-21.j,
-8.82300052 -2.74224471j, -5.39553478 -0.88032262j,
-3.93430780 -0.39995498j, -3.17398594 -0.20012574j ]
'''
wechsel = [ self[ i ]*self[ i + 1 ] < 0
for i in range( len( self ) - 1 )]
if wechsel.count( True ) == 0:
raise ValueError( 'kein Nulldurchgang!!!' )
else:
i0 = -1
di = len( self )
imin = ir
while 1:
try:
i1 = wechsel.index( True, i0+1 )
if abs( i1 - ir ) < di:
imin = i1
di = abs( i1 - ir )
i0 = i1
except:
break
# self[ imin ] + t*( self[ imin+1 ] - self[ imin ] )=0
# t*( self[ imin+1 ] - self[ imin ] )=-self[ imin ]
t = -self[ imin ].real/( self[ imin+1 ] - self[ imin ] ).real
w0 = self.freq[ imin ] + t*( self.freq[ imin+1 ] - self.freq[ imin ] )
x0 = self[ imin ] + t*( self[ imin+1 ] - self[ imin ] )
return ( w0, x0, imin + t )
#--------------------------------------------------------------------------
def numRhs( self, a, b, c ):
eps = 1.e-3
h = max( eps, eps*a )
da = sum( [ 0.5*( self.dk( a + h, b, c, self.freq[ i ], self[ i ] )
- self.dk( a - h, b, c, self.freq[ i ], self[ i ] ) )/h
for i in range( len( self ) ) ] )
h = eps*b
db = sum( [ 0.5*( self.dk( a, b + h, c, self.freq[ i ], self[ i ] )
- self.dk( a, b - h, c, self.freq[ i ], self[ i ] ) )/h
for i in range( len( self ) ) ] )
h = eps*c
dc = sum( [ 0.5*( self.dk( a, b, c + h, self.freq[ i ], self[ i ] )
- self.dk( a, b, c - h, self.freq[ i ], self[ i ] ) )/h
for i in range( len( self ) ) ] )
return Vec( [ da, db, dc ] )
#--------------------------------------------------------------------------
def numTang( self, a, b, c ):
"""f'+ = ( f(x+2h) - f(x) )/2h
f'- = ( f(x) - f(x-2h) )/2h
f'' = ( f'+ - f'- )/2h
= ( f(x+2h) - f(x) )/2h - ( f(x) - f(x-2h) )/2h )/2h
f''= ( f(x+h) - 2 f(x) + f(x-h) )/h**2
"""
eps = 1.e-5
h = max( eps, eps*a )
dp = self.numRhs( a+h, b, c )
dm = self.numRhs( a-h, b, c )
ta = 0.5*( dp - dm )/h
h = eps*b
dp = self.numRhs( a, b+h, c )
dm = self.numRhs( a, b-h, c )
tb = 0.5*( dp - dm )/h
h = eps*c
dp = self.numRhs( a, b, c+h )
dm = self.numRhs( a, b, c-h )
tc = 0.5*( dp - dm )/h
return Tensor( [ ta, tb, tc ] )
#--------------------------------------------------------------------------
def numTang1( self, a, w, D ):
"""f'+ = ( f(x+2h) - f(x) )/2h
f'- = ( f(x) - f(x-2h) )/2h
f'' = ( f'+ - f'- )/2h
= ( f(x+2h) - f(x) )/2h - ( f(x) - f(x-2h) )/2h )/2h
f''= ( f(x+h) - 2 f(x) + f(x-h) )/h**2
"""
eps = 1.e-5
h = eps*a
d0 = sum( [ self.dk( a, w, D, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len( self ) ) ] )
dp = sum( [ self.dk( a+h, w, D, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len( self ) ) ] )
dm = sum( [ self.dk( a-h, w, D, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len( self ) ) ] )
ta = 0.5*( dp - dm )/h
#ta[ 0 ] = ( dp[ 0 ] + dm[ 0 ] - 2.*d0[ 0 ] )/h
h = eps*w
dp = sum( [ self.dk( a, w+h, D, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len( self ) ) ] )
dm = sum( [ self.dk( a, w-h, D, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len( self ) ) ] )
tw = 0.5*( dp - dm )/h
#tw[ 1 ] = ( dp[ 1 ] + dm[ 1 ] - 2.*d0[ 1 ] )/h
h = eps*D
dp = sum( [ self.dk( a, w, D+h, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len( self ) ) ] )
dm = sum( [ self.dk( a, w, D-h, self.freq[ i ], self[ i ], fabl=1 )
for i in range( len( self ) ) ] )
tD = 0.5*( dp - dm )/h
#tD[ 2 ] = ( dp[ 2 ] + dm[ 2 ] - 2.*d0[ 2 ] )/h
return Tensor( [ ta, tw, tD ] )
#--------------------------------------------------------------------------
#==============================================================================
class Vec_V2R( numpy.ndarray ):
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0. ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data, dtype=dtype, subok=1 )
).view( cls )
#--------------------------------------------------------------------------
def norm( self ):
return math.sqrt( self.dot( self ).real )
#--------------------------------------------------------------------------
#==============================================================================
class Vec_V2V( numpy.ndarray ):
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0. ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data, dtype=dtype, subok=1 )
).view( cls )
#--------------------------------------------------------------------------
def normalize( self ):
'gibt den richtungsvektor'
a = self.norm()
return Vec( self/a, dtype=self.dtype ).view( self.__class__ )
#--------------------------------------------------------------------------
#==============================================================================
class Vec_VV2R( numpy.ndarray ):
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0. ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data, dtype=dtype, subok=1 )
).view( cls )
#--------------------------------------------------------------------------
def dot( self, other ):
'komplexes Skalar-Produkt'
return sum( self.conjugate()*other )
#--------------------------------------------------------------------------
def dist( u, v ):
return ( u - v ).norm()
#--------------------------------------------------------------------------
def mac( u, v ):
return u.dot( v )/math.sqrt( ( u.dot( u )*v.dot( v ) ).real )
#--------------------------------------------------------------------------
def angle( u, v ):
'Winkel zwischen diesem Vektor und einem anderen other'
return math.acos( min( 1., max( -1., u.mac( v ) ) ) )
#--------------------------------------------------------------------------
#==============================================================================
class Vec_VV2V( numpy.ndarray ):
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0. ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data, dtype=dtype, subok=1 )
).view( cls )
#--------------------------------------------------------------------------
def cross( self, *veks ):
'''Kreuz-Produkt in n=2,3,4,... Dimensionen benoetigt n-2 Argumente
v x A = v x ( ai o ei ) = ( v x ai ) o ei
= ( ei o ( v x ai ) )t = -( ei o ( ai x v ) )t
= -( ( ei o ai ) x v )t
'''
if self.size == 2:
return Vec( [ -self[ 1 ], self[ 0 ] ] )
elif self.size == 3:
other = veks[ 0 ]
if other.size == 3:
ans = [ self[ 1 ]*other[ 2 ] - self[ 2 ]*other[ 1 ],
self[ 2 ]*other[ 0 ] - self[ 0 ]*other[ 2 ],
self[ 0 ]*other[ 1 ] - self[ 1 ]*other[ 0 ] ]
else: # v x A
( a1, a2, a3 ) = other.tr() # Spaltenvektoren
return Tensor( [ self.cross( a1 ), self.cross( a2 ),
self.cross( a3 ) ], dtype=self.dtype ).tr()
else:
mat = numpy.array( [ self ] + list( veks ) )
ans = []
sgn = -( -1 )**self.size
for i in range( self.size ):
inx = range( self.size )
inx.remove( i )
ans.append( sgn*numpy.linalg.det( mat[ :, inx ] ) )
sgn = -sgn
return Vec( ans, dtype=self.dtype )
#--------------------------------------------------------------------------
#==============================================================================
class Vec_V2T( numpy.ndarray ):
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0. ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data, dtype=dtype, subok=1 )
).view( cls )
#--------------------------------------------------------------------------
def diag( self, k=0 ):
'''Gibt die Diagonal-Matrix mit dieser k-ten-Diagonalen.
k=0: Haupt-Diagonale'''
return Tensor( numpy.diag( self, k ), dtype=self.dtype )
#--------------------------------------------------------------------------
def voigt( self ):
'gibt den symmetrischen Tensor mit dieser voigt-Notation'
T = Tensor()
for k in range( 6 ):
( i,j ) = self.VOIGT[ k ]
T[ i,j ] = T[ j,i ] = self[ k ]
return T
#--------------------------------------------------------------------------
def axialerTensor( u, v=[] ):
'gibt axialen Tensor A von u so, dass A(u)*w = u x w fuer alle w gilt'
if u.size == 3:
( x, y, z ) = u
ans = Tensor( [ [ 0., -z, y ],
[ z, 0., -x ],
[ -y, x, 0. ] ], dtype=u.dtype )
elif u.size == 4:
( ux, uy, uz, ut ) = u
( vx, vy, vz, vt ) = v
ans = Tensor( [[ 0., uz*vt - ut*vz, ut*vy - uy*vt, uy*vz - uz*vy ],
[ ut*vz - uz*vt, 0., ux*vt - ut*vx, uz*vx - ux*vz ],
[ uy*vt - ut*vy, ut*vx - ux*vt, 0., ux*vy - uy*vx ],
[ uz*vy - uy*vz, ux*vz - uz*vx, uy*vx - ux*vy, 0. ]
], dtype=u.dtype )
else:
raise ValueError('axialerTensor nur in 3-4 Dimensionen definiert.')
return ans
#--------------------------------------------------------------------------
def drehTensor( self ):
'''gibt den orthogonalen Tensor, der um diese Achse mit diesem Betrag
dreht. Ist nur in 3 Dimensionen definiert.'''
t = Eye( dtype=self.dtype )
alp = self.norm()
if alp > 0.:
ud = self.normalize().axialerTensor()
t = ( t + math.sin( alp )*ud
+ ( 1. - math.cos( alp ) )*ud.mmul( ud ) )
return t
#--------------------------------------------------------------------------
def drehSpiegelTensor( self ):
'''gibt den orthogonalen Tensor, der um diese Achse mit diesem Betrag
dreht. Ist nur in 3 Dimensionen definiert.'''
t = Eye( -1., dtype=self.dtype )
alp = self.norm()
if alp > 0.:
ud = self.normalize().axialerTensor()
t = ( t + math.sin( alp )*ud
- ( 1. - math.cos( alp ) )*ud.mmul( ud ) )
return t
#--------------------------------------------------------------------------
#==============================================================================
class Vec_VV2T( numpy.ndarray ):
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0. ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data, dtype=dtype, subok=1 )
).view( cls )
#--------------------------------------------------------------------------
def mmul( self, other ):
'gibt das Matrixprodukt dieses Vekors mit other ( Tensor oder Vec )'
ans = numpy.dot( self, other )
if other.size == self.size**2: # Vec.mmul( Tensor )
ans = Vec( ans, dtype=self.dtype )
return ans
#--------------------------------------------------------------------------
def outer( u, v ):
'komplexes dyadisches Produkt'
return Tensor( numpy.outer( u, v.conjugate() ), dtype=u.dtype )
#--------------------------------------------------------------------------
def wedge( u, v ):
'gibt den Tensor u o v - v o u'
t = u.outer( v )
return t - t.H()
#--------------------------------------------------------------------------
def proj( u, v ):
'''gibt den Tensor T, fuer den gilt
Tu = u, Tv = v, Tw = 0 fuer alle w mit w*u=w*v=0'''
t = u.wedge( v )
t = 0.5*( t.dot( t ) )
uv = u.outer( ( u.dot( v )/t )*v )
return ( u.outer( ( v.dot( v )/t )*u ) - uv - uv.tr()
+ v.outer( ( u.dot( u )/t )*v ) )
#--------------------------------------------------------------------------
def exp( u, v, alp ):
'''gibt den orthogonalen Tensor, der in der u-v-Ebene mit alp
dreht. Ist in beliebig vielen Dimensionen definiert.'''
return ( Eye( n=len( u ), dtype=u.dtype )
- ( ( 1. - math.cos( alp ) )*u ).proj( v )
+ math.sin( alp )*u.wedge( v ) )
#--------------------------------------------------------------------------
def basis( self, other ):
'komplettiert self=v und other=[v] zu einem Basis-System'
b = [ self.normalize() ] + [ u.normalize() for u in other ]
n = len( self )
e0 = n*[ 0. ]
no = []
s = sum( [ abs( v ) for v in b ] )
for j in range( n - len( b ) ):
t = [ ( s[ i ], i ) for i in range( n ) ]
t.sort()
for ( x, i ) in t:
if not i in no:
break
no.append( i )
e = Vec( e0[ : ] )
e[ i ] = 1.
e = ( e - sum( [ ei.dot( e )*ei for ei in b ] ) ).normalize()
b.append( e )
s = s + abs( e )
return [ self ] + other + b[ len( other ) + 1: ]
#--------------------------------------------------------------------------
#==============================================================================
class Vec( Vec_V2R, Vec_V2V, Vec_V2T, Vec_VV2R, Vec_VV2V, Vec_VV2T ):
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0. ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data, dtype=dtype, subok=1 )
).view( cls )
#--------------------------------------------------------------------------
def __array_finalize__(self, obj):
pass
#--------------------------------------------------------------------------
#==============================================================================
class Vec2( Vec ):
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0. ], dtype=numpy.float ):
return Vec( data=data, dtype=dtype )
#--------------------------------------------------------------------------
def orth( self, dum=0. ):
return Vec2( [ -self[ 1 ], self[ 0 ] ] )
#--------------------------------------------------------------------------
#==============================================================================
class CVec( Vec ):
'komplexer Vektor'
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0. ], dtype=numpy.complex ):
return Vec( data=data, dtype=numpy.complex )
#--------------------------------------------------------------------------
#==============================================================================
class ZufallsVektor( Vec ):
'Ein Vec mit zufaelligen Koordinaten'
#--------------------------------------------------------------------------
def __new__( self, von=-1., bis=1., dtype=numpy.float ):
return Vec( data=[ random.uniform( von, bis ) for i in (1,2,3) ],
dtype=dtype )
#--------------------------------------------------------------------------
#==============================================================================
class Tensor_T2R( numpy.ndarray ):
'komplexer Tensor'
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0. ],
[ 0., 0., 0. ],
[ 0., 0., 0. ] ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def spur( self ):
'''gibt die Spur dieses Tensors Txx + Tyy + Tzz'''
return numpy.trace( self )
#--------------------------------------------------------------------------
def invar2( self ):
"""invar2() -> float: Gibt die zweite Invariante
Txx*Tyy + Tyy*Tzz + Txx*Tzz - Txy**2 - Txz**2 - Tyz**2
"""
sp = self.spur()
return 0.5*( sp*sp - self.mmul( self ).spur() )
#--------------------------------------------------------------------------
def det( self ):
"""det() -> float: Gibt die Determinante dieses Tensors
Txx*( Tyy*Tzz - Tyz**2 ) + Txy*( Tyz*Txz - Txy*Tzz ) + (
Txz*( Txy*Tyz - Txz*Tyy ) )
"""
return numpy.linalg.det( self )
#--------------------------------------------------------------------------
def invar3( self ):
"invar3() -> float: Gibt die Determinante dieses Tensors"
return self.det()
#--------------------------------------------------------------------------
def vMises( self ):
"vMises() -> float: Gibt die von Mises Vergleichsspannung"
Td = self.dev()
return math.sqrt( 1.5*Td.dot( Td ) )
#--------------------------------------------------------------------------
def norm( self ):
'gibt die 2-Norm sum( |Tij|**2 ) dieses Tensors'
return math.sqrt( self.dot( self ).real )
#--------------------------------------------------------------------------
#==============================================================================
class Tensor_T2V( numpy.ndarray ):
'komplexer Tensor'
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0. ],
[ 0., 0., 0. ],
[ 0., 0., 0. ] ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def diag( self ):
'gibt den Vektor mit den Diagonal-Elementen'
return Vec( self.diagonal() )
#--------------------------------------------------------------------------
def voigt( self ):
'gibt fuer einen symmetrischen Tensor dessen Voigt-Notation'
return Vec( [ self[ i,j ] for ( i,j ) in self.VOIGT ] )
#--------------------------------------------------------------------------
def voigt2( self ):
'gibt fuer einen Dehungstensor dessen Voigt-Notation'
return Vec( [ ( 1. + ( i!=j ) )*self[ i,j ]
for ( i,j ) in self.VOIGT ] )
#--------------------------------------------------------------------------
def drehAchse( self ):
'''gibt die Dreh-Vektor (Winkel*Achse) fuer den Fall, dass dies ein
orthogonaler Tensor ist'''
# nicht sin(alp)*\hat{n}!
ax = Vec( [ self[ 2, 1 ] - self[ 1, 2 ],
self[ 0, 2 ] - self[ 2, 0 ],
self[ 1, 0 ] - self[ 0, 1 ] ] )
if ax.norm() == 0.:
ans = ax
else:
alp = math.acos( min( 1., max( -1., 0.5*self.spur() - 0.5 ) ) )
ans = ( alp/ax.norm() )*ax
return ans
#--------------------------------------------------------------------------
def drehVektor( self ):
'''gibt die Dreh-Vektor (Winkel*Achse) fuer den Fall, dass dies ein
orthogonaler Tensor ist'''
return self.drehAchse()
#--------------------------------------------------------------------------
def axialerVektor( A ):
'''gibt den axialen Vektor des schiefsymmetrischen Anteils
diese Tensors'''
As = A - A.tr()
return Vec( [ 0.5*As[ 2, 1 ], 0.5*As[ 0, 2 ], 0.5*As[ 1, 0 ] ] )
#--------------------------------------------------------------------------
def av( A ): return A.axialerVektor()
#--------------------------------------------------------------------------
def i( A ):
'''gibt die Vektorinvariante dieses Tensors'''
return Vec( [ A[1,2] - A[2,1], A[2,0] - A[0,2], A[0,1] - A[1,0] ] )
#--------------------------------------------------------------------------
#==============================================================================
class Tensor_T2T( numpy.ndarray ):
'komplexer Tensor'
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0. ],
[ 0., 0., 0. ],
[ 0., 0., 0. ] ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def submatrix( self, i0, j0=None ):
'gibt diesen tensor ohne zeile i0 und spalte j0'
if j0 == None:
j0 = i0
ans = []
for i in range( len( self ) ):
if i != i0:
zeile = self[ i ]
inx = range( len( zeile ) )
del inx[ j0 ]
ans.append( zeile[ inx ] )
return self.__class__( ans )
#--------------------------------------------------------------------------
def tr( self ):
'gibt diesen Tensor transponiert'
return self.transpose()
#--------------------------------------------------------------------------
def H( self ):
'gibt diesen Tensor conjugiert transponiert'
return self.transpose().conjugate()
#--------------------------------------------------------------------------
def normalize( self ):
'gibt eine Kopie dieses Tensors mit Betrag 1'
return self/self.norm()
#--------------------------------------------------------------------------
def inv( self ):
'Invertierung'
return numpy.linalg.inv( self )
#--------------------------------------------------------------------------
def cof( A ):
'Kofaktor'
return 0.5*A.touter( A )
#--------------------------------------------------------------------------
def adj( A ):
'Adjunkte'
return A.cof().tr()
#--------------------------------------------------------------------------
def sym( self ):
'gibt den symmetrischen Anteil dieses Tensors'
return 0.5*( self + self.tr() )
#--------------------------------------------------------------------------
def skw( self ):
'gibt den schief-symmetrischen Anteil dieses Tensors'
return 0.5*( self - self.tr() )
#--------------------------------------------------------------------------
def sph( self ):
'gibt den Kugel-Anteil dieses Tensors'
return Eye( self.spur()/len( self ), dtype=self.dtype )
#--------------------------------------------------------------------------
def dev( self ):
'gibt den Deviator dieses Tensors'
return self - self.sph()
#--------------------------------------------------------------------------
def wurzel( self, vz=[1,1,1] ):
( a, b, c ) = vz
( l, v ) = self.eig()
return ( v[ 0 ].outer( a*math.sqrt( l[ 0 ] )*v[ 0 ] )
+ v[ 1 ].outer( b*math.sqrt( l[ 1 ] )*v[ 1 ] )
+ v[ 2 ].outer( c*math.sqrt( l[ 2 ] )*v[ 2 ] ) )
#--------------------------------------------------------------------------
#==============================================================================
class Tensor_TV2T( numpy.ndarray ):
'komplexer Tensor'
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0. ],
[ 0., 0., 0. ],
[ 0., 0., 0. ] ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def cross( A, v ):
'''gibt das Kreuzprodukt zwischen Tensor A und Vektor v
( ei o ai ) x v = ei o ( ai x v )'''
( a1, a2, a3 ) = A # Zeilenvektoren
return Tensor( [ a1.cross( v ), a2.cross( v ), a3.cross( v ) ] )
#--------------------------------------------------------------------------
#==============================================================================
class Tensor_TT2R( numpy.ndarray ):
'komplexer Tensor'
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0. ],
[ 0., 0., 0. ],
[ 0., 0., 0. ] ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def dot( self, other ):
'''gibt das Skalarprodukt des konjugiert komplexen dieses Tensors mit
other'''
return self.H().mmul( other ).spur()
#--------------------------------------------------------------------------
#==============================================================================
class Tensor_TT2V( numpy.ndarray ):
'komplexer Tensor'
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0. ],
[ 0., 0., 0. ],
[ 0., 0., 0. ] ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def dotcross( A, B ):
'''gibt das skalar Kreuzprodukt
A = ei o ai, B = ej o bj -> A.xB=(ai.ej)ei x bj
A = ai o ei, B = ej o bj -> A.xB=ai x bi
i(A)=E3:A=1.xA
A.xB=i(A.B)=E3:(A.B)
'''
At = A.tr()
return sum( [ At[ j ].cross( B[ j ] ) for j in ( 0,1,2 ) ], Vec() )
#--------------------------------------------------------------------------
def tcross( A, B ):
'''gibt das Kreuzprodukt von Tensoren
AxB = E3(A.Bt) = i(A.Bt) = Aik*Bjk*ei x ej = A.xBt
(ejxek) o ej o ek : Alm*Bnm*el o en
= Aik*Bjk*ei x ej = Aik*ei x Bjk*ej
= A.xBt
'''
At = A.tr()
Bt = B.tr()
return sum( [ At[ j ].cross( Bt[ j ] ) for j in ( 0,1,2 ) ], Vec() )
#--------------------------------------------------------------------------
#==============================================================================
class Tensor_TT2T( numpy.ndarray ):
'komplexer Tensor'
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0. ],
[ 0., 0., 0. ],
[ 0., 0., 0. ] ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def mmul( self, other ):
'gibt das Matrixprodukt dieses Tensors mit other ( Tensor oder Vec )'
ans = numpy.dot( self, other )
lot = len( other )
if hasattr( other, 'size' ):
lot = other.size
if lot == self.size:
ans = Tensor( ans, dtype=self.dtype )
elif lot**2 == self.size:
ans = Vec( ans, dtype=self.dtype )
return ans
#--------------------------------------------------------------------------
def touter( A, B ):
'gibt das aeussere Tensorprodukt dieses Tensors mit B'
return ( Eye( A.spur()*B.spur() - A.mmul( B ).spur(), A.shape[ 0 ] )
+ ( B.mmul( A ) + A.mmul( B ) - A.spur()*B - B.spur()*A ).tr() )
#--------------------------------------------------------------------------
def crosscross( A, B ):
'''gibt das doppelte Kreuzprodukt zwischen Tensor A und B
(ei o ai)xx(bj o ej) = (ai x bj) o ( ei x ej )
'''
( a1, a2, a3 ) = A
( b1, b2, b3 ) = B.tr()
( e1, e2, e3 ) = Eye()
return ( ( a2.cross( b3 ) - a3.cross( b2 ) ).outer( e1 )
+ ( a3.cross( b1 ) - a1.cross( b3 ) ).outer( e2 )
+ ( a1.cross( b2 ) - a2.cross( b1 ) ).outer( e3 ) )
#--------------------------------------------------------------------------
#==============================================================================
class Tensor_TT2T4( numpy.ndarray ):
'komplexer Tensor'
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0. ],
[ 0., 0., 0. ],
[ 0., 0., 0. ] ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def outer( A, B ):
'dyadisches Produkt mit anderem Tensor'
ans = Tensor4()
for i4 in range( 9 ):
( i,j ) = ans.IX[ i4 ]
for j4 in range( 9 ):
( k,l ) = ans.IX[ j4 ]
ans[ i4, j4 ] = A[ i,j ]*B[ k,l ]
return ans
#--------------------------------------------------------------------------
#==============================================================================
class Tensor( Tensor_T2R, Tensor_T2V, Tensor_T2T,
Tensor_TV2T,
Tensor_TT2R, Tensor_TT2V, Tensor_TT2T, Tensor_TT2T4 ):
'komplexer Tensor'
VOIGT = voigt
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0. ],
[ 0., 0., 0. ],
[ 0., 0., 0. ] ], dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def __array_finalize__(self, obj):
pass
#--------------------------------------------------------------------------
def koeff( self ):
'gibt die neun Koeffizienten als list'
return list( self[ 0 ] ) + list( self[ 1 ] ) + list( self[ 2 ] )
#--------------------------------------------------------------------------
def __getitem__( self, i ):
ans = numpy.ndarray.__getitem__( self, i )
if len( ans.shape ) == 1:
ans = Vec( ans, dtype=self.dtype )
return ans
#--------------------------------------------------------------------------
def solve( self, rhs ):
'gibt Loesung x von self.mmul( x ) = rhs'
return numpy.linalg.solve( self, rhs )
#--------------------------------------------------------------------------
def funktion( self, f ):
( l, v ) = self.eig()
return ( v[ 0 ].outer( f( l[ 0 ] )*v[ 0 ] )
+ v[ 1 ].outer( f( l[ 1 ] )*v[ 1 ] )
+ v[ 2 ].outer( f( l[ 2 ] )*v[ 2 ] ) )
#--------------------------------------------------------------------------
def polar( F ):
'F = R U = v R'
U = F.tr().mmul( F ).funktion( math.sqrt )
R = F.mmul( U.inverse() )
v = F.mmul( R.tr() )
return ( U, v, R )
#--------------------------------------------------------------------------
def svd( self ):
return numpy.linalg.svd( self )
#--------------------------------------------------------------------------
def eig( self ):
'''gibt die Eigen-Werte und -Vektoren'''
( ew, ev ) = numpy.linalg.eig( self )
return ( ew, ev.tr() )
#--------------------------------------------------------------------------
def eigenValues( self ):
'gibt die Eigen-Werte'
return numpy.linalg.eig( self )[ 0 ]
#--------------------------------------------------------------------------
def rechtsSystem( self ):
'Macht aus den Zeilen-Vektoren ein Rechts-System'
spalten = Tensor( [ x.normalize() for x in self ] ).tr()
zeilen = []
d = self.dtype
tol = 1.e-6*abs( self.det() )
inx = range( len( self ) )
all = inx[ : ]
for i in all:
x = abs( spalten[ i ] )
s = [ ( x[ j ], j ) for j in all ]
s.sort()
s.reverse()
for ( x, j ) in s:
if j in inx:
zeile = self[ j ]
t = Tensor( zeilen + [ zeile ],
dtype=d )[ :i+1, :i+1 ]
d = t.det()
if d > tol:
inx.remove( j )
if d < 0.:
zeile = -zeile
zeilen.append( zeile )
break
else:
print( t )
print( 'versagt', j, s )
else:
raise Exception( 'kein Rechts-System gefunne' )
print( Tensor( zeilen, dtype=d ) )
return Tensor( zeilen, dtype=d )
#--------------------------------------------------------------------------
def hessenberg( self ):
'''Householder reduction to Hessenberg Form (zeroes below the diagonal)
while keeping the same eigenvalues as self.'''
ans = self
q = one = Eye( n=self.shape[ 1 ] )
for i in range( self.shape[ 1 ] - 2 ):
st = ans.tr()
( v, beta ) = st[ i ].house( i + 1 )
q = ( one - v.outer( beta*v ) ).mmul( q )
ans = ans - v.outer( st.mmul( v )*beta ) # ans = ( 1 - b*v o v )*s
ans = ans - ans.mmul( v ).outer( v*beta )# ans = s*( 1 - b*v o v )
return ( ans, q )
#--------------------------------------------------------------------------
#==============================================================================
class ZufallsTensor( Tensor ):
'Ein Tensor mit zufaelligen Koordinaten'
#--------------------------------------------------------------------------
def __new__( self, von=-1., bis=1., dtype=numpy.float ):
return Tensor( data=[ [ random.uniform( von, bis ) for j in (1,2,3 ) ]
for i in (1,2,3 ) ], dtype=dtype )
#--------------------------------------------------------------------------
#==============================================================================
class Eye( Tensor ):
'Vielfaches des Einheitstensors, Diagonalmatrix: siehe Vec.diag, Diag'
#--------------------------------------------------------------------------
def __new__( self, x=1., m=3, dtype=numpy.float ):
return Tensor( data=[ [ x*( i == j ) for j in range( m ) ]
for i in range( m ) ], dtype=dtype )
#--------------------------------------------------------------------------
#==============================================================================
class Diag( Tensor ):
'Diagonal-Matrix'
#--------------------------------------------------------------------------
def __new__( self, d=[], dtype=numpy.float ):
num = range( len( d ) )
return Tensor( data=[ [ ( i==j and d[ i ] or 0. ) for j in num ]
for i in num ], dtype=dtype )
#--------------------------------------------------------------------------
#==============================================================================
class CTensor( Tensor ):
'komplexer Vektor'
#--------------------------------------------------------------------------
def __new__( cls, data=[ 0., 0., 0. ], dtype=numpy.complex ):
return Tensor( data=data, dtype=numpy.complex )
#--------------------------------------------------------------------------
#==============================================================================
class Tensor4( numpy.ndarray ):
'Tensor vierter Stufe'
IX=( (0,0),(0,1),(0,2), (1,0),(1,1),(1,2), (2,0),(2,1),(2,2) )
VOIGT=( 0, 4, 8, 5, 6, 1 )
#--------------------------------------------------------------------------
def __new__( cls, data=[ [ 0., 0., 0., 0., 0., 0., 0., 0., 0., ],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., ],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., ],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., ],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., ],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., ],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., ],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., ],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., ] ],
dtype=numpy.float ):
return numpy.asarray( numpy.array( data ), dtype=dtype ).view(
type=cls, dtype=dtype )
#--------------------------------------------------------------------------
def __array_finalize__(self, obj):
pass
#--------------------------------------------------------------------------
def voigt( self ):
'Voigtsche Notation symmetrischer Tensoren'
ans = Tensor4( [ [ 0. for j in range( 6 ) ] for i in range( 6 ) ] )
i = 0
for iv in self.VOIGT:
j = 0
for jv in self.VOIGT:
ans[ i, j ] = self[ iv, jv ]
ans[ j, i ] = self[ jv, iv ]
j += 1
i += 1
return ans
#--------------------------------------------------------------------------
def tr( T ):
'Transposition'
return T.transpose()
#--------------------------------------------------------------------------
def t12( T ):
'vertauschung der ersten beiden basen ijkl<->jikl'
ans = Tensor4()
for i4 in range( 9 ):
( i,j ) = T.IX[ i4 ]
i4a = T.IX.index( ( j,i ) )
for j4 in range( 9 ):
ans[ i4a, j4 ] = T[ i4,j4 ]
return ans
#--------------------------------------------------------------------------
def t13( T ):
'vertauschung der ersten und dritten basen ijkl<->kjil'
ans = Tensor4()
for i4 in range( 9 ):
( i,j ) = T.IX[ i4 ]
for j4 in range( 9 ):
( k,l ) = T.IX[ j4 ]
i4a = T.IX.index( (k,j) )
j4a = T.IX.index( (i,l) )
ans[ i4a, j4a ] = T[ i4,j4 ]
return ans
#--------------------------------------------------------------------------
def t14( T ):
'vertauschung der ersten und vierten basen ijkl<->ljki'
ans = Tensor4()
for i4 in range( 9 ):
( i,j ) = T.IX[ i4 ]
for j4 in range( 9 ):
( k,l ) = T.IX[ j4 ]
i4a = T.IX.index( (l,j) )
j4a = T.IX.index( (k,i) )
ans[ i4a, j4a ] = T[ i4,j4 ]
return ans
#--------------------------------------------------------------------------
def t23( T ):
'vertauschung der zweiten und dritten basen ijkl<->ikjl'
ans = Tensor4()
for i4 in range( 9 ):
( i,j ) = T.IX[ i4 ]
for j4 in range( 9 ):
( k,l ) = T.IX[ j4 ]
i4a = T.IX.index( (i,k) )
j4a = T.IX.index( (j,l) )
ans[ i4a, j4a ] = T[ i4,j4 ]
return ans
#--------------------------------------------------------------------------
def t24( T ):
'vertauschung der zweiten und vierten basen ijkl<->ilkj'
ans = Tensor4()
for i4 in range( 9 ):
( i,j ) = T.IX[ i4 ]
for j4 in range( 9 ):
( k,l ) = T.IX[ j4 ]
i4a = T.IX.index( (i,l) )
j4a = T.IX.index( (k,j) )
ans[ i4a, j4a ] = T[ i4,j4 ]
return ans
#--------------------------------------------------------------------------
def t34( T ):
'vertauschung der dritten und vierten basen ijkl<->ijlk'
ans = Tensor4()
for i4 in range( 9 ):
for j4 in range( 9 ):
( k,l ) = T.IX[ j4 ]
j4a = T.IX.index( (l,k) )
ans[ i4, j4a ] = T[ i4,j4 ]
return ans
#--------------------------------------------------------------------------
def inv( self ):
'Invertierung'
return numpy.linalg.inv( self )
#--------------------------------------------------------------------------
def mmul( A, B ):
if B.size == 9:
ans = Tensor( numpy.reshape(
numpy.dot( A, numpy.reshape( B, 9 ) ), ( 3,3 ) ) )
else:
ans = Tensor4( numpy.dot( A, B ) )
return ans
#--------------------------------------------------------------------------
#==============================================================================
class Eye4( Eye ):
'Vielfaches des Einheitstensors vierter Stufe'
#--------------------------------------------------------------------------
def __new__( self, x=1., dtype=numpy.float ):
return Tensor4( data=[ [ x*( i == j ) for j in range( 9 ) ]
for i in range( 9 ) ], dtype=dtype )
#--------------------------------------------------------------------------
#==============================================================================
class Deviator4( Tensor4 ):
'Tensor, der den Deviator extrahiert'
#--------------------------------------------------------------------------
def __new__( self, x=1., dtype=numpy.float ):
return Eye4() - Eye( 1./3. ).outer( Eye() )
#--------------------------------------------------------------------------
#==============================================================================
class Raum2D( UserDict ):
'soriert Vektoren ein'
#--------------------------------------------------------------------------
def __init__( self, dl=10., fhg=[ 0, 1 ] ):
'sortiert nodes in Quadrate der Kantenlaenge dl in der fhg-Ebene'
UserDict.__init__( self )
self.dl = dl
self.fhg = fhg
#--------------------------------------------------------------------------
def add( self, v ):
( i, j ) = self.getIj( v )
if self.has_key( i ):
di = self[ i ]
if di.has_key( j ):
di[ j ].append( v )
else:
self[ i ] = { j: [ v ] }
#--------------------------------------------------------------------------
def getNext( self, v ):
'''
Bsp.: di = 2
i0
|
#####
#ooo#
#oxo# - j0
#ooo#
#####
'''
( i0, j0 ) = self.getIj( v )
all = []
di = 0
dmin = 1.e36
vmin = None
ok = 1
while 1:
# pruefung der schale di
for i in range( i0-di, i0+di+1 ):
if self.has_key( i ):
si = self[ i ]
for j in ( j0-di, j0+di ):
if si.has_key( j ):
for u in si[ j ]:
d = ( u - v ).norm()
if d < dmin:
dmin = d
vmin = u
for i in ( i0-di, i0+di ):
if self.has_key( i ):
si = self[ i ]
for j in range( j0-di+1, j0+di ):
if si.has_key( j ):
for u in si[ j ]:
d = ( u - v ).norm()
if d < dmin:
dmin = d
vmin = u
if nmin < 0: # naechste schale
di += 1
elif ok: # noch eine schale
ok = 0
di += 1
else:
break
return vmin
#--------------------------------------------------------------------------
def getIj( self, v ):
ans = []
for i in self.fhg:
ans.append( int( v[ i ]/self.dl ) )
return ans
#--------------------------------------------------------------------------
#==============================================================================
class Raum( UserDict ):
'soriert Vektoren in n-dimensionen ein'
#--------------------------------------------------------------------------
def __init__( self, data={}, dl=1. ):
'sortiert nodes in Quadrate der Kantenlaenge tol in der fhg-Ebene'
UserDict.__init__( self, data )
self.dl = dl
#--------------------------------------------------------------------------
def get( self, v, tol=1.e-4 ):
#......................................................................
def getsd( v, ijk, dj ):
#print 'getsd', ijk
ans = []
j0 = ijk.pop( 0 )
ok = len( ijk ) > 0
for j in range( j0 - dj, j0 + dj + 1 ):
if v.has_key( j ):
if ok:
ans.extend( getsd( v[ j ], ijk[ : ], dj ) )
else:
ans.extend( v[ j ] )
return ans
#......................................................................
for u in getsd( self, self.getIJK( v ), int( tol/self.dl ) + 1 ):
if ( u - v ).norm() < tol:
ans = u
break
else:
ans = None
return ans
#--------------------------------------------------------------------------
def add( self, v ):
#......................................................................
def addsd( v, ijk, xyz ):
j = ijk.pop( 0 )
if len( ijk ) == 0:
if v.has_key( j ):
v[ j ].append( xyz )
else:
v[ j ] = [ xyz ]
elif v.has_key( j ):
addsd( v[ j ], ijk, xyz )
else:
while ijk:
v[ j ] = v = {}
j = ijk.pop( 0 )
v[ j ] = [ xyz ]
#......................................................................
addsd( self, self.getIJK( v ), v )
#--------------------------------------------------------------------------
def getIJK( self, v ):
return [ self.f2i( x ) for x in v ]
#--------------------------------------------------------------------------
def f2i( self, x ):
return int( x/self.dl )
#--------------------------------------------------------------------------
#==============================================================================
class Geometrie:
"Superklasse fuer alle geometrischen Objekte"
TOL=1.e-3
#--------------------------------------------------------------------------
def __init__( self, x0 ):
self.x0 = x0
self.dim = len( x0 )
#--------------------------------------------------------------------------
#==============================================================================
class Gerade( Geometrie ):
"Berechnet Geraden"
#--------------------------------------------------------------------------
def __init__(self, x0, x1 ):
"Definition mittels zweier Punkte"
Geometrie.__init__( self, x0 )
self.g = ( x1 - x0 ).normalize()
#--------------------------------------------------------------------------
def schnittPunkt( a, b ):
'''
x0 + l0*g0 = x1 + l1*g1 + l2*g0xg1 | *n, n = g1 x ( g0 x g1 )
l0*( g0,n ) = ( x1-x0, n )
a.x + x*a.g = b.x + y*b.g | .b.g.cross()
x = ( b.x - a.x ).dot( b.g.cross() )/a.g.dot( b.g.cross() )
'''
if a.dim == 3:
n = b.g.cross( a.g.cross( b.g ) )
l = ( b.x0 - a.x0 ).dot( n )/( a.g.dot( n ) )
elif a.dim == 2:
bgx = b.g.cross()
l = ( b.x0 - a.x0 ).dot( bgx )/a.g.dot( bgx )
return a.x0 + l*a.g
#--------------------------------------------------------------------------
def dist( a, p=[] ):
"Berechnet den Abstand von Vektor p zu dieser Geraden"
if isinstance( p, Vec ):
if len( p ) == 3:
ans = a.g.cross( p - a.x0 ).norm()
else:
ans = a.g.cross().normalize().dot( p - a.x0 )
elif isinstance( p, Gerade ):
# a.x0 + r*a.g + s*a.g x p.g = p.x0 + t*p.g | *( a.g x p.g )=n
# s = ( p.x0 - a.x0 )*n/( n*n )
if len( p.g == 3 ):
n = a.g.cross( p.g ).normalize()
ans = abs( n.dot( a.x0 - p.x0 ) )
else:
n = a.g.cross().normalize()
ans = n.dot( p.g )
if abs( ans ) > 1.e-6:
ans = 0.
else:
ans = ( a.x0 - p.x0 ).dot( n )
return ans
#--------------------------------------------------------------------------
def winkelhalbierende( a, b ):
P = a.schnittPunkt( b )
k1 = Kreis( P + a.g, 1. )
k2 = Kreis( P + b.g, 1. )
( I, Q ) = k1.schnittPunktKreis2D( k2 )
if Q.dist( P ) > I.dist( P ):
I = Q
return Gerade( P, I )
#--------------------------------------------------------------------------
def fusspunkt( g, p ):
'''Projektion von p auf g
p = x0 + l*g + e*n |.g
( p - x0 ).g = l
l = ( p - x0 ).g
'''
return g.x0 + g.g.dot( p - g.x0 )*g.g
#--------------------------------------------------------------------------
#==============================================================================
class Ebene( Geometrie ):
'Berechnungen einer Ebene'
#--------------------------------------------------------------------------
def __init__( self, x0=[], x1=[], x2=[] ):
Geometrie.__init__( self, x0 )
if len( x0 ) == 3:
self.n = ( x1 - x0 ).cross( x2 - x0 ).normalize()
else:
self.n = Vec2( [ 0., 0. ] )
#--------------------------------------------------------------------------
def dist( self, p ):
return ( p - self.x0 ).dot( self.n )
#--------------------------------------------------------------------------
def schnittGerade( a, b ):
'''gibt die Gerade, die sowohl in a als auch b liegt
a.x0 + r (a.n x g) = b.x0 + s (b.n x g) + t g | *a.n
s = ( a.x0 - b.x0 )*a.n/( (b.n x g)*a.n )
'''
g = a.n.cross( b.n ).normalize() # Richtung
s = ( a.x0 - b.x0 ).dot( a.n )/b.n.cross( g ).dot( a.n )
p = b.x0 + s*b.n.cross( g )
return Gerade( p, p + g )
#--------------------------------------------------------------------------
def durchstossPunkt( a, g ):
'''gibt den Durchstosspunkt der Geraden g und dieser Ebene'''
return g.x0 + a.n.dot( a.x0 - g.x0 )/a.n.dot( g.g )*g.g
#--------------------------------------------------------------------------
#==============================================================================
class Dreieck( Ebene ):
'Berechnungen eines Dreiecks'
#--------------------------------------------------------------------------
def __init__( self, x0=[], x1=[], x2=[] ):
Ebene.__init__( self, x0, x1, x2 )
self.x1 = x1
self.x2 = x2
#--------------------------------------------------------------------------
def inhalt( d ):
'gibt den flaecheninhalt'
a = d.x0.dist( d.x1 )
b = d.x1.dist( d.x2 )
c = d.x2.dist( d.x0 )
s = 0.5*( a + b + c )
return math.sqrt( s*( s - a )*( s - b )*( s - c ) )
#--------------------------------------------------------------------------
def umkreis( self ):
'gibt den Umkreis dieses Dreiecks'
( p0, p1, p2 ) = ( self.x0, self.x1, self.x2 )
if len( p0 ) == 2: # in der Ebene
mat = Tensor( [ p1 - p0, p2 - p0 ] )
a = p0.dot( p0 )
rhs = Vec( [ p1.dot( p1 ) - a, p2.dot( p2 ) - a ] )
m = mat.solve( 0.5*rhs )
else:
if ( p1 - p0 ).cross( p2 - p0 ).norm() > 1.e-9:
n = ( p1 - p0 ).cross( p2 - p0 )
mat = Tensor( [ p1 - p0, p2 - p0, n ] )
a = p0.dot( p0 )
rhs = Vec( [ p1.dot( p1 ) - a, p2.dot( p2 ) - a,
2.*n.dot( p0 ) ] )
m = mat.solve( 0.5*rhs )
else:
m = Vec( [ 1.e9, 1.e9, 1.e9 ] )
return Kreis( m, ( p0 - m ).norm() )
#--------------------------------------------------------------------------
def inkreis( self ):
'Gibt den Inkreis dieses Dreiecks'
( A, B, C ) = ( self.x0, self.x1, self.x2 )
gBA = ( B - A ).normalize()
gCB = ( C - B ).normalize()
gCA = ( C - A ).normalize()
if len( A ) == 3:
N = gBA.cross( gCA )
mat = numpy.array( [ gBA - gCA, gCB + gBA, N ] )
rhs = Vec( [ A.dot( gBA - gCA ), B.dot( gCB + gBA ), N.dot( A ) ] )
else:
mat = numpy.array( [ gBA - gCA, gCB + gBA ] )
rhs = Vec( [ A.dot( gBA - gCA ), B.dot( gCB + gBA ) ] )
S = Vec( numpy.linalg.solve( mat, rhs ) )
r = Gerade( A, B ).dist( S )
return Kreis( S, r )
#--------------------------------------------------------------------------
def schwerpunkt( self ):
return ( self.x0 + self.x1 + self.x2 )/3.
#--------------------------------------------------------------------------
#==============================================================================
class Viereck( Dreieck ):
'Berechnungen eines Dreiecks'
#--------------------------------------------------------------------------
def __init__( self, x0=[], x1=[], x2=[], x3=[] ):
Dreieck.__init__( self, x0, x1, x2 )
self.x3 = x3
#--------------------------------------------------------------------------
def umkreis( self ):
'gibt einen Umkreis dieses Vierecks'
( p0, p1, p2 ) = ( self.x0, self.x1, self.x2 )
if len( p0 ) == 2: # in der Ebene
mat = Tensor( [ p1 - p0, p2 - p0 ] )
a = p0.dot( p0 )
rhs = Vec( [ p1.dot( p1 ) - a, p2.dot( p2 ) - a ] )
m = mat.solve( 0.5*rhs )
else:
if ( p1 - p0 ).cross( p2 - p0 ).norm() > 1.e-9:
n = ( p1 - p0 ).cross( p2 - p0 )
mat = Tensor( [ p1 - p0, p2 - p0, n ] )
a = p0.dot( p0 )
rhs = Vec( [ p1.dot( p1 ) - a, p2.dot( p2 ) - a,
2.*n.dot( p0 ) ] )
m = mat.solve( 0.5*rhs )
else:
m = Vec( [ 1.e9, 1.e9, 1.e9 ] )
return Kreis( m, ( p0 - m ).norm() )
#--------------------------------------------------------------------------
def inkreis( self ):
'Gibt einen Inkreis dieses Vierecks'
( A, B, C ) = ( self.x0, self.x1, self.x2 )
gBA = ( B - A ).normalize()
gCB = ( C - B ).normalize()
gCA = ( C - A ).normalize()
if len( A ) == 3:
N = gBA.cross( gCA )
mat = numpy.array( [ gBA - gCA, gCB + gBA, N ] )
rhs = Vec( [ A.dot( gBA - gCA ), B.dot( gCB + gBA ), N.dot( A ) ] )
else:
mat = numpy.array( [ gBA - gCA, gCB + gBA ] )
rhs = Vec( [ A.dot( gBA - gCA ), B.dot( gCB + gBA ) ] )
S = Vec( numpy.linalg.solve( mat, rhs ) )
r = Gerade( A, B ).dist( S )
return Kreis( S, r )
#--------------------------------------------------------------------------
def schwerpunkt( self ):
( A, B, C ) = ( self.x0, self.x1, self.x2 )
return A + ( 2./3. )*( 0.5*( B + C ) - A )
#--------------------------------------------------------------------------
#==============================================================================
class Polygon:
#--------------------------------------------------------------------------
def __init__( self, points ):
self.points = points
self.nump = len( points )
#--------------------------------------------------------------------------
def inhalt( self ):
( x0, y0 ) = self.points[ -1 ]
A = 0.
for ( x1, y1 ) in self.points:
A += x0*y1 - x1*y0
( x0, y0 ) = ( x1, y1 )
return 0.5*A
#--------------------------------------------------------------------------
def schwerpunkt( self, k=0 ):
'findet schwerpunkt'
ps = self.points[ k: ] + self.points[ :k ]
( p0, p1 ) = ps[ :2 ]
s = Vec( [ 0,0 ] )
ages = 0.
for p2 in ps[ 2: ]:
d = Dreieck( p0, p1, p2 )
a = d.inhalt()
s = s + d.schwerpunkt()*a
ages += a
p1 = p2
return s/ages
#--------------------------------------------------------------------------
#==============================================================================
class Tetraeder( Dreieck ):
'Berechnungen eines Tetraeders'
#--------------------------------------------------------------------------
def __init__( self, x0=[], x1=[], x2=[], x3=[] ):
Dreieck.__init__( self, x0, x1, x2 )
self.x3 = x3
#--------------------------------------------------------------------------
def umKugelMittelpunkt( self ):
( A, B, C, D ) = ( self.x0, self.x1, self.x2, self.x3 )
mat = Tensor( [ B - A, C - A, D - A ] )
aa = A.dot( A )
rhs = Vec( [ B.dot( B ) - aa, C.dot( C ) - aa, D.dot( D ) - aa ] )
m = 0.5*( mat.solve( rhs ) )
return m
#--------------------------------------------------------------------------
#==============================================================================
class Kreis( Geometrie ):
'Berechnungen eines Kreises'
#--------------------------------------------------------------------------
def __init__( self, x0=[ 0., 0. ], r=0. ):
'Kreis mit Radius r um x0 senkrecht zu n'
Geometrie.__init__( self, x0 )
self.r = r
#--------------------------------------------------------------------------
def umkreis( self, p0, p1, p2 ):
return Dreieck( p0, p1, p2 ).umkreis()
#--------------------------------------------------------------------------
def schnittPunktGerade( self, gerade ):
'''gibt die beiden Schnittpunkte der Geraden gerade mit diesem Kreis,
wobei die Gerade auf die Ebene dieses Kreises projiziert wird.
Passiert die Gerade diesen Kreis, wird ( None, None ) gegeben.
|xg+z*g-x0|=|dx+z*g|=r
(xg-x0 + z*g)*(xg-x0 + z*g) = (xg-x0)**2 + 2*z*(xg-x0)*g + z**2=r^2
z = -(xg-x0)*g +- sqrt( ((xg-x0)*g)**2 + r**2 - (xg-x0)**2 )
'''
ans = [ None, None ]
txt = '''
xg = gerade.x0 = self.x0
( xg + l*g ).( xg + l*g ) == r**2
xg.xg + 2*l*xg.g + l**2*g.g = r**2
l**2 + 2*xg.g/g.g*l + ( xg.xg - r**2 )/g.g = 0
l = -xg.g/g.g +- sqrt( ( xg.g/g.g )**2 + ( r**2 - xg.xg )/g.g )
'''
xg = gerade.x0 - self.x0
gg = gerade.g.dot( gerade.g )
xgg = xg.dot( gerade.g )
xg2 = xg.dot( xg )
w = ( xgg/gg )**2 + ( self.r**2 - xg2 )/gg
if w >= 0.:
w = math.sqrt( w )
l1 = w - xgg/gg
ans = ( gerade.x0 + ( w - xgg/gg )*gerade.g,
gerade.x0 - ( w + xgg/gg )*gerade.g )
return ans
#--------------------------------------------------------------------------
def schnittPunktKreis( k, b ):
'schnittpunkte von Kreis k und b'
v = b.x0 - k.x0
d = v.norm() # Mittelpunktsabstand
ans = [ None, None ]
if d <= k.r + b.r and d >= abs( k.r - b.r ):
f = 0.5*( k.r**2 - b.r**2 + d**2 )/d
p = k.x0 + f/d*v
w = k.r**2 - f**2
if w >= 0.:
h = math.sqrt( w )/d*v.cross()
ans = [ Vec2( p + h ), Vec2( p - h ) ]
return ans
#--------------------------------------------------------------------------
def optKreis( self, P ):
'''minimiert sum( [ ( ( Pj - M ).dot( Pj - M ) - r**2 )**2 ] )'''
N = len( P )
o = Vec( len( P[ 0 ] )*[ 0. ] )
O = o.outer( o )
a = sum( P, o )
z = 2./N*a
b = sum( [ Pj.dot( Pj ) for Pj in P ] )/N
c = sum( [ Pj.dot( Pj )*Pj for Pj in P ], o )
A = sum( [ Pj.outer( Pj + Pj ) for Pj in P ], O ) - a.outer( z )
M = A.inv().mmul( c - b*a )
r = math.sqrt( b - M.dot( z ) + M.dot( M ) )
#r = math.sqrt( sum( [ ( Pj - M ).dot( Pj - M ) for Pj in P ] )/N )
return Kreis( M, r )
#--------------------------------------------------------------------------
#==============================================================================
class Ellipse( Geometrie ):
'Berechnungen einer Ellipse'
#--------------------------------------------------------------------------
def __init__( self, x0=[ 0., 0. ], a=[ 1., 0. ], b=[ 0., 1. ] ):
'Ellipse in x0 mit Hauptachen a und b'
Geometrie.__init__( self, Vec( x0 ) )
self.a = a
self.b = b
#--------------------------------------------------------------------------
def vierPunkte2D( self, p, q, r, s ):
'''Ellipse in der Ebene durch p, q, r und s mit t in der Naehe
u: x0, y0, e0, e1, z
x0 = ( x0, y0 ) ist Zentrum
e1 = ( e0, e1 ) ist erste Halbachse
a = (e0-x0)**2 + (e1-y0)**2
( x, y )x = ( y, -x )
e2 = z*( e0-x0, e1-y0 )x = z*( e1-y0, x0-e0 )
'''
self.p = ( p, q, r, s )
x0 = sum( self.p )/4.
pmin = pmax = None
dmin = 1.e36
dmax = 0.
for p in self.p:
d = p.dist( x0 )
if d < dmin:
dmin = d
pmin = p
if d > dmax:
dmax = d
pmax = p
g1 = pmax - x0
g2 = pmin - x0
( x0, y0, e0, e1, b ) = scipy.optimize.fmin(
self.abstand, [ x0[ 0 ], x0[ 1 ], g1[ 0 ], g1[ 1 ], g2.norm() ] )
self.x0 = Vec( [ x0, y0 ] )
self.a = Vec( [ e0 - x0, e1 - y0 ] )
g1 = Vec( [ e1 - y0, x0 - e0 ] )
self.b = g1*( b/g1.norm() )
#--------------------------------------------------------------------------
def abstand( self, param ):
( x0, y0, e0, e1, b ) = param
p0 = Vec( [ x0, y0 ] )
v1 = Vec( [ e0 - x0, e1 - y0 ] )
a = math.sqrt( v1.dot( v1 ) )
g1 = v1/a
g2 = Vec( [ e1 - y0, x0 - e0 ] ).normalize()
ans = 0.
for p in self.p:
x = ( p - p0 ).dot( g1 )
y = ( p - p0 ).dot( g2 )
ans += ( ( x*b )**2 + ( y*a )**2 - ( a*b )**2 )**2
return ans
#--------------------------------------------------------------------------
#==============================================================================