/** 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.htm
l
    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 gamma 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 is now used.
*/
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 low velocity, use arbitrary-precision sqrt
   {
      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]
      }
   }
}

/** Relativistic Z factor.   The z factor is the ratio by which a frequency
    is shifted due to movement.  Z = deltaLambda / lambda.

    The simplified equation is
       sqrt[ ( 1 + v/c ) / ( 1 - v/c) ] - 1

    This is the same as calling
       relativisticDopplerRatio[v is velocity] - 1
    which generates more digits for small velocities.

   See:
    http://hyperphysics.phy-astr.gsu.edu/hbase/Astro/redshf.html#c2

   TODO:  Make this return adequate digits for v very close to 0.
*/
relativisticZ[v is velocity] :=
{
   sqrt[ ( 1 + v/c ) / ( 1 - v/c) ] - 1
}

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


/** Inverse relativistic Z factor.  The z factor is the ratio by which a
    frequency is shifted due to movement.  Z = deltaLambda / lambda.

    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)

    TODO:  Improve Frink's implementation of arcsinh around 0 to give better
    results for small arguments, or use a Taylor series for arcsinh which is
    arcsinh[x]=x + -1/6 x^3 + 3/40 x^5 + -5/112 x^7 + 35/1152 x^9 + -63/2816 x^11...
    This means we probably need 3 times as many significant figures as we get
    small.
*/
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)

    This works hard to make this return enough digits for very small arguments.
    This should  approximate the Newtonian limit
      d = 1/2 a t^2
    for small accelerations and times.

    Alternately,

    cosh[x] - 1 = 2 sinh[x/2]^2

    which is the tactic used here to get more correct digits for small
    velocities.
*/
rocketDistanceFromProperTime[a is acceleration, properTime is time] :=
{
   // pseudovelocity as fraction of c (can be larger than c)
   v = a properTime / c

   if v > 0.01
      return c^2 / a (cosh[v] - 1)       // Normal direct calculation
   else
      return c^2 / a (2 sinh[v/2]^2)     // Low velocity trick
}


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

    TODO:  This should approximate the Newtonian result as d = 1/2 a t^2 for
    small velocities.   The Taylor series of sqrt[1+v^2] - 1 is:

    sqrt(1+x)-1 = (1/2)*x - (1/8)*x^2 + (1/16)*x^3 - (5/128)*x^4 + (7/256)*x^5 - (21/1024)*x^6 + (33/2048)*x^7 + ...
 
    Coefficients are 1/2, -1/8, 1/16, -5/128, 7/256, -21/1024, 33/2048, -429/32768, 715/65536, -2431/262144, 4199/524288, -29393/4194304, 52003/8388608, ...

    Se OEIS https://oeis.org/A002596

    This currently just uses an arbitrary-precision sqrt with lots of digits.
*/
rocketDistanceFromEarthTime[a is acceleration, earthTime is time] :=
{
   origPrec = getPrecision[]
   try
   {
      workPrec = 110
      setPrecision[workPrec]
      x = sqrt[1 + (a earthTime / c)^2, workPrec] - 1
      setPrecision[origPrec]
      return x c^2 / a 
   }
   finally
      setPrecision[origPrec]
}


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

    TODO:  Make this return enough digits for very small velocities.  Do we
    need an arbitrary-precision cosh function?  Or expand the Taylor series?
    (ArbitraryPrecision.frink now has an arbitraryCosh function.)

    The direct equation is simply:
    cosh[a properTime / c]
    
    The Taylor series of cosh is:
    1 + 1/2 x^2 + 1/24 x^4 + 1/720 x^6 + 1/40320 x^8 + 1/3628800 x^10 + 1/479001600 x^12 + 1/87178291200 x^14 + 1/20922789888000 x^16

    The denominator is (2x)!      https://oeis.org/A010050

    cosh[x] = 1 + x^2/2! + x^4/4! + x^6/6! + ...

    For small velocities, cosh[x] is a number barely above 1 (but always above)
    so we have to potentially have *many* significant digits to return a
    result that is not just (incorrectly) 1.0.
    
    Alternately,

    cosh[x] - 1 = 2 sinh[x/2]^2

    so you could calculate using sinh and just add 1 to the result.  Which is
    probably the easiest way to preserve significant digits.  This
    necessitated fixing Frink's native sinh function to use a Taylor series
    for small x, as it lost precision and always returned 0 for small
    velocities.  This is fixed in the 2026-10-16 release.

    if v = a properTime / c
    It appears that you need at least (2 * -log[v] + d) digits to get d
    significant figures.
