膨胀球壳延缓引力势中的运动物体/Java程序
以下 Java 程序可用于数值计算不同线性质量密度下的史瓦西距离。它实现了有效质量的数值积分以及史瓦西距离的数值求解。
/*
Source file: SchwarzschildDistance.java
Program: Computation of the Schwarzschild Distance within the Hubble sphere
Autor: Markus Bautsch
Location: Berlin, Germany
Licence: public domain
Date: 18th July 2024
Version: 1.3
Programming language: Java
*/
/*
This Java-Programm computes Schwarzschild distances for any black shell masses.
The black shell is at the outer rim of the Hubble space.
*/
public class SchwarzschildDistance
{
// Class constants
final static double G = 6.67430e-11; // Gravitational constant in cube metres per kilogram and square second
final static double c = 2.99792458e8; // Speed of light in metres per second
final static double TropicalYear = 365.24219052 * 24 * 3600; // Tropical year in seconds
final static double LightYear = c * TropicalYear; // Light-year in metres
final static double R = 1.36e26; // Hubble length in metres
final static double MUniverse = 2.97e53; // Mass of the visible universe in kilograms
final static double lambdaMfirst = 3.36663e26; // start value for linear mass density of black shell
final static double lambdaMlast = 3.36670e26; // startop value for linear mass density of black shell
final static double deltaLambda = 0.00001e26; // step size for linear mass density of black shell
// This method computes and returns the distance s of a mass m to a mass point dM in the black shell
// distance is the closest distance of the mass m to to the mass point dM in the black shell in m
// alpha is the angle beweteen the mass m and the mass point dM in the black shell as seen from the centre of the universe in rad
private static double distanceToMassElementShell (double distance, double alpha)
{
double cosAlpha = java.lang.Math.cos (alpha);
double squareDistance = distance * distance;
double argument = 2*R*R - 2*R*distance + squareDistance - 2*R*(R-distance)*cosAlpha;
double s;
if (argument <= squareDistance) // because of possible rounding errors
{
s = distance;
}
else
{
s = java.lang.Math.sqrt (argument);
}
return s;
}
// This method computes and returns the angle of a half chord x in m as seen from the centre of the universe in rad
private static double alpha (double x)
{
double alpha = java.lang.Math.asin (x / R);
return alpha;
}
// This method computes and returns the integrand of the integral
// distance is the closest distance of the mass m to to the mass point dM in the black shell in m
// alpha is the angle beweteen the mass m and the mass point dM in the black shell as seen from the centre of the universe in rad
private static double integrand (double distance, double alpha)
{
double s = distanceToMassElementShell (distance, alpha);
double e = R * (1 - java.lang.Math.cos (alpha)); // the auxiliary sagitta e
double cosBeta = (distance - e) / s;
double integrand = cosBeta / s / s;
return integrand;
}
// This method computes and returns the effective mass from 0 to alphaR for a given lambdaM and a given distance
// lambdaM is the linear mass density of the black shell in kg/m
// distance is the closest distance of the mass m to to the mass point dM in the black shell in m
private static double integrateEffectiveMass (double lambdaM, double distance)
{
final long steps = 1000000; // number of steps for numerical integration
double squareDistance = distance * distance;
double xMax = java.lang.Math.sqrt (2*R*distance - squareDistance); // maximum half chord for a given maximum angle alphaR
double x = xMax; // x runs from xMax to 0 in steps of deltaX
double alpha = alpha (x); // alpha runs from alphaR to 0
double integrand = integrand (distance, alpha);
double integral = 0;
long step = 0;
while (step <= steps)
{
x = xMax * (1 - (double) step / steps);
double alphaNext = alpha (x);
double integrandNext = integrand (distance, alphaNext);
double deltaAlpha = alpha - alphaNext;
double averageIntegrand = (integrandNext + integrand) / 2;
double area = averageIntegrand * deltaAlpha; // numerical integration
integral = integral + area;
alpha = alphaNext;
integrand = integrandNext;
step++;
}
double effectiveMass = R * lambdaM * squareDistance * (2*integral);
return effectiveMass;
}
// This method computes and returns the Schwarzschild distance d_S for a given effective mass mEff
private static double schwarzschildDistance (double effectiveMass)
{
double schwarzschildDistance = 2 * G * effectiveMass / c / c;
return schwarzschildDistance;
}
// This method computes and returns the effective mass mEff for a given Schwarzschild distance d_S
private static double effectiveMass (double schwarzschildDistance)
{
double effectiveMass = schwarzschildDistance * c * c / 2 / G;
return effectiveMass;
}
// This method solves and returns the Schwarzschild distance d_S for a given lambdaM
// lambdaM is the linear mass density of the black shell in kg/m
private static double solve (double lambdaM)
{
final double limit = 5e-16; // for relative precision of determination of Schwarzschild distance
double lowerDistance = 1e20; // first low guess for Schwarzschild distance
double upperDistance = 1e25; // first high guess for Schwarzschild distance
double delta; // delta shall become smaller than limit
double distance;
double schwarzschildDistance;
double deltaDistance;
do
{
distance = (lowerDistance + upperDistance) / 2;
double effectiveMass = integrateEffectiveMass (lambdaM, distance);
schwarzschildDistance = schwarzschildDistance (effectiveMass);
deltaDistance = (distance - schwarzschildDistance) / schwarzschildDistance;
if (deltaDistance < 0) // Schwarzschild distance too large
{
lowerDistance = (upperDistance + lowerDistance) / 2;
}
else // Schwarzschild distance too small
{
upperDistance = (upperDistance + lowerDistance) / 2;
}
delta = (upperDistance - lowerDistance) / lowerDistance;
} while (delta > limit);
if (java.lang.Math.abs (deltaDistance) > 1e-6)
{
java.lang.System.out.println ("unstable numerical result:");
java.lang.System.out.println ("Computed distance = " + distance);
java.lang.System.out.println ("Schwarzschild distance = " + schwarzschildDistance);
}
return schwarzschildDistance;
}
// This method outputs the table header
private static void printHeader ()
{
java.lang.System.out.print ("lambda_M in kg/m;");
java.lang.System.out.print (" M_S in kg ;");
java.lang.System.out.print (" M_S/M ;");
java.lang.System.out.print (" M_eff in kg ;");
java.lang.System.out.print (" d_S im m ;");
java.lang.System.out.print (" d_S/R ;");
java.lang.System.out.print (" d_S in ly ;");
java.lang.System.out.print (" (R-d)/R ");
java.lang.System.out.println ();
}
// This method outputs a table result line
// lambdaM is the linear mass density of the black shell in kg/m
// distance is the closest distance of the mass m to to the mass point dM in the black shell in m
// mEff is the effective mass of a sphere in the black shell
private static void printResults (double lambdaM, double distance, double effectiveMass)
{
double mShell = lambdaM * R * 2 * java.lang.Math.PI; // mass of black shell
double mRatio = mShell / MUniverse;
double dRatio = distance / R;
double dLightYears = distance / LightYear;
double deltaRatio = (R - distance) / R;
java.lang.System.out.printf ("%15.7e", lambdaM);
java.lang.System.out.print (" ;");
java.lang.System.out.printf ("%15.6e", mShell);
java.lang.System.out.print (" ;");
java.lang.System.out.printf ("%15.6f", mRatio);
java.lang.System.out.print (" ;");
java.lang.System.out.printf ("%15.2e", effectiveMass);
java.lang.System.out.print (" ;");
java.lang.System.out.printf ("%15.2e", distance);
java.lang.System.out.print (" ;");
java.lang.System.out.printf ("%15.2e", dRatio);
java.lang.System.out.print (" ;");
java.lang.System.out.printf ("%15.0f", dLightYears);
java.lang.System.out.print (" ;");
java.lang.System.out.printf ("%15.9f", + deltaRatio);
java.lang.System.out.println ();
}
// This method has a loop for computing a sequence of effective masses and Schwarzschild distances for different lambdaM values
// lambdaM is iterated from lambdaMfirst to lambdaMlast in steps of deltaLambda
private static void computeScharzschildRadii (double lambdaMfirst, double lambdaMlast, double deltaLambda)
{
double lambdaM = lambdaMfirst; // linear mass density of black shell at Hubble radius in kg/m
lambdaMlast = lambdaMlast * 1.00000001; // for the reason of numerical rounding errors
printHeader ();
do
{
double schwarzschildDistance = solve (lambdaM);
double effectiveMass = effectiveMass (schwarzschildDistance);
printResults (lambdaM, schwarzschildDistance, effectiveMass);
lambdaM = lambdaM + deltaLambda;
} while (lambdaM <= lambdaMlast);
}
// Main program
public static void main (java.lang.String [] arguments)
{
computeScharzschildRadii (lambdaMfirst, lambdaMlast, deltaLambda);
}
}