/** 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)] }