Download or view relativity.frink in plain text format
/** This contains functions for relativity-related issues, including
black holes, relativistic travel, etc.
Good resources:
The Relativistic Rocket by John Baez.
http://math.ucr.edu/home/baez/physics/Relativity/SR/Rocket/rocket.html
Ned Wright's redshift resources and FAQ:
https://www.astro.ucla.edu/~wright/CosmoCalc.html
TODO: Implement his cosmological redshift calculations.
*/
/** This is the Lorentz factor, often known as gamma[v].
It represents several transformations, such as length contraction or mass
increase when traveling close to the speed of light. The result is
a dimensionless number from 1 (at zero velocity) to infinity (at
lightspeed) that represents the contraction.
The simple equation is
1 / sqrt[1 - v^2/c^2]
However, for low velocities, the naive evaluation may be indistinguishable
from 1. For correct results, we extend the precision of the
calculation to give a meaningful answer. This could be done by the Taylor
series
1 + 1/2 c^-2 v^2 + 3/8 c^-4 v^4 + 5/16 c^-6 v^6 + 35/128 c^-8 v^8 ...
or by using an arbitrary-precision sqrt (which is how we do it.)
*/
Lorentz[v is velocity] :=
{
v = abs[v]
if v > 1000 m/s
return 1 / sqrt[1 - v^2/c^2]
else
{
if v == 0 m/s
return 1
// Use adaptive precision.
oldPrec = getPrecision[]
try
{
newPrec = 25 + max[5, round[-2 log[v / (m/s)]]]
setPrecision[newPrec]
return 1 / sqrt[1 - v^2/c^2, newPrec]
}
finally
{
setPrecision[oldPrec]
}
}
}
/** This is the inverse of the Lorentz factor. It takes
a dimensionless number from 1 (at zero velocity) to infinity (at
lightspeed) representing the contraction factor and returns a velocity.
The simple equation is:
c sqrt[1 - 1/gamma^2]
but d might be very close to 1 for low velocities, in which case we need
to extend precision to get a meaningful result. We can either do this
with the Taylor series or with arbitrary-precision square root, which is
the tactic that this is taking.
*/
inverseLorentz[gamma is dimensionless] :=
{
if gamma == 1
return 0 m/s
//println["gamma is $gamma"]
gammaMinus1 = gamma-1
if gammaMinus1 == 0 // We've lost all significant figures
loggamma = -100
else
loggamma = log[gammaMinus1]
// println["loggamma is $loggamma"]
if loggamma > 1e-11
return c sqrt[1 - 1/gamma^2] // Normal calculation
else
{
// Lorentz gamma is very close to 1; use adaptive precision.
oldPrec = getPrecision[]
try
{
newPrec = max[17, round[-2 loggamma]]
setPrecision[newPrec]
// println["new precision is $newPrec"]
r = c sqrt[1 - 1/gamma^2, newPrec]
if newPrec >= 30
{
setPrecision[newPrec div 2]
return 1. r
}
return r
}
finally
{
setPrecision[oldPrec]
}
}
}
/** This is the integral of the Lorentz transformation from v0 to v1.
Not quite sure how this is interpreted yet. */
LorentzIntegral[v0 is velocity, v1 is velocity] := c (arcsin[v1/c] - arcsin[v0/c])
/** Relativistic Doppler ratio. This returns a dimensionless number by which
the source frequency is to be multiplied to get the observed frequency.
v is positive when the source is approaching.
If v is very low, we need to extend precision to get a meaningful result.
We can either do this with the Taylor series or with arbitrary-precision
square root.
Taylor series is:
b = v/c
sum = b + b^2 / 2 + b^3 / 2 + 3 b^4 / 8 + 3 b^5 / 8 + 5 b^6 / 16 + 5 b^7 / 16 + 35 b^8 / 128 + 35 b^9 / 128 + 63 b^10 / 256 + 63 b^11 / 256
*/
relativisticDopplerRatio[v is velocity] :=
{
if v == 0 m/s
return 1
if abs[v] > 1000 m/s
return sqrt[ (1 + v/c) / (1-v/c) ] // Normal case
else // Very slow, use Taylor Series
{
b = v/c
digits = ceil[- log[b]] + 18
// println["digits is $digits"]
orig = getPrecision[]
try
{
setPrecision[digits]
return sqrt[(1 + v/c) / (1-v/c), digits]
}
finally
{
setPrecision[orig]
}
}
}
/** Inverse relativistic Doppler ratio. This takes a dimensionless number
by which the source frequency is to be multiplied to get the observed
frequency and returns a velocity. The velocity is positive when
approaching. */
inverseRelativisticDopplerRatio[r is dimensionless] :=
{
c (r^2 - 1) / (r^2 + 1)
}
/** Relativistic Z factor. The z factor is the ratio by which a frequency
is shifted due to movement. Z = deltaLambda / lambda.
See:
http://hyperphysics.phy-astr.gsu.edu/hbase/Astro/redshf.html#c2
*/
relativisticZ[v is velocity] :=
{
sqrt[ ( 1 + v/c ) / ( 1 - v/c) ] - 1
}
/** Inverse relativistic Z factor.
Returns: a velocity
See:
http://hyperphysics.phy-astr.gsu.edu/hbase/Astro/redshf.html#c2
*/
inverseRelativisticZ[z is dimensionless] :=
{
(c z^2 + 2 c z) / (z^2 + 2z + 2)
}
/** Relativistic velocity addition. This is the simplest case where the motions
lie along a line.
A B -> v C ->
where v is the relative velocity between A (assumed stationary) and B.
B launches a projectile (C).
u is the velocity of the projectile (C) as seen by A
uprime is the velocity of the projectile (C) as seen by B.
To use this function, pass *TWO* of the arguments to this function, leaving
the other as the special value undef. It will solve for the undefined
variable and return the results in the same order [v, u, uprime] with the
unspecified value filled in.
For example,
[v, u, uprime] = addVelocity[.5 c, undef, .5 c]
println[u -> "c"]
returns:
[v, u, uprime]
*/
addVelocity[v, u, uprime] :=
{
if uprime == undef
uprime = (u - v) / (1 - u v / c^2)
else
if v == undef
v = (u - uprime) / (1 - u uprime / c^2)
else
if u == undef
u = (uprime + v) / (1 + uprime v / c^2)
return [v, u, uprime]
}
/** Relativistic rocket equations. See
"The Relativistic Rocket" by John Baez.
http://math.ucr.edu/home/baez/physics/Relativity/SR/Rocket/rocket.html
and Gravitation, Misner, Thorne, and Wheeler, Section 6.2.
This defines the equations using the following conventions:
"First of all, we need to be clear what we mean by continuous acceleration
at 1g. The acceleration of the rocket must be measured at any given
instant in a non-accelerating frame of reference travelling at the same
instantaneous speed as the rocket (see the FAQ article on accelerating
clocks, http://math.ucr.edu/home/baez/physics/Relativity/SR/clock.html
), because this is the acceleration that its occupants would physically
feel—and we want them to accelerate at a comfortable rate that has the
effect of mimicking their weight on Earth. We'll call this acceleration a.
The proper time measured by the crew of the rocket (i.e. how much they age)
is called T, and the time measured in the non-accelerating frame of
reference in which they started (e.g. Earth) is t. We assume that the
stars are essentially at rest in this Earth frame. The distance covered by
the rocket as measured in this frame of reference is d, and the rocket's
velocity is v. The time dilation or length-contraction factor at any
instant is the gamma factor gamma.
*/
/** Relativistic rocket. For a constant acceleration, calculate the "fixed"
(Earth) time elapsed given the time experienced by the astronauts.
args:
[a, properTime]
where
a is acceleration (as felt by people in the rocket)
properTime is time (T, the time experienced on board.)
returns:
time elapsed on earth (t)
(This is the first half of Baez eq. 4)
*/
rocketEarthTimeFromProperTime[a is acceleration, properTime is time] :=
{
c / a sinh[a properTime / c]
}
/** Relativistic rocket. For a constant acceleration, calculate the "fixed"
(Earth) time elapsed given the time experienced by the astronauts.
args:
[a, d]
where
a is acceleration (as felt by people in the (rocket)
d is distance (as measured on earth)
returns:
time elapsed on earth (t)
(This is the second half of Baez eq. 4)
*/
rocketEarthTimeFromDistance[a is acceleration, d is length] :=
{
sqrt[(d/c)^2 + 2d/a]
}
/** Relativistic rocket. For a constant acceleration, calculate the "proper"
time experienced by the astronauts given the time elapsed on Earth.
args:
[a, earthTime]
where
a is acceleration (as felt by people in the rocket)
earthTime is time (t, the time experienced on earth.)
returns:
"proper" time experienced by astronauts. (T)
(This is the first half of Baez eq. 5)
*/
rocketProperTimeFromEarthTime[a is acceleration, earthTime is time] :=
{
c / a arcsinh[a earthTime / c]
}
/** Relativistic rocket. For a constant acceleration, calculate the "proper"
time (that is, the time experienced by the astronauts) given the distance
traveled (as seen from Earth.)
args:
[a, d]
where
a is acceleration (as experienced by people in the rocket)
d is distance (as seen from Earth)
returns:
"proper" time experienced by astronauts. (T)
(This is the second half of Baez eq. 5)
*/
rocketProperTimeFromDistance[a is acceleration, d is length] :=
{
c / a arccosh[a d / c^2 + 1]
}
/** Relativistic rocket. For a constant acceleration, calculate the distance
traveled (as seen from Earth) given the "proper" time, that is, the time
as experienced by the astronauts.
args:
[a, properTime]
where
a is acceleration (as experienced by people in the rocket)
properTime is time (T, the time experienced on board.)
returns:
distance traveled, as measured from Earth (d)
(This is the first half of Baez eq. 6)
*/
rocketDistanceFromProperTime[a is acceleration, properTime is time] :=
{
c^2 / a (cosh[a properTime / c] - 1)
}
/** Relativistic rocket. For a constant acceleration, calculate the distance
(as seen from Earth) given the time elapsed on earth.
args:
[a, earthTime]
where
a is acceleration (as experienced by people in the rocket)
earthTime is time (t)
returns:
distance, as measured from Earth (d)
(This is the second half of Baez eq. 6)
*/
rocketDistanceFromEarthTime[a is acceleration, earthTime is time] :=
{
c^2 / a (sqrt[1 + (a earthTime / c)^2] - 1)
}
/** Relativistic rocket. For a constant acceleration, calculate the velocity as
seen from earth, given the "proper time", that is, the time experienced by
the astronauts.
args:
[a, properTime]
where
a is acceleration (as experienced by people in the rocket)
properTime is time (T, the time experienced on board.)
returns:
velocity, as measured from Earth (v)
(This is the first half of Baez eq. 7)
*/
rocketVelocityFromProperTime[a is acceleration, properTime is time] :=
{
c tanh[a properTime / c]
}
/** Relativistic rocket. For a constant acceleration, calculate the velocity as
seen from earth, given the time elapsed on earth.
args:
[a, earthTime]
where
a is acceleration (as experienced by people in the rocket)
earthTime is time (t)
returns:
velocity, as measured from Earth (v)
(This is the second half of Baez eq. 7)
*/
rocketVelocityFromEarthTime[a is acceleration, earthTime is time] :=
{
a earthTime / sqrt[1 + (a earthTime/c)^2]
}
/** Relativistic rocket. For a constant acceleration, calculate the Lorentz
factor gamma given the "proper time", that is, the time experienced by
the astronauts.
args:
[a, properTime]
where
a is acceleration (as experienced by people in the rocket)
properTime is time (T, the time experienced on board.)
returns:
Lorentz factor gamma
(This is the first third of Baez eq. 8)
*/
rocketGammaFromProperTime[a is acceleration, properTime is time] :=
{
cosh[a properTime / c]
}
/** Relativistic rocket. For a constant acceleration, calculate the Lorentz
factor gamma, given the time elapsed on earth.
args:
[a, earthTime]
where
a is acceleration (as experienced by people in the rocket)
earthTime is time (t)
returns:
Lorentz factor gamma
(This is the second third of Baez eq. 8)
*/
rocketGammaFromEarthTime[a is acceleration, earthTime is time] :=
{
sqrt[1 + (a earthTime / c)^2]
}
/** Relativistic rocket. For a constant acceleration, calculate the Lorentz
factor gamma, given the distance traveled as seen from Earth.
args:
[a, d]
where
a is acceleration (as experienced by people in the rocket)
d is distance (as seen from Earth)
returns:
Lorentz factor gamma
(This is the third third of Baez eq. 8)
*/
rocketGammaFromDistance[a is acceleration, d is length] :=
{
a d / c^2 + 1
}
/** Relativistic rocket. For a constant acceleration, calculate the "proper
time", that is, the time experienced by the astronauts, required to travel
to a distant object, and *stop*.
args:
[a, d]
where
a is acceleration (as experienced by people in the rocket)
d is distance (as measured from earth)
returns:
proper time T experienced by the astronauts.
*/
rocketProperTimeToStop[a is acceleration, d is length] :=
{
return 2 rocketProperTimeFromDistance[a, d/2]
}
/** Relativistic rocket. For a constant acceleration, calculate the time as
seen from earth required to travel to a distant object, and *stop*.
args:
[a, d]
where
a is acceleration (as experienced by people in the rocket)
d is distance (as measured from earth)
returns:
earth time t
*/
rocketEarthTimeToStop[a is acceleration, d is length] :=
{
return 2 rocketEarthTimeFromDistance[a, d/2]
}
/** Relativistic rocket. This is the relativistic equivalent of the Tsiolkovsky
rocket equation also known as the "ideal rocket equation".
This calculates:
deltaV = c tanh[(vexhaust / c) ln[massratio]]
Compare this to the non-relativistic version of the ideal rocket equation:
deltaV = vexhaust ln[massratio]
To compare this to a formulation where the exhaust velocity is not known,
but the specific impulse is known and specified in seconds (ugh,) use the
equivalence
vexhaust = Isp gravity
where gravity is the standard acceleration of gravity (9.80665 m/s^2, which
is known to Frink as the unit "gravity")
args:
[vexhaust, massratio]
where
vexhaust is the exhaust velocity
massratio is the ratio between the total initial mass (including payload)
and payload mass (masstotal / masspayload)
returns:
the deltaV (change in velocity) experienced by the final payload.
*/
rocketIdealRocketEquation[vexhaust is velocity, massratio is dimensionless] :=
{
c tanh[(vexhaust / c) ln[massratio]]
}
/** This section contains black hole calculations. Much of it is based on
the nice resource (collected in one place) here:
https://www.vttoth.com/CMS/physics-notes/311-hawking-radiation-calculator
*/
/** Calculate the radius of a black hole given its mass. */
blackHoleRadius[M is mass] := 2 M G / c^2
/** Calculate the mass of a black hole given its radius. */
blackHoleMassFromRadius[r is length] := 1/2 r c^2 / G
/** Calculate the surface area of a black hole given its mass. */
blackHoleArea[M is mass] := M^2 16 pi G^2 / c^4
/** Calculate the mass of a black hole given its surface area. */
blackHoleMassFromArea[a is area] := 1/4 G^-1 a^(1/2) c^2 pi^(-1/2)
/** Calculate the "surface" gravity of a black hole given its mass. */
blackHoleGravity[M is mass] := c^4 / (4 G M)
/** Calculate the mass of a black hole given its surface gravity. */
blackHoleMassFromGravity[a is acceleration] := 1/4 c^4 / (a G)
/** Calculate the tides at the surface of a black hole given its mass. The
results have dimensions of acceleration/distance, that is, how much
gravity changes as you move away. */
blackHoleTides[M is mass] := c^6 / (M^2 4 G^2)
/** Calculate the mass of a black hole given tides at its surface. The
argument has dimensions of acceleration/distance or s^-2 */
blackHoleMassFromTide[tide] := 1/2 c^3 / (G sqrt[tide])
/** Calculate the (dimensionless) Bekenstein-Hawking entropy of a black hole
given its mass. To get the "chemist's" entropy, this is multiplied by
the Boltzmann constant k. */
blackHoleDimensionlessEntropy[M is mass] := M^2 4 pi G / (hbar c)
/** Calculate the "chemist's" entropy of a black hole
given its mass. This is just the dimensionless entropy multiplied by
the Boltzmann constant k. */
blackHoleEntropy[M is mass] := k M^2 4 pi G / (hbar c)
/** Calculate the black hole's temperature given its mass. */
blackHoleTemperature[M is mass] := hbar c^3 / (M 8 pi k G)
/** Calculate the black hole's mass given its temperature. */
blackHoleMassFromTemperature[T is temperature] := 1/8 c^3 hbar / (G T k pi)
/** Calculate the black hole's luminosity (radiated power) given its mass. */
blackHolePower[M is mass] := hbar c^6 / (M^2 15360 pi G^2)
/** Calclulate the black hole's mass given its radiated power. */
blackHoleMassFromPower[p is power] := c^3 sqrt[hbar / (15 pi p)] / (32 G)
/** Calculate a black hole's lifetime given its mass. */
blackHoleLifetime[M is mass] := M^3 5120 pi G^2 / (hbar c^4)
/** Calculate a black hole's mass given its lifetime. */
blackHoleMassFromLifetime[t is time] := (c^4 hbar t / (5120 G^2 pi))^(1/3)
/** Calculate a black hole's effective density given its mass. */
blackHoleDensity[M is mass] := 3 c^6 / (32 pi G^3 M^2)
/** Calculate a black hole's mass given its effective density. */
blackHoleMassFromDensity[d is mass_density] := sqrt[3 c^6 / (32 d G^3 pi)]
/** Calculate the maximum proper time it takes to fall to the singularity from
the Schwarzschild radius. The use of any rocket decreases this time! */
blackHoleTimeToSingularity[M is mass] := pi G M / c^3
/** Calculate the mass of a black hole given the time it takes to fall to the
singularity. */
blackHoleMassFromTimeToSingularity[t is time] := c^3 t / (pi G)
/** Relativistic escape velocity for an object with a given mass and radius.
This approximates the Newtonian escape velocity which is
v = sqrt[2 G M / r]
See:
http://www.mrelativity.net/MBriefs/Relativistic%20Escape%20Velocity%20using%20Special%20Relativity.htm
args:
[mass, radius]
*/
relativisticEscapeVelocity[M is mass, r is length] :=
{
phi = G M / r // Gravitational specific potential energy
return sqrt[2 phi - (phi / c)^2]
}
/** Relativistic gravitational potential energy relative to an object (planet,
star, black hole, with mass M) and a smaller object with rest mass m0 at
radial distance r from the more massive object, traveling at velocity v.
If the velocity is zero, this turns into the Newtonian potential energy
which is given by:
PE = - G M m0 / r
See:
http://www.mrelativity.net/MBriefs/Relativistic%20Escape%20Velocity%20using%20Special%20Relativity.htm
args:
[M, m0, r, v=0 m/s]
*/
relativisticPotentialEnergy[M is mass, m0 is mass, r is length, v is velocity = 0 m/s] :=
{
- G M m0 / (r sqrt[1 - v^2/c^2])
}
/** Relativistic kinetic energy for an object with rest mass m0 at
traveling at velocity v. This is valid for all velocities but has numerical
issues if evaluated naively for low velocities in the given form
(m0 c^2) (Lorentz[v] - 1)
so, for low velocities, we expand it to a Taylor series which is closer to
the classical 1/2 m v^2 (which you can actually see in the equation) but
with interesting correction terms..
This would be better with an arbitrary-precision square root routine like
in sqrtWayne.frink. The sample program RelativisticKineticEnergyTest.frink
tests the boundary where the naive algorithm and the Taylor series should
be used (in the absence of an arbitrary-precision square root algorithm,
which is what should really be used.)
As of the 2020-04-01 release of Frink, there is now an arbitrary-precision
sqrt routine sqrt[n, digits] guaranteed to be always available in both Frink
and Frink:The Next Generation (although it will be orders of magnitude faster
in the latter) and the Taylor series can probably be retired, although it's
very interesting and shows how even the simplest relativistic calculations
can quickly split into an infinite series of higher-order terms, or can be
approximated pretty well by a single low-order term in the classical
approximation.
args:
[m0, v]
where:
m0: The rest mass of the object
v: The velocity
This could be modified to omit the Taylor series and use an
arbitrary-precision sqrt function using something like:
digits = ceil[-log[v/c]]
digits = max[18, digits]
oldprec = getPrecision[]
try
{
setPrecision[digits]
return (m0 c^2) (1/sqrt[1-v^2/c^2, digits] - 1)
}
finally
setPrecision[oldprec]
See:
http://www.mrelativity.net/MBriefs/Relativistic%20Escape%20Velocity%20using%20Special%20Relativity.htm
*/
relativisticKineticEnergy[m0 is mass, v is velocity] :=
{
if v < 0.13 c // Threshold found empirically, use Taylor series
return 1. * m0 ( (v^2 / 2) + (3 v^4 / (8 c^2)) + (5 v^6 / (16 c^4)) + (35 v^8 / (128 c^6)) + (63 v^10 / (256 c^8)) + (231 v^12 / (1024 c^10)) + (429 v^14 / (2048 c^12)) + (6435 v^16 / (32768 c^14)) )
else
return (m0 c^2) (1/sqrt[1-v^2/c^2] - 1)
}
/** Calculates gravitational time dilation factor for an object *at rest* at
radius r from a massive (nonrotating) spherical object. This factor is
less than 1 and represents the slower rate of flow of time for the object
inside the gravitational field with respect to an object at infinity. */
timeDilation[M is mass, r is length] :=
{
return sqrt[1 - 2 G M / (r c^2)]
}
/** Calculates gravitational *and* speed-related time dilation factor for an
object *in orbit* at radius r from a massive (nonrotating) spherical
object. This combines the effects of both gravity speedup from being far
from the object and special relativistic clock slowdown due to orbital
speed. Newtonian orbital velocity is vo = sqrt[G M / r]
See the graph at
https://en.wikipedia.org/wiki/Gravitational_time_dilation#Experimental_confirmation
You can reproduce the blue "net orbital time gain" curve of that graph for
a given height h above earth's center with:
orbitalTimeDilation[earthmass, h] - timeDilation[earthmass, earthradius] -> ps/s
*/
orbitalTimeDilation[M is mass, r is length] :=
{
return sqrt[1 - 3 G M / (r c^2)]
}
Download or view relativity.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19971 days, 9 hours, 5 minutes ago.