Hi list,

I've been working on using PyX to plot Minkowski spacetime diagrams for work. 
My first attempt used the basic plotting functions to do the job, but I was a 
unsatisfied with not making use of PyX's powerful graphing capabilities. So 
I've made a second attempt, this time implementing classes that extend those 
defined by PyX, so that it integrates nicer.

I'm posting this because (a) someone else might also find it useful, and (b) 
I'd appreciate feedback.

It uses scipy (for speed of light) and numpy (for dot product and array), 
neither of which is probably necessary (but I'm lazy). The way it works is 
that you define events (represented by the event class) that occur at a 
location in 4D spacetime in some inertial frame (represented by the frame 
class). (It might seem like overkill to consider all 4 dimensions, but it was 
necessary in my case.)

The events are collected in a chronology which acts as a data provider (by 
projecting the 4D events down to a time and position along a line in space). 
The minkowski class inherits graphxy and keeps track of the inertial frame in 
which the events are drawn, and the line in space onto which the events are 
projected.

I've also created a style that draws a forward and backward lightcone (which 
can even be filled; transparent fills work nicely). I've included a test file 
to 
see how it basically works.

One thing that I am a bit stumped on is how I would have this system draw the 
axes of an alternative frame that is moving at some velocity relative to the 
one in which the events are drawn (the "local" frame). The axes for such a 
frame "squeeze" together, i.e. the x and t axes rotate in opposite directions. 
I'm not sure how to approach that. One idea I had was to extend the function 
of minkowski to draw these axes manually, but even that has me unsure of how 
to proceed. Any suggestions?

Peace,
Brendon
# -*- coding: utf-8 -*-

from pyx import *
from math import *
from scipy import constants
from numpy import array, dot

# Trafos
def lorentz_alpha(v):
	return sqrt(1 - (v/constants.c)**2)

def lorentz_gamma(v):
	return 1.0/lorentz_alpha(v)

#def lorentz_space(x, t, v):
	#return lorentz_gamma(v)*(x - v*t)

#def lorentz_time(x, t, v):
	#return lorentz_gamma(v)*(t - v*x/constants.c**2)

def lorentz_trafo(v):
	s = sqrt(dot(v, v))
	nv = v/s
	g = lorentz_gamma(s)
	mgv = -g*v
	gm1nv = (g - 1.)*nv
	return array([[g, mgv[0], mgv[1], mgv[2]],
		[mgv[0], gm1nv[0]*nv[0] + 1., gm1nv[1]*nv[0], gm1nv[2]*nv[0]],
		[mgv[1], gm1nv[0]*nv[1], gm1nv[1]*nv[1] + 1., gm1nv[2]*nv[1]],
		[mgv[2], gm1nv[0]*nv[2], gm1nv[1]*nv[2], gm1nv[2]*nv[2] + 1.]])

class frame:
	"""An inertial reference frame in which events may take place. Frames have
	common origins but may be moving relative to each other."""
	
	def __init__(self):
		# These have to be initialized in the constructor otherwise every
		# instance points to the same dictionary object.
		self.relatives_v = dict()  #: Relative velocities of other frames
		self.relatives_t = dict()  #: Lorentz transformation matrices to other frames
	
	def equivalent(self, other):
		return self is other or self.relatives_v == other.relatives_v
	
	def relate(self, other, v):
		"""Relative to this frame, the other frame is moving at v."""
		self.relatives_v[other] = v
		self.relatives_t[other] = lorentz_trafo(v)
		other.relatives_v[self] = -v
		other.relatives_t[self] = lorentz_trafo(-v)
	
	def related(self, other):
		return other in self.relatives

class event:
	"""Some event that occurs in spacetime, defined in terms of some reference
	frame."""
	
	def __init__(self, f, x, label = None):
		"""An event is a point in spacetime.
		
		f: the reference frame in which the event is defined
		x: location of the event, i.e. array([t, x, y, z])
		label: The label to draw on the Minkowski diagram, or None
		"""
		
		self.f = f
		self.x = x
		self.label = label
	
	def is_local(self, other):
		"""Does the event occur within the lightcone of another, such that they could communicate locally?"""
		d = self.x - other.x
		# Even though position changes in different frames, the lightcone test is valid in all
		return (constants.c * d[0])**2 >= (d[1]**2 + d[2]**2 + d[3]**2)
	
	def in_frame(self, f):
		"""Return an event equivalent to this defined in the given reference
		frame f. Will return self if the frame is the one used for definition.
		Will fail if the frames are not related.
		"""
		if f.equivalent(self.f):
			return self
		return event(f, dot(self.f.relatives_t[f], self.x))
	
	def lightcone_timelead(self, other):
		"""Returns how many seconds this event occurs before the lightcone of
		an other event reaches this location. Time is given in self's frame.
		May be negative if this event is within the lightcone of the other
		event."""
		d = self.x - other.in_frame(self.f).x
		dist = sqrt(d[1]**2 + d[2]**2 + d[3]**2)
		return dist/constants.c - d[0]

# An issue:
# 1) How to draw the slanted axes of moving frames?
# Old code to illustrate:
#	def axes(self, frame = None, xlabel = '$x$', tlabel = '$t$', attrs = [deco.earrow.normal]):
#		alpha1 = 0.0
#		alpha2 = pi/2.0
#		if frame is not None and not frame.equivalent(self.frame):
#			# Determine the velocity in the projected direction
#			v = dot(self.frame.relatives_v[frame], self.space)
#			alpha1 = atan(v/constants.c*self.tfactor/constants.c/self.xfactor)
#			alpha2 = atan(constants.c/v*self.tfactor/constants.c/self.xfactor)
#		hsize = self.size/2.0
#		hsca1 = hsize*cos(alpha1)
#		hssa1 = hsize*sin(alpha1)
#		hsca2 = hsize*cos(alpha2)
#		hssa2 = hsize*sin(alpha2)
#		self.canvas.insert(vector(-hsca1, -hssa1, hsca1, hssa1, xlabel, pos = 1, attrs = attrs))
#		self.canvas.insert(vector(-hsca2, -hssa2, hsca2, hssa2, tlabel, pos = 1, attrs = attrs))

