from numpy import * import matplotlib.pyplot as plt # This code computes the apparent travel time between two galaxies in an accelerating universe. # All calculations are done in terms of a_turnaround, when the ship decelerates. # The assumption is that the ship accelerates at g for some time, and decelerates at -g thereafter. # All units in this problem are in terms of [yr]^n, with c=1 h=0.68 H0=h/9.78e9 # inverse years g=1.03 # inverse years gamc=g/H0 a0=1 def v(gamma): return sqrt(1-1/pow(gamma,2)) def gettau(a1): a=1 Nbins=1000 damax=(a1-1)/Nbins tau=0 if (a1 < 1.001): tau=2.*log(1+g*x/2+sqrt(g*x/2*(g*x/2+2)))/g elif (a1 < e): #print ("Terminal velocity not reached") while(a < a1): da=damax dtau1=1./(a*a*sqrt(1+pow(gamc*log(a),2))) af=a+da dtau2=1./(af*af*sqrt(1+pow(gamc*log(af),2))) while ((dtau1-dtau2)/(dtau1+dtau2) >1e-2): da=da/2. af=a+da dtau2=1./(af*af*sqrt(1+pow(gamc*log(af),2))) tau=tau+da*(dtau1+dtau2)/2. a=a+da tau=2*tau/H0 else: #print ("Terminal Velocity is reached") while(a < exp(1)): da=damax dtau1=1./(a*a*sqrt(1+pow(gamc*log(a),2))) af=a+da dtau2=1./(af*af*sqrt(1+pow(gamc*log(af),2))) while ((dtau1-dtau2)/(dtau1+dtau2) >1e-2): da=da/2. af=a+da dtau2=1./(af*af*sqrt(1+pow(gamc*log(af),2))) tau=tau+da*(dtau1+dtau2)/2. a=a+da tau=2*tau/H0+log(a2/exp(2))/g return tau #loop over a2 (and expansion) da=1e-12 xarr=zeros(100) tarr=zeros(100) tauarr=zeros(100) i=0 for i in range(100): a2=1+da a1=sqrt(a2) x=1./H0*(1-1./a2) if (x < 100): t=2*sqrt(x/2*(g*x/2+2)/g) else: t=-1/H0*log(1-x*H0) tau=gettau(a1) #print(a2,x,t,tau) xarr[i]=x tarr[i]=t tauarr[i]=tau da=da*1.5 i=i+1 plt.loglog(xarr,tarr,lw=2,color='red') plt.loglog(xarr,tauarr,lw=3,color='blue') plt.xlabel('Distance in Lightyears (Now)') plt.ylabel('Time for Trip in Years') plt.title("Travel Time in an Accelerating Universe") plt.annotate('Universe Clock',xy=(xarr[40],tarr[40]),xytext=(xarr[20],100*tarr[40]),arrowprops=dict(facecolor='red',shrink=0.1)) plt.annotate('Rocket Clock',xy=(xarr[50],tauarr[50]),xytext=(xarr[60],1000*tauarr[60]),arrowprops=dict(facecolor='blue',shrink=0.1)) plt.show()