*/
rocketGammaFromProperTime[a is acceleration, properTime is time] :=
{
   // pseudovelocity as fraction of c (can be larger than c)
   v = a properTime / c
   // println["v is $v"]
   
   if v > 0.13
      return cosh[v]    // Direct calculation for higher pseudovelocity
   else
   {
      // Lower velocity which requires more digits to be calculated to give
      // a significant return value.
      //
      // Calculate using the trick:
      // cosh[x] = 1 + 2 sinh[x/2]^2
      origPrec = getPrecision[]
      digits = 17
      try
      {
         prec = max[17, ceil[2 * -log[v] + digits]]    // max is for safety
         //println["Setting precision to $prec"]
         setPrecision[prec]
         res = 2 sinh[v/2]^2
         //println["res is $res"]
         return  1. + res
      }
      finally
         setPrecision[origPrec]
   }
}


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

    The simplistic calculation is:
       sqrt[1 + (a earthTime / c)^2]

    but we work harder to make this return significant digits even for very
    small velocities.

    TODO:  Make this return enough digits for very small velocities
*/
rocketGammaFromEarthTime[a is acceleration, earthTime is time] :=
{
   // pseudovelocity as fraction of c (can be larger than c)
   v = a earthTime / c
   
   origPrec = getPrecision[]
   digits = 17        // Target number of actual significant digits
   prec = max[17, ceil[2 * -log[v] + digits]]    // For safety
   try
   {
      setPrecision[prec]
      return sqrt[1 + v^2, prec]
   }
   finally
      setPrecision[origPrec]
}


/** 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 given distance (as seen from earth) and time
    (as seen from earth,) calculate the necessary (constant) acceleration
    experienced by the astronauts to travel and pass by the given distance
    (not to stop.)

    This is the right-hand-side of Baez eq. 6 solved for a.

    args:
      [d, earthTime]
        where
      d is distance (as seen from Earth)
      earthTime is time (as seen from Earth)

    returns:
       the constant acceleration experienced by the astronauts.

    (This is the last part of Baez eq. 6 solved for a.)

    For low pseudovelocities, this should approximate the Newtonian value
       a = 2 d / t^2
*/
rocketAccelerationFromEarthTime[d is length, earthTime is time] :=
{
   if earthTime <sqrt[d^2 / c^2]
   {
      println["rocketAccelerationFromEarthTime:  impossible to complete in less than " + sqrt[d^2 / c^2] + " (time given was $earthTime)"]
      return undef
   }

   return (2 c^2 d)/(c^2 earthTime^2 - d^2)
}


/** Relativistic rocket.  For a given distance (as seen from earth) and time
    (as seen from earth,) calculate the necessary (constant) acceleration
    experienced by the astronauts to travel AND STOP AT the given distance.

    This is the right-hand-side of Baez eq. 6 solved for a.

    args:
      [d, earthTime]
        where
      d is distance (as seen from Earth)
      earthTime is time (as seen from Earth)

    returns:
       the constant acceleration experienced by the astronauts.

    (This is the last part of Baez eq. 6 solved for a.)

    This is the acceleration needed to travel half the distance in half the
    (earth) time.
*/
rocketAccelerationFromEarthTimeToStop[d is length, earthTime is time] :=
{
   return rocketAccelerationFromEarthTime[d/2, earthTime/2]
}


/** 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.   This is also known
    as the Schwarzschild radius. */
blackHoleRadius[M is mass] := 2 G M / 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.

    This takes effort to return meaningul digits even when the time dilation
    factor is very close to 1.  The simple equation without these efforts is
    just:

    sqrt[1 - 2 G M / (r c^2)]
*/
timeDilation[M is mass, r is length] :=
{
   arg1 = 2 G M / (r c^2)
   // println["arg1 is $arg1"]
   digits = ceil[max[-1.8 log[arg1], 20]]

   origPrec = getPrecision[]
   try
   {
      setPrecision[digits]
      arg = 1 - arg1
      // println["arg is $arg"]
      // println["digits is $digits"]
      return sqrt[arg, digits]
   }
   finally
   {
      setPrecision[origPrec]
   }
}

/** 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

    This takes effort to return meaningul digits even when the time dilation
    factor is very close to 1.  The simple equation without these efforts is
    just:

    sqrt[1 - 3 G M / (r c^2)]
*/
orbitalTimeDilation[M is mass, r is length] :=
{
   arg1 = 3 G M / (r c^2)
   // println["arg1 is $arg1"]
   digits = ceil[max[-1.8 log[arg1], 20]]

   origPrec = getPrecision[]
   try
   {
      setPrecision[digits]
      arg = 1 - arg1
      // println["arg is $arg"]
      // println["digits is $digits"]
      return sqrt[arg, digits]
   }
   finally
   {
      setPrecision[origPrec]
   }
}


/** Gravitational wave power emitted by two masses mass1 and mass2 in a circular
    orbit at the specified distance from each other.

    See:
      https://astronomy.swin.edu.au/cosmos/G/Gravitational+Waves
*/
gravitationalWavePower[mass1 is mass, mass2 is mass, distance is length] :=
{
   return 32 G^4 / (5 c^5 distance^5) (mass1 mass2)^2 (mass1 + mass2)
}

