645 lines
43 KiB
C++
645 lines
43 KiB
C++
//This function contains the equations of motion requiered for the integration (u0Dot - u3Dot)
|
||
function matrix3 getDerivatives(matrix3 FunctionIncoming; float a)
|
||
{
|
||
float t,r,theta,phi,u0,u1,u2,u3,Null;
|
||
assign(t,r,theta,phi,u0,u1,u2,u3,Null,FunctionIncoming);
|
||
float u0Dot = -2*u0*u1*(-1.0*1.0*a*r*(4*1.0*a*pow(r, 2)*pow(sin(theta), 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) - 2*1.0*a*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4)) + 0.5*(-4*1.0*pow(r, 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 2*1.0/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*(-2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - pow(a, 4)*pow(cos(theta), 2) - pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) - pow(a, 2)*pow(r, 2) - pow(r, 4))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4))) - 2*u0*u2*(2.0*1.0*pow(a, 2)*r*(-2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - pow(a, 4)*pow(cos(theta), 2) - pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) - pow(a, 2)*pow(r, 2) - pow(r, 4))*sin(theta)*cos(theta)/(pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2)*(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4))) - 1.0*1.0*a*r*(-4*1.0*pow(a, 3)*r*pow(sin(theta), 3)*cos(theta)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) - 4*1.0*a*r*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4))) - 2*u1*u3*(-1.0*1.0*a*r*(-4*1.0*pow(a, 2)*pow(r, 2)*pow(sin(theta), 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 2*1.0*pow(a, 2)*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) + 2*r)*pow(sin(theta), 2)/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4)) + 0.5*(4*1.0*a*pow(r, 2)*pow(sin(theta), 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) - 2*1.0*a*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*(-2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - pow(a, 4)*pow(cos(theta), 2) - pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) - pow(a, 2)*pow(r, 2) - pow(r, 4))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4))) - 2*u2*u3*(-1.0*1.0*a*r*((4*1.0*pow(a, 4)*r*pow(sin(theta), 3)*cos(theta)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 4*1.0*pow(a, 2)*r*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*pow(sin(theta), 2) + 2*(2*1.0*pow(a, 2)*r*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) + pow(a, 2) + pow(r, 2))*sin(theta)*cos(theta))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4)) + 0.5*(-4*1.0*pow(a, 3)*r*pow(sin(theta), 3)*cos(theta)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) - 4*1.0*a*r*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*(-2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - pow(a, 4)*pow(cos(theta), 2) - pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) - pow(a, 2)*pow(r, 2) - pow(r, 4))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4)));
|
||
float u1Dot = 2.0*pow(a, 2)*u1*u2*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) + 1.0*r*pow(u2, 2)*(-2*1.0*r + pow(a, 2) + pow(r, 2))/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) - 0.5*pow(u0, 2)*(4*1.0*pow(r, 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) - 2*1.0/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*(-2*1.0*r + pow(a, 2) + pow(r, 2))/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) - 1.0*u0*u3*(-4*1.0*a*pow(r, 2)*pow(sin(theta), 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 2*1.0*a*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*(-2*1.0*r + pow(a, 2) + pow(r, 2))/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) - 0.5*pow(u1, 2)*(2*r/(-2*1.0*r + pow(a, 2) + pow(r, 2)) + (2*1.0 - 2*r)*(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2))/pow(-2*1.0*r + pow(a, 2) + pow(r, 2), 2))*(-2*1.0*r + pow(a, 2) + pow(r, 2))/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) + 0.5*pow(u3, 2)*(-2*1.0*r + pow(a, 2) + pow(r, 2))*(-4*1.0*pow(a, 2)*pow(r, 2)*pow(sin(theta), 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 2*1.0*pow(a, 2)*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) + 2*r)*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2));
|
||
float u2Dot = 2.0*1.0*pow(a, 2)*r*pow(u0, 2)*sin(theta)*cos(theta)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 3) - 1.0*pow(a, 2)*pow(u1, 2)*sin(theta)*cos(theta)/((pow(a, 2)*pow(cos(theta), 2) + pow(r, 2))*(-2*1.0*r + pow(a, 2) + pow(r, 2))) + 1.0*pow(a, 2)*pow(u2, 2)*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) - 2.0*r*u1*u2/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) - 1.0*u0*u3*(4*1.0*pow(a, 3)*r*pow(sin(theta), 3)*cos(theta)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 4*1.0*a*r*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) - 0.5*pow(u3, 2)*(-(4*1.0*pow(a, 4)*r*pow(sin(theta), 3)*cos(theta)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 4*1.0*pow(a, 2)*r*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*pow(sin(theta), 2) - 2*(2*1.0*pow(a, 2)*r*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) + pow(a, 2) + pow(r, 2))*sin(theta)*cos(theta))/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2));
|
||
float u3Dot = -2*u0*u1*(-1.0*1.0*a*r*(-4*1.0*pow(r, 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 2*1.0/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4)) + 0.5*(4*1.0*a*pow(r, 2)*pow(sin(theta), 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) - 2*1.0*a*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*(-2*1.0*r + pow(a, 2)*pow(cos(theta), 2) + pow(r, 2))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 4) - 2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(r, 3)*pow(sin(theta), 2) + pow(a, 4)*pow(sin(theta), 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(sin(theta), 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(sin(theta), 2) + pow(r, 4)*pow(sin(theta), 2))) - 2*u0*u2*(-4.0*pow(1.0, 2)*pow(a, 3)*pow(r, 2)*sin(theta)*cos(theta)/(pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2)*(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4))) + 0.5*(-4*1.0*pow(a, 3)*r*pow(sin(theta), 3)*cos(theta)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) - 4*1.0*a*r*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*(-2*1.0*r + pow(a, 2)*pow(cos(theta), 2) + pow(r, 2))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 4) - 2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(r, 3)*pow(sin(theta), 2) + pow(a, 4)*pow(sin(theta), 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(sin(theta), 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(sin(theta), 2) + pow(r, 4)*pow(sin(theta), 2))) - 2*u1*u3*(-1.0*1.0*a*r*(4*1.0*a*pow(r, 2)*pow(sin(theta), 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) - 2*1.0*a*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4)) + 0.5*(-2*1.0*r + pow(a, 2)*pow(cos(theta), 2) + pow(r, 2))*(-4*1.0*pow(a, 2)*pow(r, 2)*pow(sin(theta), 2)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 2*1.0*pow(a, 2)*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) + 2*r)*pow(sin(theta), 2)/(2*1.0*pow(a, 2)*r*pow(sin(theta), 4) - 2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(r, 3)*pow(sin(theta), 2) + pow(a, 4)*pow(sin(theta), 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(sin(theta), 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(sin(theta), 2) + pow(r, 4)*pow(sin(theta), 2))) - 2*u2*u3*(-1.0*1.0*a*r*(-4*1.0*pow(a, 3)*r*pow(sin(theta), 3)*cos(theta)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) - 4*1.0*a*r*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(a, 2)*r - 2*1.0*pow(r, 3) + pow(a, 4)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2) + pow(r, 4)) + 0.5*((4*1.0*pow(a, 4)*r*pow(sin(theta), 3)*cos(theta)/pow(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2), 2) + 4*1.0*pow(a, 2)*r*sin(theta)*cos(theta)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)))*pow(sin(theta), 2) + 2*(2*1.0*pow(a, 2)*r*pow(sin(theta), 2)/(pow(a, 2)*pow(cos(theta), 2) + pow(r, 2)) + pow(a, 2) + pow(r, 2))*sin(theta)*cos(theta))*(-2*1.0*r + pow(a, 2)*pow(cos(theta), 2) + pow(r, 2))/(2*1.0*pow(a, 2)*r*pow(sin(theta), 4) - 2*1.0*pow(a, 2)*r*pow(sin(theta), 2) - 2*1.0*pow(r, 3)*pow(sin(theta), 2) + pow(a, 4)*pow(sin(theta), 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(sin(theta), 2)*pow(cos(theta), 2) + pow(a, 2)*pow(r, 2)*pow(sin(theta), 2) + pow(r, 4)*pow(sin(theta), 2)));
|
||
return set(u0,u1,u2,u3,u0Dot,u1Dot,u2Dot,u3Dot,0);
|
||
}
|
||
|
||
//Runge–Kutta–Fehlberg Scheme
|
||
function matrix3 getRKF45step(float it; float ir; float itheta; float iphi; float iu0; float iu1; float iu2; float iu3; float hU; float a; float Medium; int Bound)
|
||
{
|
||
float absError = 1e10, tol = 1e-5, newhV = 0, newhC = 0, newh = 0, e1,e2,e3,e4,e5,e6,e7,e8,e9,deltaStep,ts = 0.1;
|
||
matrix3 FunctionIncoming = set(it, ir, itheta, iphi, iu0, iu1, iu2, iu3, 0), newy = 0, error = 0, k1 = 0, k2 = 0, k3 = 0, k4 = 0, k5 = 0, k6 = 0;
|
||
|
||
while(absError > tol)
|
||
{
|
||
k1 = getDerivatives(FunctionIncoming,a);
|
||
k2 = getDerivatives(FunctionIncoming + (1.0/4.0)*k1*hU,a);
|
||
k3 = getDerivatives(FunctionIncoming + (3.0/32.0)*k1*hU + (9.0/32.0)*k2*hU,a);
|
||
k4 = getDerivatives(FunctionIncoming + (1932.0/2197.0)*k1*hU - (7200.0/2197.0)*k2*hU + (7296.0/2197.0)*k3*hU,a);
|
||
k5 = getDerivatives(FunctionIncoming + (439.0/216.0)*k1*hU - 8*k2*hU + (3680.0/513.0)*k3*hU - (845.0/4104.0)*k4*hU,a);
|
||
k6 = getDerivatives(FunctionIncoming - (8.0/27.0)*k1*hU + 2*k2*hU - (3544.0/2565.0)*k3*hU + (1859.0/4104.0)*k4*hU - (11.0/40.0)*k5*hU,a);
|
||
|
||
newy = FunctionIncoming + (16.0/135.0)*k1*hU + (6656.0/12825.0)*k3*hU + (28561.0/56430.0)*k4*hU - (9.0/50.0)*k5*hU + (2.0/55.0)*k6*hU;
|
||
error = (-1.0/360.0)*k1*hU + (128.0/4275.0)*k3*hU + (2197.0/75240.0)*k4*hU - (1.0/50.0)*k5*hU - (2.0/55.0)*k6*hU;
|
||
assign(e1,e2,e3,e4,e5,e6,e7,e8,e9,error);
|
||
absError = sqrt( (e1*e1) + (e2*e2) + (e3*e3) + (e4*e4) + (e5*e5) + (e6*e6) + (e7*e7) + (e8*e8) + (e9*e9));
|
||
|
||
if(Bound != 0)
|
||
{
|
||
ts = lerp(0.9,0.05,Medium);
|
||
}
|
||
|
||
newhV = ts*hU*pow((tol/absError),(1.0/5.0));
|
||
newhC = 0.9*hU*pow((tol/absError),(1.0/5.0));
|
||
newh = min(newhV,newhC);
|
||
|
||
if(absError > tol)
|
||
{
|
||
hU = newh;
|
||
}
|
||
}
|
||
|
||
float tDeriv, rDeriv, thetaDeriv, phiDeriv, u0Deriv, u1Deriv, u2Deriv, u3Deriv, a9;
|
||
assign(tDeriv, rDeriv, thetaDeriv, phiDeriv, u0Deriv, u1Deriv, u2Deriv, u3Deriv, a9, newy);
|
||
return set(tDeriv, rDeriv, thetaDeriv, phiDeriv, u0Deriv, u1Deriv, u2Deriv, u3Deriv, newh);
|
||
}
|
||
|
||
//Table
|
||
function vector getMatchingFunctions(float Wavelength)
|
||
{
|
||
vector my_Array[] = { {0.0014, 0.0000, 0.0065},{0.0022, 0.0001, 0.0105},{0.0042, 0.0001, 0.0201},{0.0076, 0.0002, 0.0362},{0.0143, 0.0004, 0.0679},{0.0232, 0.0006, 0.1102},{0.0435, 0.0012, 0.2074},{0.0776, 0.0022, 0.3713},{0.1344, 0.0040, 0.6456},{0.2148, 0.0073, 1.0391},{0.2839, 0.0116, 1.3856},{0.3285, 0.0168, 1.6230},{0.3483, 0.0230, 1.7471},{0.3481, 0.0298, 1.7826},{0.3362, 0.0380, 1.7721},{0.3187, 0.0480, 1.7441},{0.2908, 0.0600, 1.6692},{0.2511, 0.0739, 1.5281},{0.1954, 0.0910, 1.2876},{0.1421, 0.1126, 1.0419},{0.0956, 0.1390, 0.8130},{0.0580, 0.1693, 0.6162},{0.0320, 0.2080, 0.4652},{0.0147, 0.2586, 0.3533},{0.0049, 0.3230, 0.2720},{0.0024, 0.4073, 0.2123},{0.0093, 0.5030, 0.1582},{0.0291, 0.6082, 0.1117},{0.0633, 0.7100, 0.0782},{0.1096, 0.7932, 0.0573},{0.1655, 0.8620, 0.0422},{0.2257, 0.9149, 0.0298},{0.2904, 0.9540, 0.0203},{0.3597, 0.9803, 0.0134},{0.4334, 0.9950, 0.0087},{0.5121, 1.0000, 0.0057},{0.5945, 0.9950, 0.0039},{0.6784, 0.9786, 0.0027},{0.7621, 0.9520, 0.0021},{0.8425, 0.9154, 0.0018},{0.9163, 0.8700, 0.0017},{0.9786, 0.8163, 0.0014},{1.0263, 0.7570, 0.0011},{1.0567, 0.6949, 0.0010},{1.0622, 0.6310, 0.0008},{1.0456, 0.5668, 0.0006},{1.0026, 0.5030, 0.0003},{0.9384, 0.4412, 0.0002},{0.8544, 0.3810, 0.0002},{0.7514, 0.3210, 0.0001},{0.6424, 0.2650, 0.0000},{0.5419, 0.2170, 0.0000},{0.4479, 0.1750, 0.0000},{0.3608, 0.1382, 0.0000},{0.2835, 0.1070, 0.0000},{0.2187, 0.0816, 0.0000},{0.1649, 0.0610, 0.0000},{0.1212, 0.0446, 0.0000},{0.0874, 0.0320, 0.0000},{0.0636, 0.0232, 0.0000},{0.0468, 0.0170, 0.0000},{0.0329, 0.0119, 0.0000},{0.0227, 0.0082, 0.0000},{0.0158, 0.0057, 0.0000},{0.0114, 0.0041, 0.0000},{0.0081, 0.0029, 0.0000},{0.0058, 0.0021, 0.0000},{0.0041, 0.0015, 0.0000},{0.0029, 0.0010, 0.0000},{0.0020, 0.0007, 0.0000},{0.0014, 0.0005, 0.0000},{0.0010, 0.0004, 0.0000},{0.0007, 0.0002, 0.0000},{0.0005, 0.0002, 0.0000},{0.0003, 0.0001, 0.0000},{0.0002, 0.0001, 0.0000},{0.0002, 0.0001, 0.0000},{0.0001, 0.0000, 0.0000},{0.0001, 0.0000, 0.0000},{0.0001, 0.0000, 0.0000},{0.0000, 0.0000, 0.0000} };
|
||
int index = int(Wavelength - 380) / 5;
|
||
return my_Array[index];
|
||
}
|
||
|
||
//Blackbody Function
|
||
function vector4 getBlackbodyColor(float T; float Redshiftfactor)
|
||
{
|
||
float Intensity = 0, lambda = 0, cX = 0, cY = 0, cZ = 0, Wavelength = 380, LightSpeed = 2.99792458e8, kR = 1.3806504e-23, hR = 6.62607015e-34, Temperature = T / Redshiftfactor;
|
||
vector myVector = 0;
|
||
|
||
for(int i = 0; i < 80; i++)
|
||
{
|
||
Wavelength += 5;
|
||
lambda = Wavelength * 1e-9;
|
||
Intensity = ((2*hR*pow(LightSpeed,2)) / pow(lambda,5)) / (exp((hR*LightSpeed) / (lambda*kR*Temperature)) -1);
|
||
myVector = getMatchingFunctions(Wavelength);
|
||
cX += Intensity*myVector.x*5;
|
||
cY += Intensity*myVector.y*5;
|
||
cZ += Intensity*myVector.z*5;
|
||
}
|
||
|
||
float x = cX/(cX + cY + cZ), y = cY/(cX + cY + cZ), z = 1 - x - y;
|
||
float Y = 1.0, X = x/y, Z = z/y;
|
||
float R = 3.24096994*X - 1.53738318*Y - 0.49861076*Z;
|
||
float G = - 0.96924364*X + 1.8759675*Y + 0.04155506*Z;
|
||
float B = 0.05563008*X - 0.20397696*Y + 1.05697151*Z;
|
||
float max = max(R,G,B);
|
||
|
||
cX = max(cX,0);
|
||
cY = max(cY,0);
|
||
cZ = max(cZ,0);
|
||
|
||
if(max > 1)
|
||
{
|
||
R = R / max;
|
||
G = G / max;
|
||
B = B / max;
|
||
}
|
||
|
||
if(R < 0)
|
||
{
|
||
R = 0;
|
||
}
|
||
if(G < 0)
|
||
{
|
||
G = 0;
|
||
}
|
||
if(B < 0)
|
||
{
|
||
B = 0;
|
||
}
|
||
|
||
return set(R,G,B,cY);
|
||
}
|
||
|
||
//u0 Camera
|
||
function float getU0(float u1; float u2; float u3; float r; float theta; float a; float Delta; float Sigma; float My)
|
||
{
|
||
float gtPhi = -((2*a*r) / Sigma) * pow(sin(theta),2);
|
||
float gtt = -(1-((2*r) / Sigma));
|
||
float grr = Sigma / Delta;
|
||
float gThetaTheta = Sigma;
|
||
float Lambda = pow(pow(r,2) + pow(a,2),2) - pow(a,2)*Delta*pow(sin(theta),2);
|
||
float gPhiPhi = (Lambda/Sigma) * pow(sin(theta),2);
|
||
return (-gtPhi / gtt)*u3 + sqrt( pow(gtPhi/gtt,2)*pow(u3,2) - (1.0 / gtt) * ( grr*pow(u1,2) + gThetaTheta*pow(u2,2) + gPhiPhi*pow(u3,2) - My ));
|
||
}
|
||
|
||
//Jet Velocity Distribution
|
||
function vector2 getJetDisplacement(vector cameraRThetaphi; float Time; float a; float r; float theta; float t; float Delta; float Sigma; float circularVelocity; float verticalVelocity)
|
||
{
|
||
float u1 = sqrt(pow(r, 2) + pow(a, 2))/sqrt(pow(r, 2) + pow(a*cos(theta), 2))*cos(theta)*verticalVelocity;
|
||
float u2 = -r*sin(theta)/sqrt(pow(r, 2) + pow(a*cos(theta), 2))*verticalVelocity;
|
||
float u3 = circularVelocity;
|
||
float u0 = getU0(u1,u2,u3,r,theta,a,Delta,Sigma,-1);
|
||
float absoluteAngularDisplacement = Time * circularVelocity * -sign(a);
|
||
float absoluteVerticalDisplacement = Time * verticalVelocity * -cos(cameraRThetaphi.y);
|
||
float delayAngular = (circularVelocity / u0) * t;
|
||
float delayVertical = (cos(cameraRThetaphi.y)*(verticalVelocity / u0)) * t;
|
||
float delayedAngularDisplacement = absoluteAngularDisplacement - delayAngular;
|
||
float delayedVerticalDisplacement = absoluteVerticalDisplacement - delayVertical;
|
||
return set(delayedAngularDisplacement,delayedVerticalDisplacement);
|
||
}
|
||
|
||
//Redshift Function
|
||
function float getRedshiftFactor(vector rayOrig; vector4 cameraVel; float a; float r; float theta; float phi; float Delta; float Sigma; vector cameraRThetaphi; vector FirstU; float u0; float u1; float u2; float u3; float FirstU0; float Jet; float JetVerticVel; float JetCircVel; int isFinalRay)
|
||
{
|
||
//Photon Momentum
|
||
float gtPhi = -((2*a*1.0*r) / Sigma) * pow(sin(theta),2);
|
||
float gtt = -(1-((2*1.0*r) / Sigma));
|
||
float grr = Sigma / Delta;
|
||
float gThetaTheta = Sigma;
|
||
float Lambda = pow(pow(r,2) + pow(a,2),2) - pow(a,2)*Delta*pow(sin(theta),2);
|
||
float gPhiPhi = (Lambda/Sigma) * pow(sin(theta),2);
|
||
float DiskVel, Vel0;
|
||
vector Vel;
|
||
|
||
if(isFinalRay != 1)
|
||
{
|
||
//Disk Velocity
|
||
DiskVel = sign(a)*sqrt(1.0*r) / (r * sqrt( pow(r,2) - (3*1.0*r) + (2*abs(a)*sqrt(1.0*r)) ));
|
||
Vel = set(0,0,DiskVel);
|
||
|
||
if(Jet == 1)
|
||
{
|
||
float ux = JetCircVel * sin(phi);
|
||
float uy = -JetVerticVel * sign(rayOrig.y);
|
||
float uz = -JetCircVel* cos(phi);
|
||
vector u = set(ux,uy,uz);
|
||
float ur = (1.0 / sqrt(pow(r,2)+pow(a,2)*pow(cos(theta),2))) * ( r*sin(theta)*cos(phi)*u.x + r*sin(theta)*sin(phi)*u.z + sqrt(pow(r,2)+pow(a,2))*cos(theta)*u.y );
|
||
float uTheta = (1.0 / sqrt(pow(r,2)+pow(a,2)*pow(cos(theta),2))) * ( sqrt(pow(r,2)+pow(a,2))*cos(theta)*cos(phi)*u.x + sqrt(pow(r,2)+pow(a,2))*cos(theta)*sin(phi)*u.z - r*sin(theta)*u.y );
|
||
float uPhi = JetCircVel;
|
||
|
||
Vel = set(ur,uTheta,uPhi);
|
||
}
|
||
|
||
if(gtt < 0)
|
||
{
|
||
Vel0 = -(gtPhi / gtt) * Vel.z + sqrt( pow((gtPhi / gtt),2) * pow(Vel.z,2) - (1.0/gtt) * (grr*pow(Vel.x,2)+gThetaTheta*pow(Vel.y,2)+gPhiPhi*pow(Vel.z,2) + 1));
|
||
}
|
||
else
|
||
{
|
||
Vel0 = -(gtPhi / gtt) * Vel.z - sqrt( pow((gtPhi / gtt),2) * pow(Vel.z,2) - (1.0/gtt) * (grr*pow(Vel.x,2)+gThetaTheta*pow(Vel.y,2)+gPhiPhi*pow(Vel.z,2) + 1));
|
||
}
|
||
}
|
||
else
|
||
{
|
||
Vel = 0;
|
||
Vel0 = 1;
|
||
r *= 1e50;
|
||
}
|
||
|
||
float freq_emitter = gtt*u0*Vel0 + gtPhi * (u0*Vel.z + u3*Vel0) + grr*u1*Vel.x + gThetaTheta*u2*Vel.y + gPhiPhi*u3*Vel.z;
|
||
float rCam = cameraRThetaphi.x;
|
||
float thetaCam = cameraRThetaphi.y;
|
||
float phiCam = cameraRThetaphi.z;
|
||
float deltaCam = pow(rCam,2) - (2*1.0*rCam) + pow(a,2);
|
||
float sigmaCam = pow(rCam,2) + pow(a,2) * pow(cos(thetaCam),2);
|
||
float gtPhi_cam = -((2*a*1.0*rCam) / sigmaCam) * pow(sin(thetaCam),2);
|
||
float gtt_cam = -(1-((2*1.0*rCam) / sigmaCam));
|
||
float grr_cam = sigmaCam / deltaCam;
|
||
float gThetaTheta_cam = sigmaCam;
|
||
float gPhiPhi_cam = (pow(rCam,2) + pow(a,2) + ((2*1.0*rCam*pow(a,2)) / sigmaCam) * pow(sin(thetaCam),2) ) * pow(sin(thetaCam),2);
|
||
|
||
//Final
|
||
vector CameraVelocity = -set(cameraVel.y,cameraVel.z,cameraVel.w);
|
||
float CameraV0 = cameraVel.x;
|
||
float freq_cam = gtt_cam*FirstU0*CameraV0 + gtPhi_cam * (FirstU0*CameraVelocity.z + FirstU.z * CameraV0) + grr_cam * FirstU.x * CameraVelocity.x + gThetaTheta_cam*FirstU.y*CameraVelocity.y + gPhiPhi_cam*FirstU.z*CameraVelocity.z;
|
||
return freq_emitter/freq_cam;
|
||
}
|
||
|
||
//Disk Time Dilation Factor
|
||
function float getTimeDilationFactor(float a; float r; float theta; float Sigma; float Delta; float vel3)
|
||
{
|
||
//Photon Momentum
|
||
float gtPhi = -((2*a*1.0*r) / Sigma) * pow(sin(theta),2);
|
||
float gtt = -(1-((2*1.0*r) / Sigma));
|
||
float gThetaTheta = Sigma;
|
||
float Lambda = pow(pow(r,2) + pow(a,2),2) - pow(a,2)*Delta*pow(sin(theta),2);
|
||
float gPhiPhi = (Lambda/Sigma) * pow(sin(theta),2);
|
||
|
||
//Disk velocity and time dilation
|
||
float uPhi = vel3;
|
||
float timeDilationFactor;
|
||
|
||
if(gtt < 0)
|
||
{
|
||
timeDilationFactor = (-gtPhi*uPhi - sqrt(pow(gtPhi*uPhi,2) - gtt * (1 + gPhiPhi * pow(uPhi,2)))) / gtt;
|
||
}
|
||
else
|
||
{
|
||
timeDilationFactor = (-gtPhi*uPhi + sqrt(pow(gtPhi*uPhi,2) - gtt * (1 + gPhiPhi * pow(uPhi,2)))) / gtt;
|
||
}
|
||
|
||
return timeDilationFactor;
|
||
}
|
||
|
||
//Transform from ZAMO to Global Reference frame
|
||
function vector4 getZAMOtoGlobal(vector4 rayDir; float r; float theta; float phi; float Delta; float Sigma; float a)
|
||
{
|
||
float rayr = (1.0 / sqrt((r*r) + (a*a) * pow(cos(theta),2))) * ( (r*sin(theta)*cos(phi)*rayDir.x) + (r*sin(theta)*sin(phi)*rayDir.z) + (sqrt((r*r) + (a*a)) * cos(theta) * rayDir.y) );
|
||
float raytheta = (1.0 / sqrt((r*r) + (a*a) * pow(cos(theta),2))) * ( (sqrt((r*r) + (a*a)) * cos(theta) * cos(phi) * rayDir.x) + (sqrt((r*r) + (a*a)) * cos(theta) * sin(phi) * rayDir.z) - (r*sin(theta)*rayDir.y) );
|
||
float rayphi = (-sin(phi) * rayDir.x) + (cos(phi) * rayDir.z);
|
||
float rayt = rayDir.w;
|
||
float w2 = (r*r) + (a*a);
|
||
float Lambda = (w2*w2) - (a*a * Delta * pow(sin(theta),2));
|
||
float u0 = sqrt( Lambda / (Delta*Sigma) ) * rayt + (((2*a*1.0*r) / sqrt(Lambda*Delta*Sigma)) * rayphi);
|
||
float u1 = sqrt( Delta / Sigma ) * rayr;
|
||
float u2 = (1.0 / sqrt(Sigma)) * raytheta;
|
||
float u3 = (sqrt( Sigma / Lambda ) * (1.0 / sin(theta)) * rayphi);
|
||
return set(u0,u1,u2,u3);
|
||
}
|
||
|
||
//2D Rotation Function
|
||
function vector get2DRotation(vector A; float Angle; float AngularDisplacement)
|
||
{
|
||
vector2 A2 = set(A.x,A.z);
|
||
float s = sin(Angle+AngularDisplacement);
|
||
float c = cos(Angle+AngularDisplacement);
|
||
matrix2 m = set(c,-s,s,c);
|
||
A2 *= m;
|
||
return set(A2.x,A.y,A2.y);
|
||
}
|
||
|
||
//BL Coordinates conversion
|
||
function vector getCartesiantoBL(vector rayOrig; float a)
|
||
{
|
||
float w = (pow(rayOrig.x,2) + pow(rayOrig.y,2) + pow(rayOrig.z,2)) - (pow(a,2));
|
||
float r = sqrt(0.5 * (w + sqrt(pow(w,2) + (4 * (pow(a,2) * pow(rayOrig.y,2))))));
|
||
float theta = acos(rayOrig.y / r);
|
||
float phi = atan2(rayOrig.z,rayOrig.x);
|
||
return set(r,theta,phi);
|
||
}
|
||
|
||
//Relativistc Aberration of Light
|
||
function vector4 getRestframeToZAMO(vector rayDir; float rCamera; float thetaCamera; float DeltaCamera; float SigmaCamera; float phiCamera; float u0Camera; float u1Camera; float u2Camera; float u3Camera; float a)
|
||
{
|
||
float w2 = (rCamera*rCamera) + (a*a);
|
||
float Lambda = (w2*w2) - (a*a * DeltaCamera * pow(sin(thetaCamera),2));
|
||
float ut = sqrt((DeltaCamera*SigmaCamera) / Lambda) * u0Camera;
|
||
float ur = sqrt(SigmaCamera / DeltaCamera) * u1Camera;
|
||
float uTheta = sqrt(SigmaCamera) * u2Camera;
|
||
float uPhi = (sin(thetaCamera) * sqrt(Lambda / SigmaCamera) * u3Camera) - (((2*a*1.0*rCamera*sin(thetaCamera)) / sqrt(Lambda*SigmaCamera)) * u0Camera);
|
||
float vr = ur / ut;
|
||
float vTheta = uTheta / ut;
|
||
float vPhi = uPhi / ut;
|
||
float vX = (1.0 / sqrt((rCamera*rCamera) + (a*a) * pow(cos(thetaCamera),2))) * ( ((rCamera*sin(thetaCamera)*cos(phiCamera)*vr) + (sqrt((rCamera*rCamera) + (a*a))*cos(thetaCamera)*cos(phiCamera)*vTheta)) ) - (sin(phiCamera)*vPhi);
|
||
float vY = (1.0 / sqrt((rCamera*rCamera) + (a*a) * pow(cos(thetaCamera),2))) * ( (sqrt((rCamera*rCamera) + (a*a)) * cos(thetaCamera) * vr) - (rCamera*sin(thetaCamera)*vTheta) );
|
||
float vZ = (1.0 / sqrt((rCamera*rCamera) + (a*a) * pow(cos(thetaCamera),2))) * ( (rCamera*sin(thetaCamera)*sin(phiCamera)*vr) + (sqrt((rCamera*rCamera) + (a*a)) * cos(thetaCamera)*sin(phiCamera)*vTheta) ) + (cos(phiCamera)*vPhi);
|
||
vector v = set(vX,vY,vZ);
|
||
float v2 = v.x*v.x + v.y*v.y + v.z*v.z;
|
||
float gamma = 1.0 / sqrt(1 - v2);
|
||
matrix mCamera = set(gamma , -gamma*v.x , -gamma*v.y , -gamma*v.z , -gamma*v.x , 1+(gamma - 1)*((v.x*v.x) / (v2)) , (gamma-1)*((v.x*v.y) / v2) , (gamma-1)*((v.x*v.z) / v2) , -gamma*v.y , (gamma-1)*((v.x*v.y) / v2) , 1+(gamma - 1)*((v.y*v.y) / (v2)) , (gamma-1)*((v.y*v.z) / v2) , -gamma*v.z , (gamma-1)*((v.x*v.z) / v2) , (gamma-1)*((v.y*v.z) / v2) , 1+(gamma - 1)*((v.z*v.z) / (v2)));
|
||
vector4 ray4 = set(length(rayDir) , rayDir.x, rayDir.y, rayDir.z);
|
||
ray4 *= mCamera;
|
||
return set(ray4.y,ray4.z,ray4.w,ray4.x);
|
||
}
|
||
|
||
//Jet Density Field Function
|
||
function float getJetDensity(vector Pos; float a; float AngularDisplacement; float verticalDisplacement)
|
||
{
|
||
float h1 = abs(Pos.y);
|
||
float N0 = 1.05; //Outer Parabolla Exponent
|
||
float N1 = 6; //Inner Parabolla Exponent
|
||
float N2 = 2; //Gradient Exponent
|
||
float c0 = 9; //Outer equatorial plane Radius
|
||
float c1 = 0; //Inner equatorial plane Radius
|
||
float Amp = 0.5; //Gradient Amplitude
|
||
|
||
float x = length(Pos * set(1,0,1)) + 1;
|
||
float xh0 = pow((pow(c0,N0)+h1),1.0 / N0); //Outer Radius
|
||
float xh1 = pow((pow(c1,N1)+h1),1.0 / N1); //Inner Radius
|
||
|
||
//In order to have a correct looking jet, we need to ensure that the noise pattern does not "slide" through the jet
|
||
//but follows its curvature and distorts as the jet expands. To do this we construct a parabolic coordinate system defined
|
||
//between the jets two parabolic boundaries.
|
||
vector pPos = get2DRotation(Pos,h1*sign(a)*0.1,0);
|
||
|
||
float rDot = (x-xh1) / (xh0-xh1);
|
||
float phi = atan2(pPos.x,pPos.z) - (AngularDisplacement*sign(a));
|
||
float xC = rDot*sin(phi);
|
||
float yC = (Pos.y * 0.15) - verticalDisplacement;
|
||
float zC = rDot*cos(phi);
|
||
vector xyzC = set(xC,yC,zC);
|
||
float nNoise1 = abs(hscript_turb(xyzC,5))+0.1;
|
||
|
||
float xDot = x - nNoise1 * 5 * (-exp(-h1*0.2)+1);
|
||
|
||
if(xDot < xh0 && xDot > xh1)
|
||
{
|
||
float k = abs((xh0 / 2.0) - (xh1 / 2.0));
|
||
float xh2 = (xh0 + xh1) / 2.0;
|
||
float gx = pow((-((abs(xDot-xh2))-k)),N2)*abs(1.0 / pow(k,N2))*Amp*(exp(-h1*0.075));
|
||
float heightFallOff = exp(-h1 * 0.175);
|
||
float nNoise2 = abs(hscript_turb(xyzC * 2,5))+0.1;
|
||
return nNoise1 * nNoise2 * gx * heightFallOff;
|
||
}
|
||
else
|
||
{
|
||
return 0;
|
||
}
|
||
}
|
||
|
||
//Disk Density Field Function
|
||
function vector2 getDiskDensity(vector sampleOrig; float r; float variableNoise; float constantNoise)
|
||
{
|
||
//Constant Texture
|
||
float constNoise = constantNoise;
|
||
|
||
//Slope displacement
|
||
vector UpVector = set(0,sign(sampleOrig.y),0);
|
||
vector DirVector = normalize(sampleOrig);
|
||
float sampleAngle = degrees(dot(UpVector,DirVector));
|
||
float angleFallOff = pow(-exp(-1*r) + 1,2);
|
||
float diskSlope = 5;
|
||
float diskSlopeDot = (diskSlope + constNoise*7) * angleFallOff;
|
||
|
||
//In-Exact Bounds test
|
||
if(sampleAngle < diskSlopeDot)
|
||
{
|
||
//Invert Variable Noise
|
||
float inverseVariable = fit(variableNoise,0,1,1,0);
|
||
|
||
//Gradients
|
||
float OutterEmissionFalloff = 1.0 / (1 + exp(0.1*(r-44)));
|
||
float OutterAbsorptionFalloff = 1.0 / (1 + exp(0.1*(r-64)));
|
||
float InnerAbsorptionFalloff = 1.0 / (1 + exp(-0.1*(r-32)));
|
||
float N_0 = 5, N_1 = 2, h = tan(radians(diskSlopeDot))*r, x = abs(sampleOrig.y);
|
||
float gx = pow(-pow(x,N_0)*(1.0/pow(h,N_0))+1,N_1);
|
||
|
||
//Fields
|
||
float EmissionField = (pow(variableNoise,2) + variableNoise + 0.1) * 0.5 * OutterEmissionFalloff * gx;
|
||
float AbsorptionField = (pow(inverseVariable,3)*2 + inverseVariable + 0.1) * OutterAbsorptionFalloff * gx * 0.1;
|
||
return set(EmissionField,AbsorptionField);
|
||
}
|
||
else
|
||
{
|
||
return 0;
|
||
}
|
||
}
|
||
|
||
//Temperature Distribution
|
||
function float getTemperatureDistribution(float T1; float T2; float r; float r1; float r2; float N)
|
||
{
|
||
return pow(abs(r-r2),N)*((T1-T2) / pow(abs(r1-r2),N))+T2;
|
||
}
|
||
|
||
//Integrator and Render Settings
|
||
int integrationDepth = 2400;
|
||
int superSamples = 1;
|
||
vector superSampleCol = 0;
|
||
float drawDistance = 400;
|
||
|
||
//Supersampling Loop
|
||
for(int j = 0; j <= superSamples; j++)
|
||
{
|
||
//Random vector per pixel & sample which stays within the boundary of each pixel, so
|
||
//samples dont overlap entirly and each sample adds more information to the render.
|
||
vector sampleNear_1 = point(0,"P",1), sampleNear_2 = point(0,"P",2);
|
||
float Ndistance = distance(sampleNear_1,sampleNear_2);
|
||
|
||
//Error handeling for the last pixel, i.e. the one with the highest index (ptnum)
|
||
if(@ptnum == npoints(0))
|
||
{
|
||
Ndistance = 1.0 / Ndistance;
|
||
}
|
||
|
||
vector temp_randVec = nrandom(), randVec = fit(temp_randVec,0,1,-1,1);
|
||
randVec = normalize(set(randVec.x,randVec.y,abs(randVec.z)));
|
||
|
||
//Black Hole Settings
|
||
float a = -0.999; //Spin Parameter, this is variable and can go between -1 and 1 safely.
|
||
|
||
//Integration Related (stepSize is the base step distance, hStep is related to the adaptive step size. Change neither)
|
||
float hStep = 1e-2, hCamera = 1e-3, integrationMomenta = 0.05;
|
||
|
||
//Time
|
||
float ShutterSpeed = 1.0 / 24.0; //How long it takes the virtual camera to fully expose the image, small means less motion blur
|
||
float deltaFrameTime = (ShutterSpeed / float(superSamples)) * (float(j)+1); //Physical Motion Blur
|
||
float frameTime = 0 + deltaFrameTime; //Local time measured by the camera, if at 0 the camera will be quasi-stationary
|
||
float globalTime = $F + deltaFrameTime; //Time measured by an observer infinitely far away, used for disk and jet animation
|
||
|
||
//The Camera is not stationary in spacetime and will follow a Geodesic as time advances
|
||
float FocalLength = 36;
|
||
vector fragCoord = v@P;
|
||
vector camPos = set(0,60,300);
|
||
|
||
//Initials for camera Integration
|
||
vector getBLCam = getCartesiantoBL(camPos,a);
|
||
float rCamera = getBLCam.x, thetaCamera = getBLCam.y, phiCamera = getBLCam.z, geodesicTime = 0, DeltaCamera = pow(rCamera,2) - (2*rCamera) + pow(a,2), SigmaCamera = pow(rCamera,2) + pow(a,2) * pow(cos(thetaCamera),2);
|
||
vector InitialR_Theta_phi_Camera = set(rCamera,thetaCamera,phiCamera);
|
||
|
||
//Camera momentum
|
||
float u1Camera = 0; //controls the radial motion. Positive is away from the BH and negative is towards it.
|
||
float u2Camera = 0; //controls the motion in the theta direction, which is the up and down directions when in the equatorial plane
|
||
float u3Camera = 0;//controls the motion in the phi direction, which is the tangentian directions in the equatorial plane
|
||
u2Camera *= 1.0 / (2*M_PI*rCamera); u3Camera *= 1.0 / (2*M_PI*rCamera);
|
||
float u0Camera = getU0(u1Camera, u2Camera, u3Camera, rCamera, thetaCamera, a, DeltaCamera, SigmaCamera, -1), u0InitialCamera = u0Camera;
|
||
vector FirstUCamera = set(u1Camera,u2Camera,u3Camera), ExportUCamera = 0, LastPosCamera = 0, lastOrigCamera = 0;
|
||
|
||
//Variables
|
||
matrix3 updateVariablesCamera = 0;
|
||
float tUCamera,rUCamera,thetaUCamera,phiUCamera,u0UCamera,u1UCamera,u2UCamera,u3UCamera,hStepUCamera,XCamera,YCamera,ZCamera;
|
||
|
||
//Due to time dilation and the changing step size, one simulation step =/= to a time step. This while loop simulates the camera free falling and only stops when the local measured time (geodesicTime) is greater than the framTime. This ensures the temporal distance between frames is always eqal
|
||
while(frameTime > geodesicTime)
|
||
{
|
||
updateVariablesCamera = getRKF45step(geodesicTime,rCamera,thetaCamera,phiCamera,u0Camera,u1Camera,u2Camera,u3Camera,hCamera,a,0,0);
|
||
assign(tUCamera,rUCamera,thetaUCamera,phiUCamera,u0UCamera,u1UCamera,u2UCamera,u3UCamera,hStepUCamera,updateVariablesCamera);
|
||
geodesicTime = tUCamera; rCamera = rUCamera; thetaCamera = thetaUCamera; phiCamera = phiUCamera; u0Camera = u0UCamera; u1Camera = u1UCamera; u2Camera = u2UCamera; u3Camera = u3UCamera; hCamera = hStepUCamera*0.025;
|
||
DeltaCamera = pow(rCamera,2) - (2*rCamera) + pow(a,2); SigmaCamera = pow(rCamera,2) + pow(a,2) * pow(cos(thetaCamera),2);
|
||
XCamera = sqrt(pow(rCamera,2) + pow(a,2)) * sin(thetaCamera) * cos(phiCamera); YCamera = rCamera * cos(thetaCamera); ZCamera = sqrt(pow(rCamera,2) + pow(a,2)) * sin(thetaCamera) * sin(phiCamera);
|
||
lastOrigCamera = camPos; camPos = set(XCamera,YCamera,ZCamera); LastPosCamera = camPos;
|
||
}
|
||
|
||
//Ray Origin, Ray Dir calculation and orientation
|
||
vector rayOrig = set(0,0,(FocalLength/1000.0)) + camPos;
|
||
vector rayDir = normalize(((fragCoord + camPos) - rayOrig) + (randVec * Ndistance * 0.5)) * integrationMomenta;
|
||
vector viewDir = normalize(camPos - rayOrig);
|
||
vector rotAxisX = set(1,0,0);
|
||
vector rotAxisY = set(0,1,0);
|
||
vector rotAxisZ = set(0,0,1);
|
||
float AngleX = radians(-11);
|
||
float AngleY = radians(0);
|
||
float AngleZ = radians(0);
|
||
matrix3 mRotX = ident(), mRotY = ident(), mRotZ = ident();
|
||
rotate(mRotX,AngleX,rotAxisX); rotate(mRotY,AngleY,rotAxisY); rotate(mRotZ,AngleZ,rotAxisZ);
|
||
rayDir *= mRotX * mRotY * mRotZ;
|
||
|
||
//Convert initial coordinates into Boyer-Lindquist
|
||
vector getBl = getCartesiantoBL(rayOrig,a);
|
||
float r = getBl.x, theta = getBl.y, phi = getBl.z, t = 0, Delta = pow(r,2) - (2*r) + pow(a,2), Sigma = pow(r,2) + pow(a,2) * pow(cos(theta),2);
|
||
|
||
//Relativistic Aberration
|
||
vector4 ZAMOrayDir = getRestframeToZAMO(rayDir, rCamera, thetaCamera, DeltaCamera, SigmaCamera, phiCamera, u0Camera, u1Camera, u2Camera, u3Camera, a);
|
||
|
||
/*v@P = rayOrig;
|
||
v@N = set(ZAMOrayDir.x,ZAMOrayDir.y,ZAMOrayDir.z);
|
||
break;*/
|
||
|
||
//Convert Initial rayDir into Photon Momentum
|
||
vector4 photonMomentum = getZAMOtoGlobal(ZAMOrayDir,r,theta,phi,Delta,Sigma,a);
|
||
float u0 = photonMomentum.x, u1 = photonMomentum.y, u2 = photonMomentum.z, u3 = photonMomentum.w, u0Initial = u0;
|
||
|
||
//Volume related (Variable)
|
||
float sigma_s = 0.5, sigma_a = 0.5;
|
||
vector2 diskDensity = 0;
|
||
vector DiskBounds = set(120,0,radians(7)), JetBounds = set(36,56,0), sampleOrig = 0, L = 0, Le = 0, sBl = 0, xPos = 0, xBl = 0, xPosDot = 0, sPosDot = 0;
|
||
|
||
//Volume related (Constant)
|
||
float DiskEmissionDensity = 0, JetEmissionDensity = 0, DiskAbsorptionDensity = 0, JetAbsorptionDensity = 0, sigma_t = sigma_a + sigma_s, TSample = 0, Tvol = 1, totalMedium = 1, sXF1 = 0, sXF2 = 0, sXF3 = 0, VF1 = 0, VF2 = 0, VF3 = 0, xr = 0, weight = 0, sWeightF1 = 0, sWeightF2 = 0, sWeightF3 = 0;
|
||
int Octaves = 8;
|
||
|
||
//Render related
|
||
int inBB = 0;
|
||
float jitter = 0, dx = 1, lightTravelDelayExaggeration = 4, animationScale = 10;
|
||
|
||
//Physics related (Variable)
|
||
vector2 jetVelocity = set(1,1.5) * 5; //jetVelocity Refers to the anmation speed
|
||
|
||
//Physics related (Constant)
|
||
float timeDilationFactor = 0, X = 0, Y = 0, Z = 0, tU, rU, thetaU, phiU, u0U, u1U, u2U, u3U, hStepU, vel3, Temperature = 0, redshift = 0, PhysicalPhiVelocity = 0, UnphysicalPhiVelocity = 0, JetCircularVelocity = 0, JetVerticalVelocity = 0; //These velocities refer to the redshift
|
||
vector2 tDelta = 0, jetMotion;
|
||
vector LastPos = 0, FirstU = set(u1,u2,u3), InitialRThetaPhi = set(r,theta,phi);
|
||
vector4 spectralBlackbody = 0;
|
||
matrix3 updateVariables = 0;
|
||
|
||
//Render Metadata
|
||
f@FocalLength = FocalLength; v@camPos = camPos; f@cameraU0 = u0Camera; f@cameraU2 = u1Camera; f@cameraU3 = u2Camera; f@cameraU4 = u3Camera; v@orientation = set(AngleX,AngleY,AngleZ);
|
||
|
||
//Main Render Loop
|
||
for(int i = 0; i <= integrationDepth; i++)
|
||
{
|
||
//Break Condition - Event Horizion (1+sqrt(etc) is the size of the Horizon, if a ray is this close it will fall into the Horizon)
|
||
if(r < (1+sqrt(1-pow(abs(a),2)))*1.0001)
|
||
{
|
||
superSampleCol += L;
|
||
f@T = 0; v@rayOrig = rayOrig; v@rayDir = rayDir; f@Z = 1;
|
||
break;
|
||
}
|
||
|
||
//Break Condition - Celestial Sphere (i.e the ray is far enough away for the BH to have no influence anymore)
|
||
if(r > drawDistance)
|
||
{
|
||
superSampleCol += L;
|
||
f@T = Tvol; v@rayOrig = rayOrig; v@rayDir = rayDir; f@Z = getRedshiftFactor(sampleOrig,set(u0Camera,u1Camera,u2Camera,u3Camera),a,r,theta,phi,Delta,Sigma+1,InitialRThetaPhi,FirstU,u0,u1,u2,u3,u0Initial,0,0,0,1);
|
||
break;
|
||
}
|
||
|
||
//Error Handeling - Terminates loop if the integration does not finish
|
||
if(i+1 == integrationDepth)
|
||
{
|
||
superSampleCol += L;
|
||
f@Error = 1; f@T = Tvol; v@rayOrig = rayOrig; v@rayDir = rayDir; f@Z = 1;
|
||
break;
|
||
}
|
||
|
||
//In Volume Bounding region check
|
||
if(r < DiskBounds.x + 5)
|
||
{
|
||
//In-jet bounds region check
|
||
if(sqrt(r*r + a*a)*sin(theta) < JetBounds.x && abs(rayOrig.y) < JetBounds.y)
|
||
{
|
||
//In Bound Tag
|
||
inBB = 1;
|
||
|
||
//Sample Location
|
||
jitter = (nrandom()-0.5) * dx;
|
||
sampleOrig = rayOrig + (normalize(rayDir)*jitter);
|
||
|
||
//Do Volume Calculations
|
||
jetMotion = getJetDisplacement(InitialRThetaPhi, globalTime, a, r, theta, t, Delta, Sigma, jetVelocity.x, jetVelocity.y);
|
||
JetEmissionDensity = getJetDensity(sampleOrig,a,jetMotion.x,jetMotion.y);
|
||
JetAbsorptionDensity = JetEmissionDensity * 0.5 * (exp(-1*abs(rayOrig.y))+1);
|
||
}
|
||
|
||
//In-disk bounds region check
|
||
if(sqrt(r*r + a*a)*sin(theta) < DiskBounds.x && sqrt(r*r + a*a)*sin(theta) > DiskBounds.y && abs(rayOrig.y) < tan(DiskBounds.z)*(r+24))
|
||
{
|
||
//In Bound Tag
|
||
inBB = 1;
|
||
|
||
//Sample Location
|
||
jitter = (nrandom()-0.5) * dx;
|
||
sampleOrig = rayOrig + (normalize(rayDir)*jitter);
|
||
sPosDot = get2DRotation(sampleOrig*set(1,3,1),sqrt(1.0 / pow(r+3,2))*-sign(a)*128,10);
|
||
|
||
//Vortex Noise Algorithm
|
||
float rotAngle = 0, constWeight, variableWeight, vWeight, cWeight, vValue, cValue, vortexNoiseValue, constantVortexNoiseValue, NthHarmonic = 0;
|
||
|
||
//For number of octaves
|
||
for(int o = 1; o <= Octaves; o++)
|
||
{
|
||
vWeight = 0; cWeight = 0; vValue = 0; cValue = 0;
|
||
|
||
//For N points in octave
|
||
for(int g = 0; g <= npoints(-o); g++)
|
||
{
|
||
xPos = point(-o,"P",g);
|
||
X = point(-o,"V",g);
|
||
X = (X - 0.5) * 2;
|
||
xBl = getCartesiantoBL(xPos,a);
|
||
vel3 = sign(a)*sqrt(xBl.x) / (xBl.x * sqrt( pow(xBl.x,2) - (3*xBl.x) + (2*abs(a)*sqrt(xBl.x))));
|
||
timeDilationFactor = getTimeDilationFactor(a,r,theta,Sigma,Delta,vel3 * lightTravelDelayExaggeration);
|
||
rotAngle = ((globalTime*animationScale) - (t/timeDilationFactor)) * vel3 * lightTravelDelayExaggeration;
|
||
xPosDot = get2DRotation(xPos,rotAngle,0);
|
||
variableWeight = 1.0 / pow(distance(xPosDot,sPosDot),1+o);
|
||
constWeight = 1.0 / pow(distance(set(xPosDot.x,0,xPosDot.z),set(sPosDot.x,0,sPosDot.z)),1+o);
|
||
vWeight += variableWeight; cWeight += constWeight; vValue += X*variableWeight; cValue += X*constWeight;
|
||
}
|
||
vortexNoiseValue += (vValue/vWeight)*(1.0/o);
|
||
constantVortexNoiseValue += (cValue/cWeight)*(1.0/o);
|
||
NthHarmonic += 1.0/o;
|
||
}
|
||
|
||
vortexNoiseValue /= NthHarmonic; constantVortexNoiseValue /= NthHarmonic;
|
||
diskDensity = getDiskDensity(sampleOrig, r, abs(vortexNoiseValue), abs(constantVortexNoiseValue));
|
||
DiskEmissionDensity = diskDensity.x;
|
||
DiskAbsorptionDensity = diskDensity.y;
|
||
}
|
||
|
||
//Render Volume Disk
|
||
if(DiskEmissionDensity > 0 || DiskAbsorptionDensity > 0)
|
||
{
|
||
Temperature = getTemperatureDistribution(3600, 1800, r, 0, DiskBounds.x, 5);
|
||
redshift = getRedshiftFactor(sampleOrig,set(u0Camera,u1Camera,u2Camera,u3Camera),a,r,theta,phi,Delta,Sigma+1,InitialRThetaPhi,FirstU,u0,u1,u2,u3,u0Initial,0,0,0,0);
|
||
spectralBlackbody = getBlackbodyColor(Temperature, redshift);
|
||
TSample = exp((-1.0/(redshift))*dx*DiskAbsorptionDensity*sigma_t);
|
||
Tvol *= TSample;
|
||
Le = DiskEmissionDensity * set(spectralBlackbody.x,spectralBlackbody.y,spectralBlackbody.z) * spectralBlackbody.w * sigma_a;
|
||
L += Le*(Tvol/redshift)*dx*sigma_a;
|
||
|
||
//Rest Variables
|
||
totalMedium *= exp((DiskEmissionDensity + DiskAbsorptionDensity)*-dx*0.05);
|
||
DiskEmissionDensity = 0; DiskAbsorptionDensity = 0;
|
||
}
|
||
|
||
//Render Volume Jet
|
||
if(JetEmissionDensity > 0 || JetAbsorptionDensity > 0)
|
||
{
|
||
Temperature = 3600;
|
||
PhysicalPhiVelocity = (sqrt(r) / (r * sqrt( pow(r,2) - (3*r) + (2*abs(a)*sqrt(r)) )));
|
||
UnphysicalPhiVelocity = (1.0 / (pow(r-1,64)));
|
||
JetCircularVelocity = max(PhysicalPhiVelocity,UnphysicalPhiVelocity) * sign(a);
|
||
JetVerticalVelocity = 0.03;
|
||
redshift = getRedshiftFactor(sampleOrig,set(u0Camera,u1Camera,u2Camera,u3Camera),a,r,theta,phi,Delta,Sigma+1,InitialRThetaPhi,FirstU,u0,u1,u2,u3,u0Initial,1,JetVerticalVelocity,JetCircularVelocity,0);
|
||
spectralBlackbody = getBlackbodyColor(Temperature, redshift);
|
||
TSample = exp((-1.0/(redshift))*dx*JetAbsorptionDensity*sigma_t);
|
||
Tvol *= TSample;
|
||
Le = JetEmissionDensity * set(spectralBlackbody.x,spectralBlackbody.y,spectralBlackbody.z) * spectralBlackbody.w * sigma_a;
|
||
L += Le*(Tvol/redshift)*dx*sigma_a;
|
||
|
||
//Rest Variables
|
||
totalMedium *= exp((JetEmissionDensity + JetAbsorptionDensity)*-dx);
|
||
JetEmissionDensity = 0; JetAbsorptionDensity = 0;
|
||
}
|
||
}
|
||
|
||
//March Step using Runge–Kutta–Fehlberg method
|
||
updateVariables = getRKF45step(t,r,theta,phi,u0,u1,u2,u3,hStep,a,totalMedium,inBB);
|
||
assign(tU,rU,thetaU,phiU,u0U,u1U,u2U,u3U,hStepU,updateVariables);
|
||
t = tU; r = rU; theta = thetaU; phi = phiU; u0 = u0U; u1 = u1U; u2 = u2U; u3 = u3U; hStep = hStepU;
|
||
Delta = pow(r,2) - (2*r) + pow(a,2); Sigma = pow(r,2) + pow(a,2) * pow(cos(theta),2);
|
||
X = sqrt(pow(r,2) + pow(a,2)) * sin(theta) * cos(phi); Y = r * cos(theta); Z = sqrt(pow(r,2) + pow(a,2)) * sin(theta) * sin(phi);
|
||
rayOrig = set(X,Y,Z); rayDir = normalize(rayOrig - LastPos); dx = distance(rayOrig,LastPos); LastPos = rayOrig;
|
||
}
|
||
}
|
||
|
||
//Anti Aliasing
|
||
v@Cd = superSampleCol / superSamples; |