# super.py # # by D. Homan, homan + d + "at" + denison.edu # # This script is provided free with ABSOLUTELY NO WARRANTY # of any kind. # # The script has only been tested on my local macintosh computer # system and may still contain bugs and/or errors. # # This is a python script to simulate superluminal motion using # the vpython library. http://www.vpython.org/ You must run this # program under the "VIDLE" interface to get it to work... just start-up # VIDLE, load this file, then select 'Run Module' from the Run # menu. # # NOTE: The script assumes a screen size of at least 1300x820 pixels # -- slightly smaller displays will probably work OK, but for # significantly smaller displays, one will have to edit the # window sizes and offsets given below vpython_window_width = 800 # pixels vpython_window_height = 800 # pixels xy_plot_offset = 500 # horizontal pixel offset xy_plot_width = 800 # pixels xy_plot_height = 600 # pixels # # The parts of this script that allow one to make a movie come # from the example program movie.py by Kelvin Chu on the vpython site... # http://vpython.org/contributed/movie.py # # Right now the movie creation only works on a mac, but it should be # possible to do something similar on a linux box. The frames # of the movie are placed in a 'frames' sub-directory. They are # cropped by the 'convert' command from the imagemagick package. # And I use quicktime to put the sequence together into a movie. # # Change Log # ----------- # June 2010 -- simple version created, just 3-D graphics # Oct 6-7, 2010 -- Added comments, 2-D graphs of motion, movie code # Oct 8, 2010 -- Fixed emitted photon distance to actually be the distance # along the jet rather than the distance from the origin. # -- Improved scaling of X,Y plot for a variety of different # initial conditions # -- Improved quantity and quality of comments # -- fixed a number of bugs # -- Added simple bends and accelerations to the parameters # at the top of the script # May 26, 2011 -- small modifications for release on website from visual import * from visual.graph import * from commands import * ############################################################ # Define some basic parameters # Obs_plane_distance = 80 # distance to place the plane of observation Length_of_Jet = 50 # total distance 'blob' travels down jet # before stopping Beta = 0.95 # speed of jet 'blob' in units of speed of light Theta = 18 # angle to the line of sight in degrees Phi = 0.0 # azmuthal angle on the sky in degrees Number_of_Photons = 10 # Number of photons = number of observations make_movie = False # Make a movie of the motion Number_of_Frames = 300 # Number of Movie frames ############################################################ # Define some more advanced parameters # # parameters for a simple bend in jet # Delta_theta = 0 # jet bend in degrees Delta_phi = 0 # jet bend in degrees Delta_theta = Delta_theta*3.14159/180.0 # conversion to radians Delta_phi = Delta_phi*3.14159/180.0 # # parameters for simple linear acceleration of moving component # Gamma_dot = 0.0 # time derivative of the Lorentz factor, # should be pretty small, especially if # negative to avoid getting # to some sort of negative speed ############################################################ ################################################ # # Function to grab a screen shot and trim it # to the appropriate size... only used in # movie making mode # ################################################# def GetScreenShot(FrameNumber): # Take a screenshot and write # it to a numbered file. tmp = getoutput('/usr/sbin/screencapture frames/temp.png') tmp = getoutput('/sw/bin/convert -crop 1300x820+0+0 frames/temp.png frames/super.%03d.png' % FrameNumber) tmp = getoutput('/sw/bin/convert -resize 1040x656 frames/super.%03d.png frames/super_small.%03d.png' % (FrameNumber, FrameNumber)) tmp = getoutput('/sw/bin/convert -resize 650x410 frames/super.%03d.png frames/super_tiny.%03d.png' % (FrameNumber , FrameNumber)) print 'Frame: %d' % (FrameNumber) return ################################################ # # Function to setup some basic objects in # the 3-D space # ################################################ def setup_scene(): xaxis = arrow(pos=(0,0,0), axis=(20,0,0), shaftwidth=1) yaxis = arrow(pos=(0,0,0), axis=(0,20,0), shaftwidth=1) zaxis = arrow(pos=(0,0,0), axis=(0,0,20), shaftwidth=1) xaxislabel = label(pos=(10,0,0), text='x axis', xoffset=0, yoffset=+20) yaxislabel = label(pos=(0,10,0), text='y axis', xoffset=-20, yoffset=0) zaxislabel = label(pos=(0,0,10), text='toward observer', xoffset=-30, yoffset=0, zoffset=-20) obsplane = box(pos=(0,0,Obs_plane_distance),axis=(1,0,0), length=Obs_plane_distance,height=Obs_plane_distance, width=0.2, opacity=0.2, color=color.blue) obsxaxis = arrow(pos=(0,0,Obs_plane_distance), axis=(10,0,0), shaftwidth=0.5, opacity=0.3) obsyaxis = arrow(pos=(0,0,Obs_plane_distance), axis=(0,10,0), shaftwidth=0.5, opacity=0.3) obsplanelabel = label(pos=(Obs_plane_distance/2,-Obs_plane_distance/4,Obs_plane_distance), text="Plane of Observation") return ################################################ # # Main function that executes the program # ################################################ def Main(): ################################################ # # Define the "scene" display for vpython # ################################################ scene = display(title='SuperLuminal Motion', x=0, y=0, width=vpython_window_width, height=vpython_window_height, center=(15,0,25), forward=(0,-1,-1), scale=(0.03,0.03,0.03)) setup_scene() ################################################ # # Define array for photon objects # ################################################ photons = [] ################################################ # # Define jet component or 'blob' that will # travel down the jet # ################################################ blob = sphere(color=color.orange,opacity=0.3) blob.beta = Beta # speed in units of the speed of light blob.theta = Theta*3.141519/180.0 # angle to line of sight in radians blob.phi = Phi*3.141519/180.0 # azimuthial velocity angle in radians blob.distance = 0 # total distance traveled down the jet blob.path = curve(color=color.orange) # curve to trace path of the blob in 3-D # vectors for the position, velocity, and acceleration on the component # -- these vectors have (x,y,z) components. # x is to the right, y is up, and z is out of the screen blob.pos = vector(0,0,0) blob.vel = vector(blob.beta*sin(blob.theta)*cos(blob.phi), blob.beta*sin(blob.theta)*sin(blob.phi), blob.beta*cos(blob.theta)) blob.Gamma = 1.0/sqrt(1-blob.beta*blob.beta) blob.beta_dot = Gamma_dot/(blob.beta*blob.Gamma*blob.Gamma*blob.Gamma) blob.acc = vector(blob.beta_dot*sin(blob.theta)*cos(blob.phi), blob.beta_dot*sin(blob.theta)*sin(blob.phi), blob.beta_dot*cos(blob.theta)) ################################################ # # Define the direction and geometical object # that will serve as the jet in our simulation # This is done in two pieces to allow the possibility # of a simple bend in the jet. More complicated # jets are possible, but require reprogramming # this script # ################################################ jet_vector=vector((Length_of_Jet/2)*sin(blob.theta)*cos(blob.phi), (Length_of_Jet/2)*sin(blob.theta)*sin(blob.phi), (Length_of_Jet/2)*cos(blob.theta)) jet_vector2=vector((Length_of_Jet/2)*sin(blob.theta+Delta_theta)*cos(blob.phi+Delta_phi), (Length_of_Jet/2)*sin(blob.theta+Delta_theta)*sin(blob.phi+Delta_phi), (Length_of_Jet/2)*cos(blob.theta+Delta_theta)) if Delta_theta == 0 and Delta_phi == 0: jet_vector *= 2.0 jet = cylinder(color=color.green,opacity=0.3, pos=(0,0,0), radius=0.5, axis=jet_vector) else: jet = cylinder(color=color.green,opacity=0.3, pos=(0,0,0), radius=0.5, axis=jet_vector) jet2 = cylinder(color=color.green,opacity=0.3, pos=jet_vector, radius=0.5, axis=jet_vector2) ################################################ # # Setup an X,Y graph for displaying the position # versus time as 'photons' are emitted and # later observed # ################################################ # set x,y plot limits y_plot_max = Length_of_Jet*1.3 x_plot_max = (Length_of_Jet*(1-cos(blob.theta))+Obs_plane_distance)*1.1 extra_offset_factor = 1.0 if x_plot_max > 2.0*y_plot_max: extra_offset_factor = 1.5 # # setup the x,y plot for distance versus time # emit_graph = gdisplay(x=500, y=0, width=800, height=600, title='Distance vs. Time', xtitle='Time', ytitle='Distance (light-years)', xmax=x_plot_max, xmin=0., ymax=y_plot_max, ymin=0) emit_graph.display.uniform=True # foreground=color.black, background=color.white) emit_plotting = gdots(display=emit_graph.display, color=color.blue, shape='square') obs_plotting = gdots(display=emit_graph.display, color=color.green, shape='square') label(display=emit_graph.display, pos=(x_plot_max*0.3,y_plot_max*0.9*extra_offset_factor), color=color.blue, height=20, box=0, text='Blue = Distance Along Jet') label(display=emit_graph.display, pos=(x_plot_max*0.8,y_plot_max*0.9*extra_offset_factor), color=color.green, height=20, box=0, text='Green = Distance Across Sky') ############################################################ # # What follows below is a loop with very small time steps, dt # -- At each step in the loop, the acceleration (if any) makes # a small change in the velocity of the jet component, # and the velocity makes a small change in the position. # The component motion stops when it reaches the defined # end of the jet. # -- photons are emitted at regular intervals and travel # toward the observer at a constant speed of 1.0 (speed of light) # The positions of the photons are also updated at each # step in the loop # -- photons accumulate on plane which defines the plane of # observation and illustrates how the jet component # appears to move on the plane of the sky when observed # from far away # ############################################################ # Setup parameters for Loop for Jet motion and Photon Travel ############################################################ dt = 0.01 # a small time step total_time = 0.0 # time since start of simulation loop_count = 0 # count of number of times loop executes zoomed_in = True # is our viewpoint zoomed in on the jet? # -- used to tell us it is OK to zoom out # once the component reaches the end # the jet update_path = True # Flag to tell us whether we still need # update the path after the 'bend' point # point in the jet # interval for emitting photons photon_step_interval = int(Length_of_Jet/(Number_of_Photons*dt)) # interval for making movie frames movie_loop_interval = int(Obs_plane_distance/int((Number_of_Frames-2)*dt)) FrameNumber = 1 # initialize frame number # Get a movie frame before the motion starts and update frame number if make_movie: GetScreenShot(FrameNumber) FrameNumber += 1 ################################################ # Start of Loop for Jet motion and Photon Travel ################################################# while 1: rate(500) # rate at which action unfolds #-------------------------------------------- # If jet component is still in jet # -- update its velocity and position # including any acceleration or bending # -- append a line to its path so far # Else # -- check if we need to zoom out our point of view #--------------------------------------------- if blob.distance < Length_of_Jet: # update the path if we have gone past the point # where there is a bend in the jet... note the # 'bend' might be zero degrees, in which case # this will have no effect if update_path and blob.distance > Length_of_Jet/2.0: update_path = False blob.theta += Delta_theta blob.phi += Delta_phi blob.vel = vector(blob.beta*sin(blob.theta)*cos(blob.phi), blob.beta*sin(blob.theta)*sin(blob.phi), blob.beta*cos(blob.theta)) # update the acceleration of the blob using change in Gamma blob.beta_dot = Gamma_dot/(blob.beta*blob.Gamma*blob.Gamma*blob.Gamma) blob.acc = vector(blob.beta_dot*sin(blob.theta)*cos(blob.phi), blob.beta_dot*sin(blob.theta)*sin(blob.phi), blob.beta_dot*cos(blob.theta)) # Calculate the new velocity and speed of the blob blob.vel = blob.vel + blob.acc*dt blob.beta = blob.vel.mag blob.Gamma = 1.0/sqrt(1-blob.beta*blob.beta) # Calculate the position of the blob blob.pos = blob.pos + blob.vel*dt # Calculate the total distance traveled along jet blob.distance = blob.distance + blob.beta*dt # add to the path tracking the blob blob.path.append(pos=blob.pos) else: if zoomed_in: zoomed_in = False scene.scale *= 0.7 scene.center = (30,0,50) scene.forward = (0,-1.5,-1) #-------------------------------------------- # Loop over each photon we currently have # -- check if it has reached the plane # of observation: if so, color it # green and add a point to the 2-D # of distance versus time. if not, # then update its position and time #-------------------------------------------- no_more_update = True for i in photons: if i.pos.z < Obs_plane_distance: i.pos = i.pos + i.vel*dt i.obs_time = total_time no_more_update = False else: i.color = color.green obs_plotting.plot(pos=(i.obs_time,sqrt(i.pos.x*i.pos.x+i.y*i.pos.y))) #-------------------------------------------- # if all photons have reached the obs. plane # we are done! #-------------------------------------------- if (len(photons) > 0) and no_more_update: break #-------------------------------------------- # tracking the time and loop count #-------------------------------------------- total_time = total_time + dt loop_count += 1 #-------------------------------------------- # Check if we need to emit a new photon # while the jet component is still moving #-------------------------------------------- if not(loop_count%photon_step_interval) and blob.distance < Length_of_Jet: photons.append(sphere(color=color.blue, radius=0.5)) photons[len(photons)-1].pos = blob.pos photons[len(photons)-1].vel = vector(0.,0.,1.0) photons[len(photons)-1].emit_time = total_time photons[len(photons)-1].emit_dist = blob.distance emit_plotting.plot(pos=(total_time,blob.distance)) #print photons #-------------------------------------------- # Grab a movie frame if we are making a movie #-------------------------------------------- if make_movie and not(loop_count%movie_loop_interval): GetScreenShot(FrameNumber) FrameNumber += 1 ################################################ # End of Loop for Jet motion and Photon Travel ################################################# ################################################ # # Find min and max emission and observations # times and locations and computer the actual # and apparent speeds # ################################################ min_emit_time = photons[0].emit_time max_emit_time = photons[len(photons)-1].emit_time min_emit_dist = photons[0].emit_dist max_emit_dist = photons[len(photons)-1].emit_dist min_obs_time = photons[0].obs_time max_obs_time = photons[len(photons)-1].obs_time min_obs_dist = sqrt(photons[0].pos.x*photons[0].pos.x + photons[0].pos.y*photons[0].pos.y) max_obs_dist = sqrt(photons[len(photons)-1].pos.x*photons[len(photons)-1].pos.x + photons[len(photons)-1].pos.y*photons[len(photons)-1].pos.y) actual_speed = (max_emit_dist-min_emit_dist)/(max_emit_time-min_emit_time) print "Actual speed = %.2f" % actual_speed app_speed = (max_obs_dist-min_obs_dist)/(max_obs_time-min_obs_time) print "Apparent speed = %.2f" % app_speed ################################################ # # Put a couple of lines on both the emitted # and observed distance versus time data # to show the limits of the motion # ################################################ funct1 = gcurve(color=color.blue) funct1.plot(pos=(0,min_emit_dist)) funct1.plot(pos=(min_emit_time,min_emit_dist)) funct1.plot(pos=(min_emit_time,0)) funct2 = gcurve(color=color.blue) funct2.plot(pos=(0,max_emit_dist)) funct2.plot(pos=(max_emit_time,max_emit_dist)) funct2.plot(pos=(max_emit_time,0)) label(display=emit_graph.display,pos=(min_emit_time*0.9,(max_emit_dist-min_emit_dist)/1.6), text="%.2f light-years" % (max_emit_dist-min_emit_dist)) label(display=emit_graph.display,pos=((max_emit_time-min_emit_time)/2.0+min_emit_time,-y_plot_max*0.1), text="%.2f years" % (max_emit_time-min_emit_time)) funct3 = gcurve(color=color.green) funct3.plot(pos=(min_obs_time*0.8,min_obs_dist)) funct3.plot(pos=(min_obs_time,min_obs_dist)) funct3.plot(pos=(min_obs_time,0)) funct3 = gcurve(color=color.green) funct3.plot(pos=(min_obs_time*0.8,max_obs_dist)) funct3.plot(pos=(max_obs_time,max_obs_dist)) funct3.plot(pos=(max_obs_time,0)) label(display=emit_graph.display,pos=(min_obs_time*0.8,(max_obs_dist-min_obs_dist)/1.6), text="%.2f light-years" % (max_obs_dist-min_obs_dist)) label(display=emit_graph.display,pos=((max_obs_time-min_obs_time)/2.0+min_obs_time,-y_plot_max*0.1), text="%.2f years" % (max_obs_time-min_obs_time)) ################################################ # # Add labels reporting the speeds # ################################################ label(display=emit_graph.display,pos=(x_plot_max*0.3,y_plot_max*0.8*extra_offset_factor), text="Actual Speed = %.2f x light speed" % actual_speed, height=15, color=color.blue, box=0) label(display=emit_graph.display,pos=(x_plot_max*0.8,y_plot_max*0.8*extra_offset_factor), text="Apparent Speed = %.2f x light speed" % app_speed, height=15, color=color.green, box=0) #################################################### # Grab a final movie frame if we are making a movie #################################################### if make_movie: GetScreenShot(FrameNumber) Main()