r""" Python script for drawing trajectories in the periodic wind-tree model. The script is intended to draw orbit under the billiard flow in the windtree table (regular scatterers of size a times b are displayed periodically in the plane). Here is a short sample: >>> from windtree_pictures import windtree_trajectory >>> from math import pi >>> a = b = 0.5 >>> theta = pi/2.112 >>> p = ((0.1,0.1),0) # initial point as ((x,y),square) >>> l = windtree_trajectory(a,b,theta,p,N=200) >>> l[0] (0.10000000000000001, 0.10000000000000001, (1, 1)) >>> l[1] (0.5, 4.8908202153682216, (1, -1)) >>> l[2] (0, 10.879345484578497, (1, 1)) >>> l[3] (0.5, 16.867870753788775, (1, -1)) The function returns a list of 3-tuple (x,y,j) where x,y is a point on a scatterers and j is the relative orientation of the particle. To get a tikz picture (to be compiled with latex or pdflatex), there is a dedicated command:: >>> from windtree_pictures import windtree_trajectory >>> from math import pi >>> a = b = 0.5 >>> theta = pi/2.112 >>> p = ((0.1,0.1),0) # initial point as ((x,y),square) >>> windtree_tikz_file(a,b,theta,p,N=200,filename="mypicture.tex") """ ################################################################################ # # Copyright (C) 2009, Vincent Delecroix # # Distributed under the terms of the GNU General Public License (GPL) # # The full text of the GPL is available at: # # http://www.gnu.org/licenses/ # ################################################################################ from math import cos,sin,tan class HypIET: r""" Hyperelliptic interval exchange transformation for the rectangular L-shaped surface. The parameters a, b and theta are respectively the width and the height of the rectangular obstacle and the angle of the initial state of the particle. The numbering of the wedges is as follow: .---. | 2 | .---.---. | 0 | 1 | .---.---. """ def __init__(self,a,b,theta): self.a = a self.b = b self.theta = theta self.lambda_l = [(1-b)*cos(theta),(1-b)*cos(theta),b*cos(theta)] self.lambda_r = [(1-a)*sin(theta),a*sin(theta),(1-a)*sin(theta)] self.pi_l = [2,1,0] self.pi_r = [1,0,2] # check the train track condition for i in xrange(3): assert (abs( self.lambda_l[i] + self.lambda_r[i] - self.lambda_r[self.pi_l[i]] - self.lambda_l[self.pi_r[i]] ) < 0.0001) def _repr_(self): return "iet for a=%f b=%f and theta=%f"%(self.a,self.b,self.theta) def __call__(self,p): r""" p should be a 2-tuple """ x,i = p if x < -self.lambda_l[i] or x > self.lambda_r[i]: raise ValueError, "%s not in the interval %d"%(x,i) xx = x+self.lambda_l[i] ii = self.pi_l[i] if xx < self.lambda_r[ii]: return xx,ii xx = self.lambda_r[i]-x ii = self.pi_r[i] if xx < self.lambda_l[ii]: return -xx,ii raise ValueError, "wrong values for p=(%s,%s)"%(x,i) def L_trajectory(a,b,theta,p,N=20,verbose=False): r""" INPUT: - ``a,b,theta`` - parameters of the wind-tree table - ``p`` - 2-tuple - initial point - ``N`` - positive integer - number of iterations """ (x,y),i = p res = [(x,y)] t = sin(theta)*x - cos(theta)*y T = HypIET(a,b,theta) shifts = [(0,0),(1-a,0),(0,1-b)] widths = [1-a,a,1-a] heights = [1-b,1-b,b] if verbose: print "data of the iet:" print " lambda_1 = (%f,%f)"%(T.lambda_l[0],T.lambda_r[0]) print " lambda_2 = (%f,%f)"%(T.lambda_l[1],T.lambda_r[1]) print " lambda_3 = (%f,%f)"%(T.lambda_l[2],T.lambda_r[2]) print "iterates:" for _ in xrange(N): t,ii = T((t,i)) if verbose: print " (%d, %f)"%(i+1,t) if t < 0: res.append((shifts[i][0]+widths[i],shifts[i][1]-t/cos(theta))) res.append((shifts[ii][0],shifts[ii][1]-t/cos(theta))) else: res.append((shifts[i][0]+t/sin(theta),shifts[i][1]+heights[i])) res.append((shifts[ii][0]+t/sin(theta),shifts[ii][1])) i = ii return res def L_tikz_file(a,b,theta,p,N=20,filename=None): if filename is None: filename = "picture.tex" f = open(filename, "w") f.write("\\begin{tikzpicture}\n") f.write("\\fill[gray] (0,0) -- (1,0) -- (1,%f) -- (%f,%f) -- (%f,1) -- (0,1) -- cycle;\n" %(1-b,1-a,1-b,1-a)) l = L_trajectory(a,b,theta,p,N) for i in xrange(len(l)//2): x0,y0 = l[2*i] x1,y1 = l[2*i+1] f.write("\\draw (%f,%f) -- (%f,%f);\n"%(x0,y0,x1,y1)) f.write("\\end{tikzpicture}\n") f.close() def L4_trajectory(a,b,theta,p,N=20,verbose=False): r""" INPUT: - ``a,b,theta`` - parameters of the wind-tree table - ``p`` - 2-tuple - initial point - ``N`` - positive integer - number of iterations """ (x,y),i = p j = (1,1) res = [(x,y)] t = sin(theta)*x - cos(theta)*y T = HypIET(a,b,theta) d = 1.5 shifts = { ( 1, 1): [(0,0),(1-a,0),(0,1-b)], ( 1,-1): [(d+a,0),(d,0),(d+a,1-a)], (-1, 1): [(0,d+b),(1-a,d+b),(0,d)], (-1,-1): [(d+a,d+b),(d,d+b),(d+a,d)] } widths = [1-a,a,1-a] heights = [1-b,1-b,b] if verbose: print "data of the iet:" print " lambda_1 = (%f,%f)"%(T.lambda_l[0],T.lambda_r[0]) print " lambda_2 = (%f,%f)"%(T.lambda_l[1],T.lambda_r[1]) print " lambda_3 = (%f,%f)"%(T.lambda_l[2],T.lambda_r[2]) print "iterates:" for _ in xrange(N): t,ii = T((t,i)) jj = j if ii == 1 and t > 0: jj = (-j[0],j[1]) elif ii == 2 and t < 0: jj = (j[0],-j[1]) if verbose: print " (%d, %f), (%d,%d)"%(ii+1,t,j[0],j[1]) if t < 0: res.append((shifts[j][i][0]+widths[i],shifts[j][i][1]-t/cos(theta))) res.append((shifts[jj][ii][0],shifts[jj][ii][1]-t/cos(theta))) else: res.append((shifts[j][i][0]+t/sin(theta),shifts[j][i][1]+heights[i])) res.append((shifts[jj][ii][0]+t/sin(theta),shifts[jj][ii][1])) i = ii j = jj return res def windtree_trajectory(a,b,theta,p,N=20,intermediate=False,verbose=False): r""" INPUT: - ``a,b,theta`` - parameters of the wind-tree table - ``p`` - 2-tuple - initial point - ``N`` - positive integer - number of iterations - ``intermediate`` - if True, print the return map on the natural rectangle subdivision of the billiard table, if False return only the time for which a scatterer is reached. OUTPUT: A list of 3-tuples (x,y,or) where or is either ``(1,1)``, or ``(1,-1)`` or ``(-1,1)`` or ``(-1,-1)``. """ (x,y),i = p # position j = (1,1) # orientation v = (0,0) # copy dx=dy=0 res = [(x,y,j)] t = sin(theta)*x - cos(theta)*y T = HypIET(a,b,theta) d = 1.5 shifts = [(0,0),(1-a,0),(0,1-b)] widths = [1-a,a,1-a] heights = [1-b,1-b,b] if verbose: print "data of the iet:" print " lambda_1 = (%f,%f)"%(T.lambda_l[0],T.lambda_r[0]) print " lambda_2 = (%f,%f)"%(T.lambda_l[1],T.lambda_r[1]) print " lambda_3 = (%f,%f)"%(T.lambda_l[2],T.lambda_r[2]) print "iterates:" print " position dir. level" for _ in xrange(N): change = False # compute new point in L (t and i) t,i = T((t,i)) # update cocycle (j and v) if i == 1 and t > 0: j = (-j[0],j[1]) if j[0] == 1: dy = 0 else: dy = 1-b change = True elif i == 2 and t < 0: j = (j[0],-j[1]) if j[1] == 1: dx = 0 else: dx = 1-a change = True if i == 0: if t < 0: v = (v[0]+j[1],v[1]) else: v = (v[0],v[1]+j[0]) if verbose: print " (%d, %f), (%d,%d), (%d,%d)"%(i+1,t,j[0],j[1],v[0],v[1]) if change or intermediate: if t < 0: res.append(( v[0]+j[1]*shifts[i][0]+dx, v[1]+j[0]*shifts[i][1]+dy-j[0]*t/cos(theta), j)) else: res.append(( v[0]+j[1]*shifts[i][0]+dx+j[1]*t/sin(theta), v[1]+j[0]*shifts[i][1]+dy, j)) return res def windtree_tikz_file(a,b,theta,p,N=20,intermediate=False,filename=None): r""" Write a latex file with name ``filename`` (default is picture.tex) which contains a tikz picture. You should include this file in the main document using an \input statement. INPUT: - ``a,b,theta,p,N`` - parameters sent to the function ``windtree_trajectory`` - ``filename`` - name of the file """ l = windtree_trajectory(a,b,theta,p,N,intermediate) if filename is None: filename = "picture.tex" f = open(filename,"w") colors = { (+1,+1): 'pp', (+1,-1): 'pm', (-1,+1): 'mp', (-1,-1): 'mm'} xmin=xmax=ymin=ymax=0 for x,y,_ in l: if x < xmin: xmin=floor(x) if x > xmax: xmax=ceil(x) if y < ymin: ymin=floor(y) if y > ymax: ymax=ceil(y) xmin -= 1 xmax -= 1 ymin -= 1 ymax -= 1 f.write("\\begin{tikzpicture}\n") f.write("\\useasboundingbox (%f,%f) rectangle (%f,%f);\n"%(xmin,ymin,xmax+1+a,ymax+1+b)) f.write("\\definecolor{pp}{rgb}{1.0,0.0,0.0};\n") f.write("\\definecolor{pm}{rgb}{0.8,0.8,0.0};\n") f.write("\\definecolor{mp}{rgb}{0.1,0.8,0.1};\n") f.write("\\definecolor{mm}{rgb}{0.5,0.0,1.0};\n") f.write("\\foreach \\x in {%d,...,%d} \\foreach \\y in {%d,...,%d}\n"%(xmin,xmax,ymin,ymax)) f.write(" \\fill[fill=blue,opacity=0.2] (\\x,\\y) +(%f,%f) rectangle +(1,1);\n"%(1-a,1-b)) for i in xrange(len(l)-1): x0,y0,j = l[i] x1,y1,_ = l[i+1] f.write("\\draw[color=%s] (%f,%f) -- (%f,%f);\n"%(colors[j],x0,y0,x1,y1)) f.write("\n") f.write("\\end{tikzpicture}\n") f.close()