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();
|
||
|
||
|
||
}
|
||
} |