#!/usr/bin/python3
'geometrie in der komplexen zahlenebene'
from random import uniform
import cmath
cum = [ 0,0 ]
i = complex( 0,1 )
fp = 0
floatfmt = '%1.4f' # format fuer reelle zahlen
#==============================================================================
def realpart( z ):
'fuer wxmaxima'
return z.real
#==============================================================================
def imagpart( z ):
'fuer wxmaxima'
return z.imag
#==============================================================================
def u( x=1,y='' ):
if not isinstance( y, str ):
ans = uniform( x,y )
else:
ans = uniform( -x,x )
return ans
#==============================================================================
class Zahl:
#--------------------------------------------------------------------------
def isNumber( self, x ):
return ( isinstance( x, complex ) or isinstance( x, float )
or isinstance( x, int ) )
#--------------------------------------------------------------------------
#==============================================================================
class Punkt( Zahl ):
#--------------------------------------------------------------------------
def __init__( self, x, y=0 ):
if isinstance( x, Punkt ):
self.z0 = x.z0
else:
self.z0 = complex( x,y )
self.real = self.z0.real
self.imag = self.z0.imag
#--------------------------------------------------------------------------
def spiegel( p, obj ):
'punktspiegelung'
if p.isNumber( obj ):
obj = Punkt( obj )
p2 = p + p
if isinstance( obj, Punkt ):
ans = p2 - obj # p + ( p - obj )
elif isinstance( obj, Kreis ):
ans = Kreis( p2 - obj.z0, obj.r )
elif isinstance( obj, Gerade ):
ans = Gerade( p2 - obj.z0, p2 - obj.z1 )
elif isinstance( obj, Dreieck ):
ans = Dreieck( p2 - obj.z0, p2 - obj.z1, p2 - obj.z2 )
elif hasattr( obj, '__len__' ): # liste
ans = [ p.spiegel( q ) for q in obj ]
else:
raise NotImplementedError( 'Punktspiegelung fuer %s'%type( obj )
+ 'nicht implementiert' )
return ans
#--------------------------------------------------------------------------
def normalize( self ):
r = abs( self )
return Punkt( self.real/r, self.imag/r )
#--------------------------------------------------------------------------
def outer( u,v ):
'u=(a,b), v=(c,d) -> ( ( a*c, a*d ), ( b*c, b*d ) )'
return Matrix( u.z0.real*v, u.z0.imag*v )
#--------------------------------------------------------------------------
def angle( self, b ):
'gibt den winkel zwischen diesen "vektoren"'
if self.isNumber( b ):
b = complex( b )
elif isinstance( b, Punkt ):
b = b.z0
return cmath.phase( b/self.z0 )
#--------------------------------------------------------------------------
def dot( a, b ):
'a.real*b.real + a.imag*b.imag'
return ( a.conjugate()*b ).real
#--------------------------------------------------------------------------
def cross( a, b ):
'( a x b ).ez = ax*by - ay*bx'
return ( a.conjugate()*b ).imag
#--------------------------------------------------------------------------
def dist( a, b ):
'abstand'
return abs( a - b )
#--------------------------------------------------------------------------
def p2t( self ):
'gibt tuple'
return ( self.real, self.imag )
#--------------------------------------------------------------------------
def toVec( self ):
'gibt mf.Vec'
return mf.Vec( [ self.real, self.imag ] )
#--------------------------------------------------------------------------
def toNpVec( self ):
'gibt nf.Vec'
return nf.Vec( [ self.real, self.imag ] )
#--------------------------------------------------------------------------
def conjugate( self ):
return Punkt( self.z0.conjugate() )
#--------------------------------------------------------------------------
def __str__( self ):
'self.z0'
fmt = 'Punkt( %s,%s )'%( floatfmt, floatfmt )
return fmt%( self.real, self.imag )
#--------------------------------------------------------------------------
def __neg__( self ):
'-self'
return Punkt( -self.z0 )
#--------------------------------------------------------------------------
def __abs__( self ):
'abs( self )'
return abs( self.z0 )
#--------------------------------------------------------------------------
def __eq__( a,b ):
'a==b'
ans = 0
if isinstance( b, Punkt ):
ans = a.z0 == b.z0
else:
ans = a.z0 == b
return ans
#--------------------------------------------------------------------------
def __add__( self, z ):
'self + z'
if isinstance( z, Punkt ):
ans = Punkt( self.z0 + z.z0 )
else:
ans = Punkt( self.z0 + z )
return ans
#--------------------------------------------------------------------------
def __sub__( self, z ):
'self - z'
return self + ( -z )
#--------------------------------------------------------------------------
def __mul__( self, z ):
'self*z'
if isinstance( z, Punkt ):
ans = Punkt( self.z0*z.z0 )
else:
ans = Punkt( self.z0*z )
return ans
#--------------------------------------------------------------------------
def __div__( self, z ):
print( 'self//z' )
return self.__truediv__( z )
#--------------------------------------------------------------------------
def __truediv__( self, z ):
'self/z'
if isinstance( z, Punkt ):
ans = Punkt( self.z0/z.z0 )
else:
ans = Punkt( self.z0/z )
return ans
#--------------------------------------------------------------------------
def __pow__( self, z ):
'self**z'
if isinstance( z, Punkt ):
ans = Punkt( self.z0**z.z0 )
else:
ans = Punkt( self.z0**z )
return ans
#--------------------------------------------------------------------------
def __radd__( self, z ):
'z + self'
return self + z
#--------------------------------------------------------------------------
def __rsub__( self, z ):
'z - self'
return -( self - z )
#--------------------------------------------------------------------------
def __rmul__( self, z ):
'z*self'
return self*z
#--------------------------------------------------------------------------
def __rdiv__( self, z ):
'z/self'
return self.__rtruediv__( z )
#--------------------------------------------------------------------------
def __rtruediv__( self, z ):
'z/self'
if isinstance( z, Punkt ):
ans = Punkt( z.z0/self.z0 )
else:
ans = Punkt( z/self.z0 )
return ans
#--------------------------------------------------------------------------
#==============================================================================
class Matrix( Zahl ):
'''zwei komplexe zahlen als 2x2 matrix.
( z0.real z0.imag )
( z1.real z1.imag )
eigensystem nur fuer symmetrische matrizen moeglich.
'''
#--------------------------------------------------------------------------
def __init__( self, a, b=0 ):
if isinstance( a, Punkt ) or self.isNumber( a ):
self.z0 = Punkt( a )
self.z1 = Punkt( b )
elif isinstance( a, Matrix ): # kopie
self.z0 = Punkt( a.z0 )
self.z1 = Punkt( a.z1 ) # c+i*d
#--------------------------------------------------------------------------
def sym( self ):
'symmetrischer anteil'
s = 0.5*( self.z0.imag + self.z1.real )
return Matrix( Punkt( self.z0.real, s ), Punkt( s, self.z1.imag ) )
#--------------------------------------------------------------------------
def inv( self ):
'inverse'
d = i/self.det()
return Matrix( -d*self.z1, d*self.z0 ).tr()
#--------------------------------------------------------------------------
def mmul( self, M ):
'matrizenprodukt'
if isinstance( M, Matrix ):
# self.z0=a+i*b=r, self.z1=c+i*d=s, M.z0=e+i*g=u, M.z1=f+i*h=v
N = M.tr()
ans = Matrix( Punkt( self.z0.dot( N.z0 ),
self.z0.dot( N.z1 ) ),
Punkt( self.z1.dot( N.z0 ),
self.z1.dot( N.z1 ) ) )
elif isinstance( M, Punkt ):
ans = Punkt( self.z0.dot( M ), self.z1.dot( M ) )
else:
raise NotImplementedError( 'mmul fuer Matrix %s not implemented'%
type( M ) )
return ans
#--------------------------------------------------------------------------
def eig(self):
'''eigensystem ( ew, ev )
| a b |(u)=(a*u+b*v)=l*(u) -> b*v = (l-a)*u
| c d |(v) (c*u+d*v) (v) -> c*u = (l-d)*v
0 = l**2 - sp*l + det = l**2 - sp*l + 0.25*sp**2 - 0.25*sp**2 + det
( l - 0.5*sp )**2 = 0.25*sp**2 - det
'''
( sp, det ) = ( self.spur(), self.det() )
if 0.25*sp**2 < det:
raise ValueError( 'eig nur fuer symmetrische matrizen moeglich' )
w = cmath.sqrt( 0.25*sp**2 - det ).real
( l1, l2 ) = ew = ( 0.5*sp + w, 0.5*sp - w )
ev = []
( a, b, c, d ) = ( self.z0.real, self.z0.imag,
self.z1.real, self.z1.imag )
for l in ( l1, l2 ):
if abs( b ) > abs( c ):
ev.append( Punkt( b, l - a ).normalize() )
else:
ev.append( Punkt( l - d, c ).normalize() )
return ( ew, ev )
#--------------------------------------------------------------------------
def tr( self ):
'transponierte'
return Matrix( Punkt( self.z0.real, self.z1.real ),
Punkt( self.z0.imag, self.z1.imag ) )
#--------------------------------------------------------------------------
def spur( self ):
'spur'
return self.z0.real + self.z1.imag
#--------------------------------------------------------------------------
def det( self ):
'determinante'
return self.z0.cross( self.z1 )
#--------------------------------------------------------------------------
def toSquare( self ):
'gibt mf.Vec'
return mf.Square( [ self.z0.toVec(), self.z1.toVec() ] )
#--------------------------------------------------------------------------
def toNpTensor( self ):
'gibt nf.Vec'
return nf.Tensor( [ self.z0.toNpVec(), self.z1.toNpVec() ] )
#--------------------------------------------------------------------------
def __str__( self ):
'self.z0,b'
fmt = 'Matrix( [ [ %s,%s ], [ %s,%s ] ] )'%tuple( 4*[ floatfmt ] )
return fmt%( self.z0.real, self.z0.imag, self.z1.real, self.z1.imag )
#--------------------------------------------------------------------------
def __neg__( self ):
'-self'
return Matrix( -self.z0, -self.z1 )
#--------------------------------------------------------------------------
def __eq__( a,b ):
'a==b'
ans = 0
if isinstance( b, Matrix ):
ans = a.z0 == b.z0 and a.z1 == b.z1
return ans
#--------------------------------------------------------------------------
def __add__( a, b ):
'a + b'
return Matrix( a.z0 + b.z0, a.z1 + b.z1 )
#--------------------------------------------------------------------------
def __sub__( a, b ):
'a - b'
return Matrix( a.z0 - b.z0, a.z1 - b.z1 )
#--------------------------------------------------------------------------
def __mul__( a, b ):
'a*b'
return Matrix( b*a.z0, b*a.z1 )
#--------------------------------------------------------------------------
def __div__( a, b ):
'a/b'
return self*( 1./z )
#--------------------------------------------------------------------------
def __rmul__( a, b ):
'b*a'
return a*b
#--------------------------------------------------------------------------
#==============================================================================
class Eye( Matrix ):
'einheitsmatrix'
#--------------------------------------------------------------------------
def __init__( self ):
Matrix.__init__( self, Punkt( 1,0 ), Punkt( 0,1 ) )
#--------------------------------------------------------------------------
#==============================================================================
class Geometrie( Zahl ):
'geometrisches objekt mit einem referenzpunkt z0'
#--------------------------------------------------------------------------
def __init__( self, z0 ):
self.z0 = Punkt( z0 )
#--------------------------------------------------------------------------
#==============================================================================
class Gerade( Geometrie ):
'berechnungen mit geraden'
#--------------------------------------------------------------------------
def __init__( self, a, b, *punkte ):
if len( punkte ) > 0:
g = self.regression( [ a, b ] + list( punkte ) )
( a,b ) = ( g.z0, g.z1 )
self.e = g.e
else:
self.e = ( a - b )/abs( a - b )
Geometrie.__init__( self, a )
self.z1 = Punkt( b )
#--------------------------------------------------------------------------
def copy( g ):
'gibt eine kopie dieser geraden'
return Gerade( self.z0, self.z1 )
#--------------------------------------------------------------------------
def regression( self, p ):
'regression durch mehr als 2 punkte in 2d'
# minimiert sum( [ |e.cross( pj - s )|**2 for pj in p ] )
# siehe /home/austausch/doku/wissen/mathe/regression
o = Punkt( 0 )
O = o.outer( o )
p = [ Punkt( pj ) for pj in p ]
N = len( p )
sump = sum( p, o )
sumd = sum( [ pj.dot( pj ) for pj in p ] )
sumo = sum( [ pj.outer( pj ) for pj in p ], O )
s = sump/N
T = sumd*Eye() - sumo + sump.outer( s )
( ew, ev ) = T.eig()
n = ev[ 0 ]
if ew[ 1 ] < ew[ 0 ]:
n = ev[ 1 ]
return Gerade( s, s+n )
#--------------------------------------------------------------------------
def trimm( g, a=0, b='' ):
'trimmt diese gerade auf den abschnitt z0-a:z0+b'
if isinstance( b, str ):
b = a
( g.z0, g.z1 ) = ( g.z0 - a*g.e, g.z0 + b*g.e )
#--------------------------------------------------------------------------
def fussPunkt( g, p ):
"""projiziert p auf g
p = g.z0 + z*g.e -> z = ( p - g.z0 )/g.e = a + i*b -> p' = g.z0 + a*g.e
"""
z = ( p - g.z0 )/g.e
return Punkt( g.z0 + z.real*g.e )
#--------------------------------------------------------------------------
def schnittPunkt( a, b ):
'gibt den schnittpunkt der geraden a und b'
ans = []
q = a.e.cross( b.e )
if q != 0: # nicht parallel
ans = [ Punkt( a.z0 + b.e.cross( a.z0 - b.z0 )/q*a.e ) ]
return ans
#--------------------------------------------------------------------------
def dist( g, z=0 ):
'Berechnet den Abstand von z zu dieser Geraden'
return abs( g.e.cross( z - g.z0 ) )
#--------------------------------------------------------------------------
def angle( a, b ):
'gibt den winkel zwischen diesen geraden'
return a.e.angle( b.e )
#--------------------------------------------------------------------------
def winkelhalbierenden( self, b ):
'gibt die winkelhalbierenden der geraden a und b'
z = self.schnittPunkt( b )[ 0 ]
g = Gerade( z, Punkt( z + self.e + b.e ) )
return ( g, g.senkrechteDurch( z ) )
#--------------------------------------------------------------------------
def paralleleDurch( self, p ):
'gibt die parallele zu dieser geraden durch p'
return Gerade( p, p + self.e )
#--------------------------------------------------------------------------
def senkrechteDurch( self, p ):
'gibt die senkrechte zu dieser geraden durch p'
return Gerade( p, p + i*self.e )
#--------------------------------------------------------------------------
def spiegel( g, obj ):
'geradenspiegelung an dieser geraden'
if g.isNumber( obj ):
obj = Punkt( obj )
if isinstance( obj, Punkt ):
# obj = g.z0 + a*g.e + b*i*g.e = g.z0 + (a+i*b)*g.e->
ans = g.z0 + ( obj - g.z0 ).conjugate()*g.e**2
elif isinstance( obj, Gerade ):
ans = Gerade( g.spiegel( obj.z0 ), g.spiegel( obj.z1 ) )
elif isinstance( obj, Kreis ):
ans = Kreis( g.spiegel( obj.z0 ), obj.r )
elif isinstance( obj, Dreieck ):
ans = Dreieck( g.spiegel( obj.z0 ), g.spiegel( obj.z1 ),
g.spiegel( obj.z2 ) )
elif hasattr( obj, '__len__' ): # liste
ans = [ g.spiegel( p ) for p in obj ]
else:
raise NotImplementedError( 'Geradenspiegelung fuer %s'%type( obj )
+ 'nicht implementiert' )
#print( 'spiegel %s -> %s'%( obj, ans ) )
return ans
#--------------------------------------------------------------------------
def __str__( self ):
fmt = 'Gerade( ( %s,%s ), ( %s,%s ) )'%tuple( 4*[ floatfmt ] )
return fmt%( self.z0.real, self.z0.imag,
( self.z0+self.e ).real, ( self.z0+self.e ).imag )
#--------------------------------------------------------------------------
#==============================================================================
class Dreieck( Geometrie ):
'Berechnungen eines Dreiecks'
#--------------------------------------------------------------------------
def __init__( self, z0, z1, z2 ):
Geometrie.__init__( self, z0 )
self.z1 = Punkt( z1 )
self.z2 = Punkt( z2 )
#--------------------------------------------------------------------------
def inhalt( d ):
'gibt den flaecheninhalt'
return 0.5*abs( ( d.z1 - d.z0 ).cross( d.z2 - d.z1 ) )
#--------------------------------------------------------------------------
def umkreis( self ):
'gibt den Umkreis dieses Dreiecks'
( A, B, C ) = ( self.z0, self.z1, self.z2 )
( a,b,c ) = [ X.dot( X ) for X in ( A, B, C ) ]
d = 2*( A.cross( B ) + B.cross( C ) + C.cross( A ) )
m = i/d*( ( b - c )*A + ( c - a )*B + ( a - b )*C )
r = abs( A - m )
return Kreis( m,r )
#--------------------------------------------------------------------------
def inkreis( self ):
'gibt den Inkreis dieses Dreiecks'
( A, B, C ) = ( self.z0, self.z1, self.z2 )
a = abs( B - C )
b = abs( C - A )
c = abs( A - B )
s = a + b + c
m = ( a*A + b*B + c*C )/s
s *= 0.5
r = cmath.sqrt( ( s - a )*( s - b )*( s - c )/s )
return Kreis( m,r )
#--------------------------------------------------------------------------
def schwerpunkt( self ):
'gibt den schwerpunkt dieses dreiecks'
return ( self.z0 + self.z1 + self.z2 )/3.
#--------------------------------------------------------------------------
def hoehenschnittpunkt( self ):
'gibt den hoehenschnittpunkt dieses dreiecks'
( A,B,C ) = ( self.z0, self.z1, self.z2 )
d = A.cross( B ) + B.cross( C ) + C.cross( A )
return i/d*( A.dot( B - C )*A + B.dot( C - A )*B + C.dot( A - B )*C )
#--------------------------------------------------------------------------
def baryPunkt( self, m1, m2, m3 ):
'gibt den punkt mit baryzentrischen Koordinaten m1, m2, m3'
return ( m1*self.z0 + m2*self.z1 + m3*self.z2 )/( m1 + m2 + m3 )
#--------------------------------------------------------------------------
def baryKoord( self, z ):
'gibt die normierten baryzentrischen Koordinaten von z'
d = ( self.z1 - self.z0 ).cross( self.z2 - self.z1 )
m1 = ( self.z1 - z ).cross( self.z2 - z )/d
m2 = ( self.z2 - z ).cross( self.z0 - z )/d
return ( m1, m2, 1 - m1 - m2 ) # summe = 1
#--------------------------------------------------------------------------
def trilinPunkt( self, t1, t2, t3 ):
'gibt den punkt mit trilinearen Koordinaten t1, t2, t3'
( a,b,c ) = ( abs( self.z1-self.z2 ), abs( self.z2-self.z0 ),
abs( self.z0-self.z1 ) )
( m1, m2, m3 ) = ( a*t1, b*t2, c*t3 )
return ( m1*self.z0 + m2*self.z1 + m3*self.z2 )/( m1 + m2 + m3 )
#--------------------------------------------------------------------------
def trilinKoord( self, z ):
'gibt die normierten trilinearen koordinaten von z'
( a,b,c ) = ( abs( self.z1-self.z2 ), abs( self.z2-self.z0 ),
abs( self.z0-self.z1 ) )
dBC = ( self.z2 - self.z1 ).cross( z - self.z1 )/a
dCA = ( self.z0 - self.z2 ).cross( z - self.z2 )/b
dAB = ( self.z1 - self.z0 ).cross( z - self.z0 )/c
d = dAB + dBC + dCA
return ( dBC/d, dCA/d, dAB/d ) # summe = 1
#--------------------------------------------------------------------------
def enthaelt( self, z ):
'prueft, ob z in diesem dreieck liegt'
return all( [ x>=0 for x in self.baryKoord( z ) ] )
#--------------------------------------------------------------------------
#==============================================================================
class Polygon:
#--------------------------------------------------------------------------
def __init__( self, *points ):
self.points = points
self.nump = len( points )
rev = list( reversed( points ) )
if ( abs( self.innenwinkelsumme( points ) )
> abs( self.innenwinkelsumme( rev ) ) ):
self.points = rev
#--------------------------------------------------------------------------
def innenwinkelsumme( self, points=[] ):
if len( points ) == 0:
points = self.points
( z0, z1 ) = points[ -2: ]
w = 0.
for z2 in points:
w += ( z2 - z1 ).angle( z0 - z1 )%( 2*cmath.pi )
z0 = z1
z1 = z2
return w
#--------------------------------------------------------------------------
def inhalt( self ):
'gibt den flaecheninhalt dieses polygons'
z0 = self.points[ -1 ]
A = 0.
for z1 in self.points:
A += z0.cross( z1 )
z0 = z1
return 0.5*A
#--------------------------------------------------------------------------
def schwerpunkt( self, k=0 ):
'findet schwerpunkt fuer konvexe polygone'
ps = self.points[ k: ] + self.points[ :k ]
( p0, p1 ) = ps[ :2 ]
s = 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 Kreis( Geometrie ):
'Berechnungen eines Kreis'
#--------------------------------------------------------------------------
def __init__( self, m, r ):
'kreis um m mit radius |r|'
Geometrie.__init__( self, m )
self.r = abs( r )
#--------------------------------------------------------------------------
def inhalt( k ):
'gibt den flaecheninhalt'
return cmath.pi*self.r**2
#--------------------------------------------------------------------------
def project( k, p ):
'projiziert p auf k'
return k.z0 + k.r/abs( p - k.z0 )*( p - k.z0 )
#--------------------------------------------------------------------------
def schnittPunkteGerade( k, g ):
'gibt 0,1 oder 2 Schnittpunkte der Geraden g mit diesem Kreis k.'
ans = ()
d = g.z0 - k.z0
p = g.e.dot( d )
q = p**2 - abs( d )**2 + k.r**2
if q > 0:
w = cmath.sqrt( q )
ans = ( g.z0 - ( w + p )*g.e, g.z0 + ( w - p )*g.e )
elif abs( q ) < 1.e-9*k.r: # w=0
ans = ( g.z0 - p*g.e, )
return [ Punkt( p ) for p in ans ]
#--------------------------------------------------------------------------
def schnittPunkteKreis( k1, k2 ):
'gibt 0,1 oder 2 schnittpunkte von Kreis k1 und k2'
ans = ()
g = k2.z0 - k1.z0
ga = abs( g )
sx = 0.5*( ( k1.r**2 - k2.r**2 )/ga**2 + 1 )
sy = k1.r**2/ga**2 - sx**2
if sy > 0:
sy = cmath.sqrt( sy ).real
s = complex( sx,sy )
sq = s.conjugate()
ans = ( ( k1.z0 + s*g ), ( k1.z0 + sq*g ) )
elif abs( sy ) < 1.e-9*k1.r: # sy=0
ans = ( ( k1.z0 + sx*g ), )
return [ Punkt( p ) for p in ans ]
#--------------------------------------------------------------------------
def tangenteIn( self, P ):
'gibt die tangente in P'
return Gerade( P, P + i*( P - self.z0 ) )
#--------------------------------------------------------------------------
def tangentenDurch( self, P ):
'gibt 0,1 oder 2 tangenten an diesen kreis durch p'
ans = []
d = abs( P - self.z0 )
if d > self.r:
k = Kreis( 0.5*( self.z0 + P ), 0.5*abs( self.z0 - P ) )
( T1, T2 ) = self.schnittPunkteKreis( k )
ans = [ Gerade( T1, T1 + i*( T1 - self.z0 ) ),
Gerade( T2, T2 + i*( T2 - self.z0 ) ) ]
elif d == self.r:
ans = [ Gerade( P, P + i*( P - self.z0 ) ) ]
return ans
#--------------------------------------------------------------------------
def tangentenAussen( K,k ):
'''die aeusseren tangenten mit kreis k
ges e: M+R*e+l*i*e-r*e=M+(R-r+i*l)*e=m, l**2+(R-r)**2=|M-m|**2
oder M-R*e+l*i*e+r*e=M+(r-R+i*l)*e=m
'''
( M, R, m, r ) = ( K.z0, K.r, k.z0, k.r )
l = cmath.sqrt( ( M - m ).dot( M - m ) - ( R - r )**2 )
e1 = ( m - M )/complex( R - r, l )
t1 = Gerade( M + R*e1, m + r*e1 )
e2 = ( m - M )/complex( r - R, l )
t2 = Gerade( M - R*e2, m - r*e2 )
return ( t1, t2 )
#--------------------------------------------------------------------------
def tangentenInnen( K,k ):
'''die inneren tangenten mit kreis k
ges e: M+R*e+l*i*e+r*e=M+(R+r+i*l)*e=m, l**2+(R+r)**2=|M-m|**2
oder M+R*e-l*i*e+r*e=M+(r+R-i*l)*e=m
'''
( M, R, m, r ) = ( K.z0, K.r, k.z0, k.r )
ans = []
l = ( M - m ).dot( M - m ) - ( R + r )**2
if l >= 0:
l = cmath.sqrt( l )
e1 = ( m - M )/complex( R + r, l )
t1 = Gerade( M + R*e1, m - r*e1 )
e2 = ( m - M )/complex( R + r, -l )
t2 = Gerade( M + R*e2, m - r*e2 )
ans = ( t1, t2 )
cumtest( cum, 'e12', ( 1,1 ), ( abs( e1 ), abs( e2 ) ), fp )
return ans
#--------------------------------------------------------------------------
def optKreis( self, P ):
'minimiert sum( [ ( ( Pj - M ).dot( Pj - M ) - r**2 )**2 ] )'
N = len( P )
o = Punkt( 0,0 )
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 = 2*sum( [ Pj.outer( Pj ) for Pj in P ], a.outer( -0.5*z ) )
M = A.inv().mmul( c - b*a )
r = cmath.sqrt( b + M.dot( M - z ) )
return Kreis( M, r )
#--------------------------------------------------------------------------
def spiegel( self, obj, tol=1.e-6 ):
'kreisspiegelung an diesem kreis'
if self.isNumber( obj ):
obj = Punkt( obj )
if isinstance( obj, Punkt ):
ans = self.z0 + self.r**2/( obj - self.z0 ).conjugate()
elif isinstance( obj, Gerade ):
if abs( obj.dist( self.z0 ) ) < tol*self.r: # gerade durch O
ans = obj.copy( )
else: # abb auf kreis durch ursprung
p = obj.fussPunkt( self.z0 ) # C'
ps = self.z0 + self.r**2/( p - self.z0 ).conjugate() # C
ans = Kreis( 0.5*( self.z0 + ps ), 0.5*abs( self.z0 - ps ) )
elif isinstance( obj, Kreis ):
g = obj.z0 - self.z0
if abs( abs( g ) - obj.r ) < tol*self.r: # durch O -> Gerade
p = self.z0 + 2*g
ps = self.z0 + self.r**2/( p - self.z0 ).conjugate()
ans = Gerade( ps, ps + i*g )
else: # abb auf kreis
s = self.r**2/( abs( g )**2 - obj.r**2 )
ans = Kreis( self.z0 + s*g, abs( s )*obj.r )
elif hasattr( obj, '__len__' ): # liste
ans = [ self.spiegel( p ) for p in obj ]
else:
raise NotImplementedError( 'Kreisspiegelung fuer %s'%type( obj )
+ 'nicht implementiert' )
#print( 'spiegel %s -> %s'%( obj, ans ) )
return ans
#--------------------------------------------------------------------------
def __str__( self ):
fmt = 'Kreis( ( %s,%s ), %s )'%tuple( 3*[ floatfmt ] )
return fmt%( self.z0.real, self.z0.imag, self.r )
#--------------------------------------------------------------------------
#==============================================================================
class Test:
'testet routinen'
#--------------------------------------------------------------------------
def __init__( self ):
global uniform, mf, nf, cumtest
from random import uniform
import matfunc2 as mf, numpyfunc as nf
cumtest = nf.cumtest
# self.punkt()
# self.angle()
# self.schnittPunkt()
# self.dist()
# self.winkelhalbierende()
self.dreieck()
# self.kreis()
# self.matrix()
# self.spiegelung()
print( 'Test: %d/%d richtig'%tuple( cum ) )
cumtest( cum, '', 0,0, fout=fp )
#--------------------------------------------------------------------------
def punkt( self ):
( p,q ) = [ Punkt( u(),u() ) for j in 'ab' ]
( j,r,c ) = ( 4, 4., 4.+0j )
cumtest( cum, '+', Punkt( p.z0 + q.z0 ).p2t(), ( p + q ).p2t(),
fout=fp )
cumtest( cum, '-', Punkt( p.z0 - q.z0 ).p2t(), ( p - q ).p2t(),
fout=fp )
cumtest( cum, '*', Punkt( p.z0 * q.z0 ).p2t(), ( p * q ).p2t(),
fout=fp )
cumtest( cum, '/', Punkt( p.z0 / q.z0 ).p2t(), ( p / q ).p2t(),
fout=fp )
for z in ( c,r,j ):
cumtest( cum, 'r+', Punkt( z + p.z0 ).p2t(), ( z + p ).p2t(),
fout=fp )
cumtest( cum, 'l+', Punkt( p.z0 + z ).p2t(), ( p + z ).p2t(),
fout=fp )
cumtest( cum, 'r-', Punkt( z - p.z0 ).p2t(), ( z - p ).p2t(),
fout=fp )
cumtest( cum, 'l-', Punkt( p.z0 - z ).p2t(), ( p - z ).p2t(),
fout=fp )
cumtest( cum, 'r*', Punkt( z * p.z0 ).p2t(), ( z * p ).p2t(),
fout=fp )
cumtest( cum, 'l*', Punkt( p.z0 * z ).p2t(), ( p * z ).p2t(),
fout=fp )
cumtest( cum, 'r/?', Punkt( z / p.z0 ).p2t(), ( z / p ).p2t(),
fout=fp )
cumtest( cum, 'l/', Punkt( p.z0 / z ).p2t(), ( p / z ).p2t(),
fout=fp )
#--------------------------------------------------------------------------
def matrix( self ):
( a,b,c,d,e,f,g,h,w,x,y,z ) = [ u() for j in 12*'a' ]
mz = Matrix( Punkt( a, b ), # r
Punkt( c, d ) ) # s
cumtest( cum, 'objs', 11*( True, ), (
isinstance( complex( c,d )*Punkt( a,b ), Punkt ),
isinstance( Punkt( a,b )*complex( c,d ), Punkt ),
isinstance( Punkt( a,b ) + Punkt( c,d ), Punkt ),
isinstance( Punkt( a,b ).outer( Punkt( c,d ) ), Matrix ),
isinstance( complex( c,d )*mz, Matrix ),
isinstance( mz*complex( c,d ), Matrix ),
isinstance( mz.mmul( mz ), Matrix ),
isinstance( ( complex( c,d )*mz ).z0, Punkt ),
isinstance( ( mz*complex( c,d ) ).z0, Punkt ),
isinstance( ( complex( c,d )*mz ).z1, Punkt ),
isinstance( ( mz*complex( c,d ) ).z1, Punkt )
), fout=fp )
mm = mz.toSquare()
cumtest( cum, 'mz2mm', mf.Square( [ [ a,b ],[ c,d ] ] ), mm, fout=fp )
cumtest( cum, 'mt2mt', mz.tr().toSquare(), mm.tr(), fout=fp )
nz = Matrix( Punkt( e, f ), # u
Punkt( g, h ) ) # v
nm = mf.Square( [ [ e,f ],[ g,h ] ] )
cumtest( cum, 'nz2nm', nz.toSquare(), nm, fout=fp )
vm = mf.Vec( [ x,y ] )
cumtest( cum, 'mm.vm', ( a*x + b*y, c*x + d*y ), mm.mmul( vm ),
fout=fp )
cumtest( cum, 'mm.nm', ( [ b*g + a*e, b*h + a*f],
[ d*g + c*e, d*h + c*f ] ), mm.mmul( nm ),
fout=fp )
vz = Punkt( x,y )
cumtest( cum, 'a o bmf', 2*( ( ( a*c, a*d ), ( b*c, b*d ) ), ), (
mf.Vec( [ a,b ] ).outer( mf.Vec( [ c,d ] ) ),
Punkt( a,b ).outer( Punkt( c,d ) ).toSquare() ), fout=fp )
#nz = Square(
cumtest( cum, 'mz.vz', Punkt( a*x + b*y, c*x + d*y ).z0,
mz.mmul( vz ).z0, fout=fp )
cumtest( cum, 'mz.nz', mz.mmul( nz ).toSquare(), mm.mmul( nm ),
fout=fp )
cumtest( cum, 'spur', mm.spur(), mz.spur(), fout=fp )
cumtest( cum, 'det', mm.det(), mz.det(), fout=fp )
mz1 = mz.inv()
e = mz.mmul( mz1 )
cumtest( cum, 'inv', Eye().toSquare(), e.toSquare(), fout=fp )
c = Punkt( u(), 0 )
l = c #Punkt( c )
dum1 = ( mz - Eye()*l ).toSquare()
dum2 = mm - c.z0*mf.eye( 3 )
cumtest( cum, 'm-lE', dum1, dum2,
fout=fp )
cumtest( cum, 'char', ( mz - Eye()*l ).det(),
Punkt( l**2 - mz.spur()*l + mz.det() ).z0, fout=fp )
mz = mz.sym()
mm = mz.toSquare()
( ( l1,l2 ), ( v1,v2 ) ) = mz.eig()
cumtest( cum, 'det1', 0, ( mz - l1*Eye() ).det(), fout=fp )
cumtest( cum, 'det2', 0,
( mz - l2*Matrix( Punkt( 1,0 ), Punkt( 0,1 ) ) ).det(),
fout=fp )
cumtest( cum, 'l1', Punkt( l1*v1 ).z0, mz.mmul( v1 ).z0, fout=fp )
cumtest( cum, 'l2', Punkt( l2*v2 ).z0, mz.mmul( v2 ).z0, fout=fp )
#--------------------------------------------------------------------------
def kreis( self ):
m1 = Punkt( u(), u() )
( r1, r2 ) = [ abs( u() ), abs( u() ) ]
g = Punkt( u(), u() ).normalize()
m2 = m1 + cmath.sqrt( r1**2 + r2**2 )*g
k1 = Kreis( m1, r1 )
k2 = Kreis( m2, r2 )
p = Punkt( u(), u() )
p = k1.z0 + u()*k1.r*p/abs( p )
for z in k1.schnittPunkteGerade( Gerade( p, Punkt( u(), u() ) ) ):
cumtest( cum, 'k^g', k1.r, abs( z - k1.z0 ), fout=fp )
( z1, z2 ) = k1.schnittPunkteKreis( k2 )
m = Punkt( 0.5*( z1 + z2 ) )
for z in ( z1, z2 ):
cumtest( cum, 'k^k1', k1.r, abs( z - k1.z0 ), fout=fp )
cumtest( cum, 'k^k2', k2.r, abs( z - k2.z0 ), fout=fp )
cumtest( cum, 'm', k1.r**2, m.dist( z )**2 + m.dist( k1.z0 )**2,
fout=fp )
cumtest( cum, 'm', k2.r**2, m.dist( z )**2 + m.dist( k2.z0 )**2,
fout=fp )
P = [ Punkt( u(), u() ) for j in 10*'a' ]
k = k1.optKreis( P )
km = nf.Kreis().optKreis( [ Pj.toNpVec() for Pj in P ] )
cumtest( cum, 'opt', ( k.z0.real, k.z0.imag, k.r ),
( km.x0[ 0 ], km.x0[ 1 ], km.r ), fout=fp )
( P0, P1 ) = P[ :2 ]
( v0, v1 ) = ( P0.toNpVec(), P1.toNpVec() )
g = Gerade( P0, P1 )
gv = nf.Gerade( v0, v1 )
k0 = Kreis( P0, 3. )
k1 = Kreis( P1, 3. )
kn0 = nf.Kreis( v0, 3. )
kn1 = nf.Kreis( v1, 3. )
cumtest( cum, 'spP', kn0.spiegel( v1 ), k0.spiegel( P1 ).p2t(),
fout=fp )
cumtest( cum, 'spK', kn0.spiegel( kn1 ).x0, k0.spiegel( k1 ).z0.p2t(),
fout=fp )
cumtest( cum, 'spG', kn0.spiegel( gv ).x0, k0.spiegel( g ).z0.p2t(),
fout=fp )
( M, R, m, r ) = ( complex( u(),u() ), abs( u() ),
complex( u(),u() ), abs( u() ) )
K = Kreis( M,R )
k = Kreis( m,r )
( t1, t2 ) = K.tangentenAussen( k )
cumtest( cum, 'tout', ( R, R, r, r ), (
t1.dist( M ), t2.dist( M ), t1.dist( m ), t2.dist( m ) ), fp )
while abs( m - M ) < R + r:
m += m - M
k = Kreis( m, r )
( t1, t2 ) = K.tangentenInnen( k )
cumtest( cum, 'tin', ( R, R, r, r ), (
t1.dist( M ), t2.dist( M ), t1.dist( m ), t2.dist( m ) ), fp )
#--------------------------------------------------------------------------
def dreieck( self ):
e = ( za, zb, zc ) = [ Punkt( u(), u() ) for j in 'abc' ]
( va, vb, vc ) = [ p.toNpVec() for p in ( za, zb, zc ) ]
d = Dreieck( za, zb, zc )
( a,b,c ) = ( abs( zb - zc ), abs( za - zc ), abs( za - zb ) )
cumtest( cum, 'inhalt', 3*( d.inhalt(), ), (
0.5*abs( ( za-zc ).cross( zb-za ) ),
0.5*abs( ( zb-za ).cross( zc-zb ) ),
0.5*abs( ( zc-zb ).cross( za-zc ) ) ), fp )
ku = d.umkreis()
cumtest( cum, 'umkreis', 3*( ku.r, ), [ abs( z - ku.z0 ) for z in e ],
fp )
ki = d.inkreis()
cumtest( cum, 'inkreis', 3*( ki.r, ), (
Gerade( za, zb ).dist( ki.z0 ), Gerade( zb, zc ).dist( ki.z0 ),
Gerade( zc, za ).dist( ki.z0 ) ), fp )
S = d.schwerpunkt()
cumtest( cum, 'S', ( 0,0,0 ), (
Gerade( za,0.5*(zb + zc ) ).dist( S ),
Gerade( zb,0.5*(zc + za ) ).dist( S ),
Gerade( zc,0.5*(za + zb ) ).dist( S ) ), fp )
H = d.hoehenschnittpunkt()
cumtest( cum, 'H', ( 0,0,0 ), (
Gerade( zb,zc ).senkrechteDurch( za ).dist( H ),
Gerade( zc,za ).senkrechteDurch( zb ).dist( H ),
Gerade( za,zb ).senkrechteDurch( zc ).dist( H ) ), fp )
cumtest( cum, 'barya', ( 1,0,0 ), d.baryKoord( za ), fp )
cumtest( cum, 'baryb', ( 0,1,0 ), d.baryKoord( zb ), fp )
cumtest( cum, 'baryc', ( 0,0,1 ), d.baryKoord( zc ), fp )
z = Punkt( u(), u() )
( ba, bb, bc ) = bcoord = d.baryKoord( z )
cumtest( cum, 'bary', z.p2t(), d.baryPunkt( *bcoord ).p2t(), fp )
cumtest( cum, 'trilina', ( 1,0,0 ), d.trilinKoord( za ), fp )
cumtest( cum, 'trilinb', ( 0,1,0 ), d.trilinKoord( zb ), fp )
cumtest( cum, 'trilinc', ( 0,0,1 ), d.trilinKoord( zc ), fp )
( ta, tb, tc ) = tcoord = d.trilinKoord( z )
cumtest( cum, 'tribary', 2*( ta*a/ba, ), ( tb*b/bb, tc*c/bc ), fp )
cumtest( cum, 'trilin', z.p2t(), d.trilinPunkt( *tcoord ).p2t(), fp )
cumtest( cum, 'norm', ( 1,1 ), ( sum( bcoord ), sum( tcoord ) ), fp )
P = 0.5*(za + zb )
g = P - zc
print( d.trilinKoord( S ), d.trilinKoord( H ),
all( [ x>=0 for x in d.baryKoord( P-1.001*g ) ] ) )
#--------------------------------------------------------------------------
def angle( self ):
( za, zb ) = [ Punkt( u(), u() ) for j in 'ab' ]
( va, vb ) = [ p.toNpVec() for p in ( za, zb ) ]
a = za.angle( zb )
b = va.angle( vb )
cumtest( cum, 'angle', 0, min( abs( cmath.pi - a - b ),
abs( a - b ) ), fout=fp )
#--------------------------------------------------------------------------
def schnittPunkt( self ):
( za, zb, zc, zd ) = [ Punkt( u(), u() ) for j in 'abcd' ]
( va, vb, vc, vd ) = [ p.toNpVec() for p in ( za, zb, zc, zd ) ]
zab = Gerade( za, zb )
zcd = Gerade( zc, zd )
zp = zab.schnittPunkt( zcd )[ 0 ]
vp = nf.Gerade( va, vb ).schnittPunkt( nf.Gerade( vc, vd ) )
cumtest( cum, 'spkt', zp.toNpVec(), vp, fout=0 )
#--------------------------------------------------------------------------
def dist( self ):
( za, zb, z ) = [ Punkt( u(), u() ) for j in 'abc' ]
( va, vb, v ) = [ p.toNpVec() for p in ( za, zb, z ) ]
cumtest( cum, 'distP', za.dist( zb ), va.dist( vb ), fout=fp )
cumtest( cum, 'dist', abs( nf.Gerade( va, vb ).dist( v ) ),
abs( Gerade( za, zb ).dist( z ) ), fout=fp )
#--------------------------------------------------------------------------
def winkelhalbierende( self ):
( za, zb, zc, zd ) = [ Punkt( u(), u() ) for j in 'abcd' ]
( va, vb, vc, vd ) = [ p.toNpVec() for p in ( za, zb, zc, zd ) ]
zab = Gerade( za, zb )
zcd = Gerade( zc, zd )
zp = zab.winkelhalbierende( zcd )
vp = nf.Gerade( va, vb ).winkelhalbierende( nf.Gerade( vc, vd ) )
cumtest( cum, 'whz', zp.z0.toNpVec(), vp.x0, fout=fp )
cumtest( cum, 'whg', zp.e.toNpVec(), vp.g, fout=fp )
#--------------------------------------------------------------------------
def spiegelung( self ):
( za, zb, zc, zd ) = [ Punkt( u(), u() ) for j in 'abcd' ]
# punkt
zb2 = za.spiegel( zb )
cumtest( cum, 'punkt', 1*( 0, ), ( za.dist( 0.5*( zb + zb2 ) ),
), fout=fp )
# gerade y=a*x+b
( a, b ) = ( zd.real, zd.imag )
g = Gerade( 0,1 )
h = Gerade( i*b, complex( 1,a+b ) )
h1 = g.spiegel( h )
hs = Gerade( -i*b, complex( 1,-a-b ) )
cumtest( cum, 'waagerecht', 2*( 0, ), (
hs.dist( h1.z0 ), hs.dist( h1.z1 ) ), fout=fp )
g = Gerade( 0,i )
h1 = g.spiegel( h )
hs = Gerade( i*b, complex( -1,a+b ) )
cumtest( cum, 'senkrecht', 2*( 0, ), (
hs.dist( h1.z0 ), hs.dist( h1.z1 ) ), fout=fp )
zab = Gerade( za, zb )
zc2 = zab.spiegel( zc )
zcc = Gerade( zc, zc2 )
cumtest( cum, 'gerade', 5*( 0, ), (
zab.dist( 0.5*( zc + zc2 ) ), zab.e.dot( zc2 - zc ),
zab.e.dot( zcc.e ), hs.dist( h1.z0 ), hs.dist( h1.z1 )
), fout=fp )
#--------------------------------------------------------------------------
#==============================================================================
if __name__ == '__main__':
Test()
