330 lines
10 KiB
C#
330 lines
10 KiB
C#
|
using System;
|
|||
|
|
using System.Collections;
|
|||
|
|
using System.Collections.Generic;
|
|||
|
|
using System.Security.Cryptography;
|
|||
|
|
using JetBrains.Annotations;
|
|||
|
|
using TMPro;
|
|||
|
|
using UnityEngine;
|
|||
|
|
using UnityEngine.Rendering;
|
|||
|
|
|
|||
|
|
public class BlackHoleRenderer : MonoBehaviour {
|
|||
|
|
void Start() {
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
struct Ray
|
|||
|
|
{
|
|||
|
|
public SphericalDir Dir;
|
|||
|
|
public SphericalCoord P; //Canonical momenta
|
|||
|
|
public double b;
|
|||
|
|
public double q;
|
|||
|
|
|
|||
|
|
public Double3 DirectionOfMotionInCameraCortesianSpace(Double3 N, double beta) { //F
|
|||
|
|
double betaSquare = beta * beta;
|
|||
|
|
double nfy = (-N.Y + beta) / (1.0 - (beta * N.Y));
|
|||
|
|
double nfx = (-Math.Sqrt(1.0 - betaSquare) * N.X) / (1.0 - beta * N.Y);
|
|||
|
|
double nfz = (-Math.Sqrt(1.0 - betaSquare) * N.Z) / (1.0 - beta * N.Y);
|
|||
|
|
return new Double3 { X = nfx, Y = nfy, Z = nfz };
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public SphericalCoord FSpherical(SphericalCoord B, Double3 nF) {
|
|||
|
|
double kappa = Math.Sqrt(B.R * B.R + B.Dir.Phi * B.Dir.Phi);
|
|||
|
|
|
|||
|
|
double nFr = ((B.Dir.Phi / kappa) * nF.X) + (B.R * nF.Y) + ((B.R * B.Dir.Theta) / kappa) * nF.Z;
|
|||
|
|
double nFTheta = B.Dir.Theta * nF.Y - kappa * nF.Z;
|
|||
|
|
double nFPhi = (-B.R / kappa) * nF.X + B.Dir.Phi * nF.Y + ((B.Dir.Phi * B.Dir.Theta) / kappa) * nF.Z;
|
|||
|
|
|
|||
|
|
SphericalCoord result;
|
|||
|
|
result.R = nFr;
|
|||
|
|
result.Dir.Theta = nFTheta;
|
|||
|
|
result.Dir.Phi = nFPhi;
|
|||
|
|
return result;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public void CanonicalMomenta(FidoCamera camera, KerrMetric metric) { //p
|
|||
|
|
var N = Dir.ToCartesian();
|
|||
|
|
var Fdir = DirectionOfMotionInCameraCortesianSpace(N, camera.Beta);
|
|||
|
|
var F = FSpherical(camera.Speed, Fdir);
|
|||
|
|
|
|||
|
|
double Ef = 1.0 / (metric.Alpha + metric.Omega * metric.OmegaUpLine * F.Dir.Phi);
|
|||
|
|
|
|||
|
|
double pr = Ef * (metric.Rho / Math.Sqrt(metric.Delta)) * F.R;
|
|||
|
|
double pTheta = Ef * metric.Rho * F.Dir.Theta;
|
|||
|
|
double pPhi = Ef * metric.OmegaUpLine * F.Dir.Phi;
|
|||
|
|
|
|||
|
|
SphericalCoord result;
|
|||
|
|
result.R = pr;
|
|||
|
|
result.Dir.Phi = pPhi;
|
|||
|
|
result.Dir.Theta = pTheta;
|
|||
|
|
|
|||
|
|
P = result;
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public void Update(FidoCamera camera, KerrMetric metric)
|
|||
|
|
{
|
|||
|
|
CanonicalMomenta(camera, metric);
|
|||
|
|
|
|||
|
|
b = P.Dir.Phi;
|
|||
|
|
double cosTheta = Math.Cos(Dir.Theta);
|
|||
|
|
double sinTheta = Math.Sin(Dir.Theta);
|
|||
|
|
q = P.Dir.Theta * P.Dir.Theta +
|
|||
|
|
cosTheta * cosTheta * (((b * b) / (sinTheta * sinTheta)) - metric.A * metric.A);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double betaZero(double r_zero, BlackHole hole)
|
|||
|
|
{
|
|||
|
|
double r_zero_square = r_zero * r_zero;
|
|||
|
|
double r_zero_cube = r_zero_square * r_zero;
|
|||
|
|
|
|||
|
|
double a_square = hole.A * hole.A;
|
|||
|
|
return -((r_zero_cube - 3 * r_zero_square + a_square * r_zero + a_square) / (hole.A * (r_zero - 1.0)));
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double qZero(double r, BlackHole hole) {
|
|||
|
|
double rSquare = r * r;
|
|||
|
|
double rCube = rSquare * r;
|
|||
|
|
double aSquare = hole.A * hole.A;
|
|||
|
|
double rm1 = r - 1;
|
|||
|
|
double rmSquare = rm1 * rm1;
|
|||
|
|
return -(rCube * (rCube - 6 * rSquare + 9 * r - 4 * aSquare)) / (aSquare * rmSquare);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double rZero_0(double b, double a) { //my
|
|||
|
|
var aSquare = a * a;
|
|||
|
|
var subExpr0 = 3 * aSquare + 3 * a * b - 9;
|
|||
|
|
var subExpr2 = (54 - 54 * aSquare);
|
|||
|
|
|
|||
|
|
var cubeRootOf2 = Math.Pow(2.0, 1.0 / 3.0);
|
|||
|
|
var subExpr1 = Math.Pow(Math.Sqrt(Math.Abs(4 * subExpr0 * subExpr0 * subExpr0 + subExpr2 * subExpr2)) + subExpr2, 1.0 / 3.0);
|
|||
|
|
|
|||
|
|
return -(cubeRootOf2 * subExpr0) / (3.0 * subExpr1) + (subExpr1 / (3.0 * cubeRootOf2)) + 1;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double rZero_1(double b, double a) { //my
|
|||
|
|
var aSquare = a * a;
|
|||
|
|
var subExpr0 = 3 * aSquare + 3 * a * b - 9;
|
|||
|
|
var subExpr2 = (54 - 54 * aSquare);
|
|||
|
|
|
|||
|
|
var cubeRootOf2 = Math.Pow(2.0, 1.0 / 3.0);
|
|||
|
|
|
|||
|
|
var subExpr1 = Math.Pow(Math.Sqrt(Math.Abs(4 * subExpr0 * subExpr0 * subExpr0 + subExpr2 * subExpr2)) + subExpr2, 1.0 / 3.0);
|
|||
|
|
return ((1 + Math.Sqrt(3.0)) * subExpr0) / (3.0 * Math.Pow(2.0, 2.0/3.0) * subExpr1) - (( (1 - Math.Sqrt(3)) * subExpr1) / (6.0 * cubeRootOf2)) + 1;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double rZero_2(double b, double a) { //my
|
|||
|
|
var aSquare = a * a;
|
|||
|
|
var subExpr0 = 3 * aSquare + 3 * a * b - 9;
|
|||
|
|
var subExpr2 = (54 - 54 * aSquare);
|
|||
|
|
|
|||
|
|
var cubeRootOf2 = Math.Pow(2.0, 1.0 / 3.0);
|
|||
|
|
|
|||
|
|
var subExpr1 = Math.Pow(Math.Sqrt(Math.Abs(4 * subExpr0 * subExpr0 * subExpr0 + subExpr2 * subExpr2)) + subExpr2, 1.0 / 3.0);
|
|||
|
|
return ((1 - Math.Sqrt(3.0)) * subExpr0) / (3.0 * Math.Pow(2.0, 2.0 / 3.0) * subExpr1) - (((1 + Math.Sqrt(3)) * subExpr1) / (6.0 * cubeRootOf2)) + 1;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
[Serializable]
|
|||
|
|
public class FidoCamera
|
|||
|
|
{
|
|||
|
|
public SphericalCoord Position;
|
|||
|
|
public SphericalCoord Speed;
|
|||
|
|
public double Beta;
|
|||
|
|
|
|||
|
|
public void UpdateMetricAndSelf(KerrMetric metric, BlackHole hole)
|
|||
|
|
{
|
|||
|
|
metric.Update(hole.A, Position.R, Position.Dir.Theta);
|
|||
|
|
double bigOmega = 1.0 / (metric.A + Math.Pow(Position.R, 3.0 / 2.0));
|
|||
|
|
Beta = (metric.OmegaUpLine / metric.Alpha) * (bigOmega - metric.Omega);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
[Serializable]
|
|||
|
|
public class BlackHole
|
|||
|
|
{
|
|||
|
|
public double A;
|
|||
|
|
public double r1() {
|
|||
|
|
return 2 * (1 + Math.Cos((2.0 / 3.0) * Math.Acos(-A)));
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public double r2() {
|
|||
|
|
return 2 * (1 + Math.Cos((2.0 / 3.0) * Math.Acos(A)));
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
[Serializable]
|
|||
|
|
public class KerrMetric {
|
|||
|
|
public double Delta;
|
|||
|
|
public double Rho;
|
|||
|
|
public double OmegaUpLine;
|
|||
|
|
public double Omega;
|
|||
|
|
public double Sigma;
|
|||
|
|
public double Alpha;
|
|||
|
|
public double A;
|
|||
|
|
|
|||
|
|
public void Update(double a, double r, double theta) {
|
|||
|
|
double cosTheta = Math.Cos(theta);
|
|||
|
|
double sinTheta = Math.Sin(theta);
|
|||
|
|
|
|||
|
|
double cosThetaSquare = cosTheta * cosTheta;
|
|||
|
|
double sinThetaSquare = sinTheta * sinTheta;
|
|||
|
|
|
|||
|
|
double aSquare = a * a;
|
|||
|
|
double rSquare = r * r;
|
|||
|
|
|
|||
|
|
Rho = Math.Sqrt(rSquare + (aSquare * cosThetaSquare));
|
|||
|
|
Delta = rSquare - (2 * r) + aSquare;
|
|||
|
|
Sigma = Math.Sqrt((rSquare + aSquare) * (rSquare + aSquare) - aSquare * Delta * sinThetaSquare);
|
|||
|
|
Alpha = (Rho * Math.Sqrt(Delta)) / Sigma;
|
|||
|
|
Omega = (2 * a * r) / (Sigma * Sigma);
|
|||
|
|
OmegaUpLine = (Sigma * sinTheta) / Rho;
|
|||
|
|
|
|||
|
|
A = a;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public FidoCamera Camera = new FidoCamera();
|
|||
|
|
public KerrMetric Metric = new KerrMetric();
|
|||
|
|
public BlackHole Hole = new BlackHole();
|
|||
|
|
|
|||
|
|
public void RecalcState()
|
|||
|
|
{
|
|||
|
|
Camera.UpdateMetricAndSelf(Metric, Hole);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double drdt(double delta, double rho, double Pr)
|
|||
|
|
{
|
|||
|
|
return delta / (rho * rho) * Pr;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
SphericalDir[] rays = new SphericalDir[256];
|
|||
|
|
|
|||
|
|
class BlackHoleRayIntegrator : Integrator
|
|||
|
|
{
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void DrawCameraGizmo() {
|
|||
|
|
var pos = Camera.Position.ToCartesian().ToVector3();
|
|||
|
|
|
|||
|
|
UpdateRays();
|
|||
|
|
|
|||
|
|
foreach (var rayDir in rays)
|
|||
|
|
{
|
|||
|
|
Ray ray = new Ray();
|
|||
|
|
ray.Dir = rayDir;
|
|||
|
|
ray.Update(Camera, Metric);
|
|||
|
|
|
|||
|
|
|
|||
|
|
var v0 = rZero_0(0.666, 0.9);
|
|||
|
|
var v1 = rZero_1(0.666, 0.9);
|
|||
|
|
var v2 = rZero_2(0.666, 0.9);
|
|||
|
|
// var b1 = betaZero(Hole.r2());
|
|||
|
|
// var b2 = betaZero(Hole.r1());
|
|||
|
|
|
|||
|
|
//-0.1184053483681316836 + 0.×10^-19 i
|
|||
|
|
//0.7514421233051754296 + 0.×10^-19 i
|
|||
|
|
//2.3669632250629562540 + 0.×10^-20 i
|
|||
|
|
var N = ray.Dir.ToCartesian();
|
|||
|
|
var Fdir = ray.DirectionOfMotionInCameraCortesianSpace(N, Camera.Beta);
|
|||
|
|
var F = ray.FSpherical(Camera.Speed, Fdir);
|
|||
|
|
|
|||
|
|
Gizmos.DrawLine(pos, pos + ray.P.ToCartesian().ToVector3());
|
|||
|
|
|
|||
|
|
/*
|
|||
|
|
var r0 = rZero2(b);
|
|||
|
|
var r0_ = rZero(b);
|
|||
|
|
//r0 = 2.1899962982234369;
|
|||
|
|
var q0 = qZero2(r0, b);
|
|||
|
|
|
|||
|
|
|
|||
|
|
Color rayColor = Color.white;
|
|||
|
|
if (((b1 < b) && (b < b2)) && (q < q0))
|
|||
|
|
{
|
|||
|
|
//there are no radial turning points for that {b, q}
|
|||
|
|
if (rayCanonicalMomenta_.R > 0)
|
|||
|
|
{
|
|||
|
|
//horizon
|
|||
|
|
rayColor = Color.black;
|
|||
|
|
horizon++;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
celestial++;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
var a = a_spinAngularMomentumPerMass;
|
|||
|
|
|
|||
|
|
double P = Math.Sqrt(Metric.Delta * ((b - a) * (b - a) + q));
|
|||
|
|
double rUp0 = P + a * b - a * a;
|
|||
|
|
double rUp1 = -P + a * b - a * a;
|
|||
|
|
|
|||
|
|
double rUp = rUp0 > rUp1 ? rUp0 : rUp1;
|
|||
|
|
|
|||
|
|
if (cameraPosition.R > rUp)
|
|||
|
|
{
|
|||
|
|
celestial++;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
horizon++;
|
|||
|
|
//horizon
|
|||
|
|
rayColor = Color.black;
|
|||
|
|
}
|
|||
|
|
}*/
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public void OnDrawGizmos() {
|
|||
|
|
|
|||
|
|
|
|||
|
|
// double accretionDiskInnerRadius = 9.26 * M_blackHoleMass;
|
|||
|
|
// double accretionDiskOuterRadius = 18.7 * M_blackHoleMass;
|
|||
|
|
|
|||
|
|
Gizmos.DrawLine(Vector3.zero, Camera.Position.ToCartesian().ToVector3());
|
|||
|
|
Gizmos.DrawWireSphere(Vector3.zero, 1);
|
|||
|
|
|
|||
|
|
DrawCameraGizmo();
|
|||
|
|
// double OmegaBig = CameraGeodesicAngularVelocity(Rc);
|
|||
|
|
//DrawCameraGizmo(cameraCoord, cameraDirectionOfMotion);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void UpdateRays()
|
|||
|
|
{
|
|||
|
|
int imageHeight = 16;
|
|||
|
|
int imageWidth = 16;
|
|||
|
|
|
|||
|
|
int totalRays = imageHeight * imageWidth;
|
|||
|
|
|
|||
|
|
if (rays == null || rays.Length != totalRays)
|
|||
|
|
{
|
|||
|
|
rays = new SphericalDir[totalRays];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
double verticalFov = Math.PI / 2;
|
|||
|
|
double horizontalFov = ((double)imageWidth / imageHeight) * verticalFov;
|
|||
|
|
|
|||
|
|
for (int y = 0; y < imageHeight; y++)
|
|||
|
|
{
|
|||
|
|
double rayTheta = (Math.PI / 2) - verticalFov / 2 + (verticalFov / (imageHeight - 1) * y);
|
|||
|
|
for (int x = 0; x < imageWidth; x++)
|
|||
|
|
{
|
|||
|
|
double rayPhi = Math.PI + horizontalFov / 2 - (horizontalFov / (imageWidth - 1) * x);
|
|||
|
|
|
|||
|
|
rays[y * imageWidth + x].Phi = rayPhi;
|
|||
|
|
rays[y * imageWidth + x].Theta = rayTheta;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public void Trace()
|
|||
|
|
{
|
|||
|
|
RecalcState();
|
|||
|
|
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
}
|