class chronology(graph.data._data):
	"""A data provider for Minkowski diagrams: A story of events in spacetime,
	automatically adjusted to the inertial frame of the diagram."""
	
	defaultstyles = graph.data.defaultsymbols
	
	def __init__(self, events, title=None):
		self.events = events
		self.title = title
		self.xname = "x"
		self.yname = "y"
		self.columns = {}
		self.columnnames = [self.xname, self.yname]
	
	def dynamiccolumns(self, graph):
		# NOTE graph must be a minkowski, so we can know what frame and
		# projection line (spaceline) to use
		frame = graph.frame
		spaceline = graph.spaceline
		dynamiccolumns = {self.xname: [], self.yname: []}
		for e in self.events:
			x = e.in_frame(frame).x
			dynamiccolumns[self.xname].append(dot(x.take([1, 2, 3]), spaceline))
			dynamiccolumns[self.yname].append(x[0])
		return dynamiccolumns

class lightcone(graph.style._styleneedingpointpos):
	"""A style to draw lightcones of events on a Minkowski diagram."""
	
	needsdata = ["vpos", "vposavailable"]
	
	def __init__(self, attrs=[deco.stroked([style.linestyle.dashed])]):
		self.attrs = attrs
	
	def initdrawpoints(self, privatedata, sharedata, graph):
		x0, y0 = graph.vpos_pt(0, 0)
		x1, y1 = graph.vpos_pt(1, 1)
		privatedata.lightconecanvas = canvas.canvas([canvas.clip(path.rect_pt(
			x0, y0, x1 - x0, y1 - y0))])
	
	def drawpoint(self, privatedata, sharedata, graph, point):
		if sharedata.vposavailable:
			x_pt, y_pt = graph.vpos_pt(*sharedata.vpos)
			# Figure out the graph extents
			gx0, gy0 = graph.vpos_pt(0, 0)
			gx1, gy1 = graph.vpos_pt(1, 1)
			# Guess the lightcone directions
			llx, lly = graph.pos_pt(point["x"] + constants.c, point["y"] - 1)
			lux, luy = graph.pos_pt(point["x"] + constants.c, point["y"] + 1)
			# Project them to the top/bottom of the graph
			llhwidth = (llx - x_pt)*(gy0 - y_pt)/(lly - y_pt)
			luhwidth = (lux - x_pt)*(gy1 - y_pt)/(luy - y_pt)
			# Draw them as open triangles
			privatedata.lightconecanvas.draw(path.path(
				path.moveto_pt(x_pt - llhwidth, gy0),
				path.lineto_pt(x_pt, y_pt),
				path.lineto_pt(x_pt + llhwidth, gy0)), self.attrs)
			privatedata.lightconecanvas.draw(path.path(
				path.moveto_pt(x_pt - luhwidth, gy1),
				path.lineto_pt(x_pt, y_pt),
				path.lineto_pt(x_pt + luhwidth, gy1)), self.attrs)
	
	def donedrawpoints(self, privatedata, sharedata, graph):
		graph.insert(privatedata.lightconecanvas)

class minkowski(graph.graphxy):
	"""Construct a Minkowski spacetime diagram."""
	
	def __init__(self, frame, spaceline, xpos=0, ypos=0, width=None, height=None,
			ratio=graph.graph.goldenmean, backgroundattrs=None,
			x=graph.axis.linear(title="$x$"), t=graph.axis.linear(title="$t$")):
		"""We use standard 0-centred axes. TODO No key, for now, but maybe it
		would be useful."""
		graph.graphxy.__init__(self, xpos=xpos, ypos=ypos, width=width, height=height,
			ratio=ratio, key=None, backgroundattrs=None, xaxisat=0, yaxisat=0,
			x=x, y=t)
		self.frame = frame
		self.spaceline = spaceline
# -*- coding: utf-8 -*-

from minkowski import *
from pyx import *
from math import *
from scipy import constants
from numpy import array

f1 = frame()

f2 = frame()

f1.relate(f2, array([0.5*constants.c, 0.0, 0.0]))

e1 = event(f1, array([0.0, 0.0, 0.0, 0.0]))
e2 = event(f1, array([0.0001, 60000, 0.0, 0.0]))

m = minkowski(f1, array([1.0, 0.0, 0.0]), width=8.5,
	x=graph.axis.linear(title="$x$", min=-120000, max=120000),
	t=graph.axis.linear(title="$t$", min=-0.0002, max=0.0002))

d = chronology([e1, e2])

m.plot(d, styles=[graph.style.symbol(graph.style.symbol.circle), lightcone()])
#m.plot(d, styles=[graph.style.symbol(graph.style.symbol.circle)])

m.writePDFfile("minkowski_test")

Attachment: signature.asc
Description: This is a digitally signed message part.

------------------------------------------------------------------------------
Benefiting from Server Virtualization: Beyond Initial Workload 
Consolidation -- Increasing the use of server virtualization is a top
priority.Virtualization can reduce costs, simplify management, and improve 
application availability and disaster protection. Learn more about boosting 
the value of server virtualization. http://p.sf.net/sfu/vmware-sfdev2dev
_______________________________________________
PyX-user mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/pyx-user

Reply via email to