/** Gravitational-wave-based period derivative (the rate of change of the
    orbital period) of two masses mass1 and mass2 in an eccentric orbit with
    eccentricity ecc and the specified period.  The return value is a
    dimensionless value, dperiod/dt.

    See:
      https://astronomy.swin.edu.au/cosmos/G/Gravitational+Waves
      https://astronomy.swin.edu.au/cosmos/P/Period+Derivative
*/
gravitationalWavePeriodDerivative[mass1 is mass, mass2 is mass, ecc is dimensionless, period is time] :=
{
   fe = (1 + 73/24 ecc^2 + 37/96 ecc^4)/((1-ecc^2)^(7/2))
   dP = -192 pi / (5 c^5) G^(5/3) (period/(2 pi))^(-5/3) (mass1 mass2) /(mass1 + mass2)^(1/3) fe
   return dP
}

/** This is the simplified Einstein equation for the bending of light near a
    massive object.  M is the mass of the object and r0 is the radius of the
    closest approach to the center of the body.  It is only valid when r0 is
    much less than the distance to the emitting object or the viewing location.

    This is about 1.75 arcsec for starlight grazing the surface of the Sun.
*/
lightBendingEinstein[M is mass, r0 is length] :=
{
   rs = blackHoleRadius[M]
   return 2 rs / r0
}

/** Calculate the light-bending angle due to gravity.  This is valid when the
    source and destination are infinitely far from the gravitating object.  It
    is more accurate than the Einstein equation above.

  See:
  APPROXIMATION OF LIGHT-RAY DEFLECTION ANGLE AND GRAVITATIONAL LENSES IN THE
  SCHWARZSCHILD METRIC. II. LENSING MAGNIFICATION IN A BINARY SYSTEM
     Petri T. Mutka and Petri Mahonen, The Astrophysical Journal, 581:1328–1336, 2002 December 20
     https://iopscience.iop.org/article/10.1086/344438/pdf

   For a more detailed equation (that I haven't gotten to work) see:

   An accurate equation for the gravitational bending of light by a static
   massive object Open Access Oscar del Barco,
   Monthly Notices of the Royal Astronomical Society, Volume 535, Issue 3,
   December 2024, Pages 2504–2510,
   https://doi.org/10.1093/mnras/stae2277

   https://academic.oup.com/mnras/article/535/3/2504/7810584
*/
lightBendingDarwin[M is mass, r0 is length] :=
{
   rs = blackHoleRadius[M]
   R = r0/(G M / c^2)  // Is this right?  WTF are "gravitationalized units?"
   Q = sqrt[(R-2)(R+6)]
   k = sqrt[(Q-R+6)/(2 Q)]
   phi = arcsin[sqrt[(2+Q-R)/(6+Q-R)]]
   f1 = ellipticFDarwin[pi/2, k]
   f2 = ellipticFDarwin[phi, k]
   return 4 sqrt[R/Q] (f1-f2) - pi
}

/** This is the elliptical integral of the first type using the convention
    expected by lightBendingDarwin above.  There are lots of different
    conventions for this integral, some with m, some with m^2, some with the
    sin of a term, and this is to avoid ambiguity.

    THINK ABOUT: Make this library depend on NumericalIntegration.frink ?
*/
ellipticFDarwin[x,m] :=
{
   unsafeEval["use NumericalIntegration.frink"]
   f = {|t,m| 1 / sqrt[1 - m^2 sin[t]^2] }
   return NIntegrateData[f, 0, x, m]
}

/** Calculate the Unruh temperature of an accelerating object in a vacuum.
    This is important to understand black hole radiation / Hawking radiation:

    https://bigthink.com/starts-with-a-bang/hawking-radiation-really-work/

    Also see:
    https://en.wikipedia.org/wiki/Unruh_effect
*/
UnruhTemperature[a is acceleration] :=
{
   (hbar a) / (2 pi c k)
}

/** Calculate the Unruh acceleration corresponding to a temperature T in a
    vacuum.

    See:
    https://en.wikipedia.org/wiki/Unruh_effect
*/
UnruhAcceleration[T is temperature] :=
{
   (2 pi c k T) / hbar
}

/** Gravitational binding energy of a symmetric uniform spherical body.  The
    negative of this is the energy that must be added to a body to make it
    cease to be in a gravitationally bound state.  This is the energy required
    to lift all of its mass to an infinite distance at zero velocity.
    (Note that almost no bodies are of uniform density.   One estimate is that
     the Earth's actual gravitational binding energy is 1.11 times this large
     given a realistic estimate of densities.)

    See:
    https://en.wikipedia.org/wiki/Gravitational_binding_energy
*/
gravitationalBindingEnergy[m is mass, radius is length] :=
{
   -3 G m^2 / (5 radius)
}