diff --git a/hypertiling/core.py b/hypertiling/core.py index f8f9455d3da5acdad31792595c4b275c5307eb7f..b7443b6e5c87863ed37b3f04f083aa469756033a 100644 --- a/hypertiling/core.py +++ b/hypertiling/core.py @@ -48,7 +48,8 @@ class HyperbolicTiling: polygon = HyperPolygon(self.p, self.q) for i in range(self.p): - z = complex(r * np.cos(i*self.phi), r * np.sin(i*self.phi)) # = r*exp(i*phi) + z = complex(r, i*self.phi) +# z = complex(r * np.cos(i*self.phi), r * np.sin(i*self.phi)) # = r*exp(i*phi) polygon.verticesP[i] = z polygon.verticesW[:, i] = p2w(z) return polygon @@ -88,7 +89,7 @@ class HyperbolicTiling: # compute center and angle center = np.round(adj_pgon.centerP, self.dgts) adj_pgon.find_angle(360) # divide the disk into 360*7 sectors - + # cut away cells outside the allowed sectors if adj_pgon.sector in sectors: # try adding to centerlist; it is a set() and takes care of duplicates @@ -128,7 +129,7 @@ class HyperbolicTiling: for p in range(1, k): for polygon in polygons: pgon = copy.deepcopy(polygon) - pgon.rotate(p*self.phi) + pgon.moeb_rotate(-p*self.phi) pgon.find_angle(1) # set polygon.sector such that the disk is divided into 1*p sectors self.polygons.append(pgon) diff --git a/hypertiling/hyperpolygon.py b/hypertiling/hyperpolygon.py index 307ffead702691022aac972a62f75f787a13aa10..d0741ec0362661b80ebed34bafffc6fd7d9d8f94 100644 --- a/hypertiling/hyperpolygon.py +++ b/hypertiling/hyperpolygon.py @@ -81,11 +81,11 @@ class HyperPolygon: self.verticesW[:, i] = p2w(self.verticesP[i]) - def moeb_inverse(self, z0, phi=0, s=0): - self.centerP = moeb_inverse_trafo(self.centerP, z0, -phi, s) + def moeb_inverse(self, z0): + self.centerP = moeb_origin_trafo_inverse(z0, self.centerP) self.centerW = p2w(self.centerP) for i in range(self.p): - z = moeb_inverse_trafo(self.verticesP[i], z0, -phi, s) + z = moeb_origin_trafo_inverse(z0, self.verticesP[i]) self.verticesP[i] = z self.verticesW[:, i] = p2w(self.verticesP[i]) @@ -120,7 +120,10 @@ class HyperPolygon: def find_angle(self, k, offset=0): # has to be called after the inverse trafo... - self.angle = np.arctan2(self.centerP.imag, self.centerP.real)*180/np.pi-offset # requires y first + angle = self.centerP.imag + angle += 2*np.pi if angle < 0 else 0 + angle -= 2*np.pi if angle > 2*np.pi else 0 + self.angle = angle * 180/np.pi - offset self.angle = np.round(self.angle, 10) self.angle += 360 if self.angle < 0 else 0 self.sector = floor(self.angle/(360/(k*self.p))) # k*p sectors; insert offset +0.1 here?! diff --git a/hypertiling/plot.py b/hypertiling/plot.py index 54bf05b03636cc9fe2e5fa18332e0fadab3a15b3..50730a9af0e3f03466dba7bcbcb516c27c37bd91 100644 --- a/hypertiling/plot.py +++ b/hypertiling/plot.py @@ -42,7 +42,8 @@ def poly2patch(polygons, colors=None, **kwargs): # extract vertex coordinates u = poly.verticesP # transform to matplotlib Polygon format - stack = np.column_stack((u.real,u.imag)) + stack = np.column_stack((u.real*np.cos(u.imag),u.real*np.sin(u.imag))) +# stack = np.column_stack((u.real,u.imag)) polygon = Polygon(stack, True) patches.append(polygon) diff --git a/hypertiling/transformation.py b/hypertiling/transformation.py index 68fa091944132d76fbbd1edc932e2865f28b30f2..a5ce9ec26dd921427528584ab17741e12fa026e2 100644 --- a/hypertiling/transformation.py +++ b/hypertiling/transformation.py @@ -1,32 +1,94 @@ import numpy as np +import math - +def myclamp(z): + angle = z.imag + angle += 2*np.pi if angle < -np.pi else 0 + angle -= 2*np.pi if angle > np.pi else 0 + return complex(z.real, angle) + def p2w(z): - x, y = z.real, z.imag - xx = x*x - yy = y*y - factor = 1 / (1-xx-yy) - return factor*np.array([(1+xx+yy), 2*x, 2*y]) + factor = 1 / (1-z.real*z.real) + return factor*np.array([(1+z.real*z.real), 2*z.real*np.cos(z.imag), 2*z.real*np.sin(z.imag)]) def w2p(point): [t, x, y] = point factor = 1 / (1+t) - return complex(x*factor, y*factor) + mag = np.hypot(x, y) + return complex(mag/(1+t), 2*np.arctan(y/(mag + x)) ) # maps all points z such that z0 -> 0, respecting the Poincare projection +# z = (|z|, phi) def moeb_origin_trafo(z0, z): - return (z-z0) / (1-z*np.conjugate(z0)) + zt = complex(z.real*np.cos(z.imag), z.real*np.sin(z.imag)) + z0t = complex(z0.real*np.cos(z0.imag), z0.real*np.sin(z0.imag)) + res = (zt-z0t)/(1-zt*np.conjugate(z0t)) + + mz0 = z0.real + az0 = z0.imag + mz = z.real + az = z.imag -def moeb_origin_trafo_inverse(z0, z): - return (z+z0) / (1+z*np.conjugate(z0)) + retval = 0 + anglediff = z.imag - z0.imag + if math.isclose(az, 0.0) and math.isclose(az0, np.pi): + retval = (mz+mz0)/(1+mz*mz0) + if math.isclose(az, np.pi) and math.isclose(az0, np.pi): + retval = (-mz+mz0)/(1-mz*mz0) + if math.isclose(az, 0.0) and math.isclose(az0, 0.0): + retval = (mz-mz0)/(1-mz*mz0) + if math.isclose(az, np.pi) and math.isclose(az0, 0.0): + retval = (-mz-mz0)/(1+mz*mz0) + if (retval != 0): + if (retval < 0): + retval = complex(np.abs(retval), np.pi) + return retval + # now comes treatment of congruent angle differences... + if math.isclose(anglediff, 0.0): + retval = (mz-mz0)/(1-mz*mz0) + if math.isclose(anglediff, np.pi): + retval = (-mz-mz0)/(1 + mz*mz0) + if (retval != 0): + if (retval < 0): + retval = complex(np.abs(retval), np.pi) + return complex(retval.real, retval.imag + az0) + + if mz == 0: + return myclamp(complex(z0.real, z0.imag + np.pi)) + if z0 == z: + return 0 + y = myclamp(complex(z.real, z.imag - az0)) + r1 = mz/mz0 + r2 = mz*mz0 + ca = np.cos(az - az0)# shift angle for simpler form of transform. + sa = np.sin(az - az0) + nn1 = r1 + 1/r1 - 2*ca + nn2 = r2 + 1/r2 - 2*ca + mag = np.sqrt(nn1/nn2) + #phi = np.arctan((np.sin(az - az0)*(1-mz0*mz0))/(ca*(1+mz0*mz0) - 1/r1 * (1-mz*mz) )) + #p1 = np.arctan(mz*sa/(mz*ca - mz0)) + #p2 = np.arctan(-mz0*mz*sa/(1-mz0*mz*ca)) + #print (np.angle(y - mz0), y, z, z0) +## p3 = 2*np.arctan(mz * sa/(mz*ca - mz0 + np.sqrt( mz*mz0 * (r1+1/r1 - 2*ca)))) + ## p4 = 2*np.arctan(-mz*mz0*sa/(np.sqrt(mz*mz0*(r2 + 1/r2 - 2*ca)) + 1 - mz*mz0*ca ) ) + n1 = np.sqrt(nn1*mz*mz0) + n2 = np.sqrt(nn2*mz*mz0) + phi = 2*np.arctan(mz*sa * (n2 + 1 + mz0*n1 - mz0*mz0) / ((n1 - mz0)*(n2 + 1 - mz*mz0*ca) + mz*ca*(n2 + 1) - mz*mz*mz0 )) + #print (z0, z, np.abs(res), np.angle(res), mag, phi + az0) +# return complex(np.abs(res), np.angle(res)) + return myclamp(complex(mag, phi + az0)) - # rotates z by phi counter-clockwise about the origin -def moeb_rotate_trafo(z, phi): - return z * np.exp(complex(0, phi)) +def moeb_origin_trafo_inverse(z0, z): + y = myclamp(complex(z0.real, z0.imag + np.pi)) + return moeb_origin_trafo(y, z) + # rotates z by phi counter-clockwise about the origin +def moeb_rotate_trafo(z, phi): +# return z * np.exp(complex(0, phi)) + return myclamp(complex(z.real, z.imag + phi)) def moeb_translate_trafo(z, s): num = z-s