MATLAB Examples


function [selftime, maxspeed, distance, selfdist] = twins(t1, a1, t2, v2, t3, a3)

Calculate characteristics of relativistic rocket

The rocket trajectory consists of three phases - acceleration, uniform movement and deceleration.



  • t1 - duration of acceleration phase (years)
  • a1 - value of acceleration ($ly/y^2$)
  • t2 - duration of uniform movement phase (years)
  • v2 - speed of uniform movement (fraction of light speed). Only used when t1 is zero
  • t3 - duration of deceleration phase (years)
  • a3 - value of deceleration ($ly/y^2$)


  • selftime - self time of rocket (traveller's proper time) (years)
  • maxspeed - maximum speed reached by the end of acceleration phase (fraction of light speed)
  • distance - distance travelled from earth's POV (light years)
  • selfdist - distance travelled from traveller's POV (light years)


  • Initial speed is zero (see another formula where v0 is considered)
  • c = 1 to make formulas tidy (then acceleration is $ly/y^2$ and speed is in light-years)


  • Results are best viewed with format longG
  • Acceleration of 1 $ly/y^2$ is approximately 9.5 $m/s^2$. 1.03 $ly/y^2$ is 1g (9.81 $m/s^2$)


  • Travel to Proxima-Centaury, a star about 4.223 light years away, at a speed of 80% that of light. No acceleration and deceleration, just uniform movement.
[selftime, maxspeed, distance, selfdist] = twins(0,0,4.223/0.8,0.8,0,0);
  • Travel to Alpha-Centaury, a star about 4.3 light years away. Accelerate half of way and decelerate half of way with 9.5 m/s^2 (no uniform movement).
[selftime, maxspeed, distance, selfdist] = twins(3,1,0,0,3,1);

Define some useful constants

Speed of light in vacuum, $m/s$:

c = 299792458; %#ok

Convert light year to meters (c*365.25*24*3600)

lightyear_m = 9460730472580800;

Seconds in a year (365.25*24*3600)

seconds_per_year = 31557600;

Convert acceleration from $ly/y^2$ to $m/s^2$

acc_mss = lightyear_m/seconds_per_year^2;

Sane check inputs

if (t1 == 0)
    % Don't account for acceleration if it is applied during zero time duration
    a1 = eps;
    % Don't calculate maxspeed if there was no acceleration.
    % Use user-supplied uniform speed in calculations (not greater than 1c).
    maxspeed = min(1,v2);
    % Calculate maximum speed by the end of acceleration
    maxspeed = (a1*t1)/sqrt(1+(a1*t1)^2);
if (t3 == 0)
    % Don't account for deceleration if it is applied during zero time duration
    a3 = eps;

Main calculations


  • Speed reached by the end of acceleration (calculated above)

$maxspeed = \frac{a1*t1}{\sqrt{1+(a1*t1)^2}}$

  • Rocket proper time

$selftime = \frac{1}{a1}asinh(a1*t1) + t2\sqrt{1-maxspeed^2} + \frac{1}{a3}asinh(a3*t3)$

selftime = 1/a1*asinh(a1*t1) + t2*(sqrt(1-maxspeed^2)) + 1/a3*asinh(a3*t3);
  • Distance travelled in Earth's frame

Account initial speed when deceleration starts:

$\pi_0 = \frac{maxspeed}{\sqrt{1-maxspeed^2}}$

pi0 = maxspeed/sqrt(1-maxspeed^2);
% Speed decrease due to deceleration
% decspeed = (pi0-a3*t3)/sqrt(1+(pi0-a3*t3)^2);

Distance travelled during acceleration phase:

$distance\_acc = \frac{1}{a1}(\sqrt{1+(a1*t1)^2} - 1)$

distance_acc = 1/a1*(sqrt(1+(a1*t1)^2) - 1);

Distance travelled during uniform movement phase:

$distance\_uni = t2*maxspeed$

distance_uni = t2*maxspeed;

Distance travelled during deceleration phase:

$distance\_dec = -\frac{1}{a3}(\sqrt{1+(\pi_0-a3*t3)^2} - \sqrt{1+\pi_0^2})$

distance_dec = -1/a3*(sqrt(1+(pi0-a3*t3)^2) - sqrt(1+pi0^2));

Total distance:

$distance = distance\_acc + distance\_uni + distance\_dec$

distance = distance_acc + distance_uni + distance_dec;
  • Rocket's proper distance

During accelerated phase length contraction is dependant on acceleration. During uniform, it's a function of speed.

$selfdist = \frac{distance\_acc}{\sqrt{1+{a1}^2}} + distance\_uni\sqrt{1-maxspeed^2} + \frac{distance\_dec}{\sqrt{1+{a3}^2}}$

selfdist = distance_acc/sqrt(1+a1^2)...
         + distance_uni*sqrt(1-maxspeed^2)...
         + distance_dec/sqrt(1+a3^2);

Plot results

% Acceleration phase
tacc = linspace(0,t1);
% Deceleration phase
tdec0 = linspace(t1+t2,t1+t2+t3);
tdec = linspace(0,t3);
xlabel('Time in earth frame, Years');
ylabel('Speed profile, Fraction of Light Speed');
legend('Acceleration','Uniform movement','Deceleration','Location','South');
out  = sprintf('Earth distance and time: {\\bf% 11.2f} ly in{\\bf% 11.2f} years',distance,t1+t2+t3);
out1 = sprintf('Rocket distance and time:{\\bf% 11.2f} ly in{\\bf% 11.2f} years',selfdist,selftime);
out2 = '';
out3 = sprintf('Speed by the end of acceleration:{\\bf% 24.17f} c',maxspeed);
out4 = sprintf('Acceleration: {\\bf%.2f} m/s^2, Deceleration: {\\bf%.2f} m/s^2',a1*acc_mss,a3*acc_mss);