Skip to content

Commit

Permalink
Merge pull request #309 from BarthPaleologue/268-feature-make-neutron…
Browse files Browse the repository at this point in the history 
…-stars-matter-jets-volumetric

Make neutron stars matter jets volumetric
  • Loading branch information
@BarthPaleologue
BarthPaleologue authored Feb 17, 2025
2 parents b7761f0 + 832c7b2 commit 86a031b
Show file tree
Hide file tree
Showing 10 changed files with 364 additions and 113 deletions.
285 changes: 206 additions & 79 deletions src/shaders/matterjet.glsl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,10 @@ uniform sampler2D depthSampler;// the depth map of the camera

uniform float time;

uniform mat4 inverseRotation;

uniform float dipoleTilt;

#include "./utils/camera.glsl";

#include "./utils/object.glsl";
Expand All @@ -36,79 +40,189 @@ uniform float time;

#include "./utils/removeAxialTilt.glsl";

// from https://www.shadertoy.com/view/MtcXWr
bool rayIntersectCone(vec3 rayOrigin, vec3 rayDir, vec3 tipPosition, vec3 orientation, float coneAngle, out float t1, out float t2) {
vec3 co = rayOrigin - tipPosition;

float a = dot(rayDir, orientation)*dot(rayDir, orientation) - coneAngle*coneAngle;
float b = 2. * (dot(rayDir, orientation)*dot(co, orientation) - dot(rayDir, co)*coneAngle*coneAngle);
float c = dot(co, orientation)*dot(co, orientation) - dot(co, co)*coneAngle*coneAngle;

float det = b*b - 4.*a*c;
if (det < 0.) return false;

det = sqrt(det);
t1 = (-b - det) / (2. * a);
t2 = (-b + det) / (2. * a);

// This is a bit messy; there ought to be a more elegant solution.
float t = t1;
if (t < 0. || t2 > 0. && t2 < t) t = t2;
if (t < 0.) return false;

vec3 cp = rayOrigin + t*rayDir - tipPosition;
float h = dot(cp, orientation);

vec3 n = normalize(cp * dot(orientation, cp) / dot(cp, cp) - orientation);

return true;
}

// see https://www.shadertoy.com/view/tslcW4
const float a=1.0;
const float b=.1759;
const float PI=3.14159265359;

float spiralSDF(float theta, float radius) {

float t=theta;
// t=(t+PI)/(2.*PI);
float r=radius;

float n=(log(r/a)/b-t)/(2.*PI);

// Cap the spiral
// float nm = (log(0.11)/b-t)/(2.0*PI);
// n = min(n,nm);
// return (n+1.0)/100.0;
float upper_r=a*exp(b*(t+2.*PI*ceil(n)));
float lower_r=a*exp(b*(t+2.*PI*floor(n)));
// float lower_r = 0.0;

return min(abs(upper_r-r), abs(r-lower_r));
float sdSpiral(vec3 p, float coneTheta, float frequency) {
// Y-axis is now the spiral direction
float spiralPos = p.y; // Using Y instead of Z
float radius = spiralPos * tan(coneTheta);
float theta = frequency * spiralPos;

// Spiral now wraps around Y-axis
vec3 spiralPoint = vec3(
radius * cos(theta),
spiralPos, // Y position is our spiral progress
radius * sin(theta) // Z component completes the circular motion
);

return length(p - spiralPoint);
}

float spiralDensity(vec3 pointOnCone, vec3 coneAxis, float coneMaxHeight) {
// Then we rotate that point so that we eliminate the axial tilt of the star from the equation
vec3 pointOnYCone = removeAxialTilt(pointOnCone, coneAxis);
float spiralDensity(vec3 p, float coneTheta, float coneHeight) {
float frequency = 0.1;
float dist = sdSpiral(p, coneTheta, frequency);

float heightFraction = abs(p.y) / coneHeight;

vec2 pointOnXZPlane = vec2(pointOnYCone.x, pointOnYCone.z);
float theta = atan(pointOnXZPlane.y, pointOnXZPlane.x) + 3.14 * min(0.0, sign(dot(pointOnCone, coneAxis)));
float heightFraction = abs(pointOnYCone.y) / coneMaxHeight;
dist /= heightFraction;

float density = 1.0;

// smoothstep fadeout when the height is too much (outside of cone)
density *= exp(-0.02 * dist * dist);

density *= 1.0 - smoothstep(0.0, 1.0, heightFraction);

return density;
}

float d = spiralSDF(theta + time, 0.2 + sqrt(heightFraction) / 2.0) / (0.3 + heightFraction * 2.0);
//d = pow(d, 4.0);

density *= smoothstep(0.6, 1.0, pow(1.0 - d, 8.0)) * 2.0; //smoothstep(0.85, 1.0, 1.0 - d) * 2.0;
// from https://www.shadertoy.com/view/MtcXWr
// Returns true if the ray intersects the closed cone.
// If so, 't' is set to the first hit along the ray (entry if outside, exit if inside)
// and 'distTrough' is set to the chord length of the ray inside the cone.
bool rayIntersectCone(vec3 rayOrigin, vec3 rayDir,
float coneHeight, float cosTheta,
out float t, out float distTrough)
{
vec3 coneUp = vec3(0.0, 1.0, 0.0);

float cosTheta2 = cosTheta * cosTheta;

// Compute tangent once.
float sinTheta = sqrt(max(0.0, 1.0 - cosTheta2));
float tanTheta = sinTheta / cosTheta;

// --- Determine if ray origin is inside the cone volume ---
bool inside = false;
{
vec3 v = rayOrigin;
float h = v.y;
if (h >= 0.0 && h <= coneHeight) {
float rAtH = h * tanTheta;
float d = length(v - coneUp * h);
if (d <= rAtH)
inside = true;
}
}

// --- Collect candidate intersection t's ---
// We will add valid intersections from the lateral surface and the base.
float tCandidates[3];
int count = 0;

// Intersection with the infinite cone's lateral surface.
vec3 co = rayOrigin;
float A = rayDir.y * rayDir.y - cosTheta2;
float B = 2.0 * (rayDir.y * co.y - dot(rayDir, co) * cosTheta2);
float C = co.y * co.y - dot(co, co) * cosTheta2;

float det = B * B - 4.0 * A * C;
if (det >= 0.0)
{
float sqrtDet = sqrt(det);
float t1 = (-B - sqrtDet) / (2.0 * A);
float t2 = (-B + sqrtDet) / (2.0 * A);

// Check t1 for validity (only consider if t>=0)
if (t1 >= 0.0)
{
vec3 cp1 = rayOrigin + t1 * rayDir;
float h1 = cp1.y;
if (h1 >= 0.0 && h1 <= coneHeight)
tCandidates[count++] = t1;
}
// Check t2
if (t2 >= 0.0)
{
vec3 cp2 = rayOrigin + t2 * rayDir;
float h2 = cp2.y;
if (h2 >= 0.0 && h2 <= coneHeight)
tCandidates[count++] = t2;
}
}

// Intersection with the base plane.
vec3 baseCenter = coneUp * coneHeight;
float denom = rayDir.y;
if (abs(denom) > 1e-6)
{
float tBase = (baseCenter - rayOrigin).y / denom;
if (tBase >= 0.0)
{
vec3 hitPoint = rayOrigin + tBase * rayDir;
// Check if hitPoint is within the circular base.
float baseRadius = coneHeight * tanTheta;
if (length(hitPoint - baseCenter) <= baseRadius)
tCandidates[count++] = tBase;
}
}

// If no valid intersections were found, return false.
if (count == 0)
return false;

// --- Sort the candidate intersections in increasing order ---
for (int i = 0; i < count - 1; i++)
{
for (int j = i + 1; j < count; j++)
{
if (tCandidates[j] < tCandidates[i])
{
float tmp = tCandidates[i];
tCandidates[i] = tCandidates[j];
tCandidates[j] = tmp;
}
}
}

// --- Determine entry and exit t values ---
float tEntry, tExit;
if (inside)
{
// If the ray origin is inside, entry is at t = 0 and the exit is the first candidate.
tEntry = 0.0;
tExit = tCandidates[0];
}
else
{
// Outside: entry is the first hit, exit is the second (if present).
tEntry = tCandidates[0];
if (count > 1)
tExit = tCandidates[1];
else
tExit = tEntry; // Tangential hit: chord length is zero.
}

// Set the outputs.
// We return tEntry as the “hit” point (for an outside ray, where the ray first enters the cone;
// for an inside ray, tEntry=0 is the starting point so we return the exit instead).
if (inside)
t = 0.0;
else
t = tEntry;

distTrough = tExit - tEntry;

return true;
}

//density *= d * 500.0;
vec3 rayMarchSpiral(vec3 rayOrigin, vec3 rayDir, float distThrough, int nbSteps, float coneTheta, float coneHeight, out float transmittance) {
vec3 col = vec3(0.0);

transmittance = 1.0;

float stepSize = distThrough / float(nbSteps);

for(int i = 0; i < nbSteps; i++) {
vec3 p = rayOrigin + float(i) * stepSize * rayDir;

float density = spiralDensity(p, coneTheta * 0.7, coneHeight);

vec3 emission = 3.0 * vec3(0.3, 0.6, 1.0) * density;
float absorption = 0.2 * density;

transmittance *= exp(-absorption * stepSize);
col += emission * transmittance * stepSize;
}

return density;
return col;
}

void main() {
Expand All @@ -123,29 +237,42 @@ void main() {

vec3 rayDir = normalize(pixelWorldPosition - camera_position);// normalized direction of the ray

vec4 finalColor = screenColor;
float scaling_factor = object_radius * 10000.0;

const float jetHeight = 10000000e3;
const vec3 jetColor = vec3(0.5, 0.5, 1.0);
// move in object's local space to simplify the calculations
vec3 rayOriginLocalSpace = mat3(inverseRotation) * (camera_position - object_position) / scaling_factor;
vec3 rayDirLocalSpace = mat3(inverseRotation) * rayDir;

float coneTheta = dipoleTilt;
float coneHeight = 100.0;

float t1, t2;
if (rayIntersectCone(camera_position, rayDir, object_position, object_rotationAxis, 0.95, t1, t2)) {
if (t2 > 0.0 && t2 < maximumDistance) {
vec3 jetPointPosition2 = camera_position + t2 * rayDir - object_position;
vec3 color = screenColor.rgb;
float finalAlpha = screenColor.a;

float density2 = spiralDensity(jetPointPosition2, object_rotationAxis, jetHeight);
float t, distThrough;
if(rayIntersectCone(rayOriginLocalSpace, rayDirLocalSpace, coneHeight, cos(coneTheta), t, distThrough) && t * scaling_factor < maximumDistance) {
vec3 startPoint = rayOriginLocalSpace + t * rayDirLocalSpace;

finalColor.rgb = mix(finalColor.rgb, jetColor, density2);
}
if (t1 > 0.0 && t1 < maximumDistance) {
vec3 jetPointPosition1 = camera_position + t1 * rayDir - object_position;
float transmittance;
vec3 jetColor = rayMarchSpiral(startPoint, rayDirLocalSpace, distThrough, 100, coneTheta, coneHeight, transmittance);

float density1 = spiralDensity(jetPointPosition1, object_rotationAxis, jetHeight);
color = mix(jetColor, color, transmittance);
finalAlpha = max(transmittance, finalAlpha);
}

finalColor.rgb = mix(finalColor.rgb, jetColor, density1);
}
// flip the coordinates to display the other cone
rayOriginLocalSpace *= -1.0;
rayDirLocalSpace *= -1.0;

if(rayIntersectCone(rayOriginLocalSpace, rayDirLocalSpace, coneHeight, cos(coneTheta), t, distThrough) && t * scaling_factor < maximumDistance) {
vec3 startPoint = rayOriginLocalSpace + t * rayDirLocalSpace;

float transmittance;
vec3 jetColor = rayMarchSpiral(startPoint, rayDirLocalSpace, distThrough, 100, coneTheta, coneHeight, transmittance);

color = mix(jetColor, color, transmittance);
finalAlpha = max(transmittance, finalAlpha);
}

gl_FragColor = finalColor;// displaying the final color
gl_FragColor = vec4(color, finalAlpha);
}
4 changes: 4 additions & 0 deletions src/ts/playground.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ import { createDebugAssetsScene } from "./playgrounds/debugAssets";
import { createSpaceStationScene } from "./playgrounds/spaceStation";
import { createXrScene } from "./playgrounds/xr";
import { createFlightDemoScene } from "./playgrounds/flightDemo";
import { createNeutronStarScene } from "./playgrounds/neutronStar";

const canvas = document.getElementById("renderer") as HTMLCanvasElement;
canvas.width = window.innerWidth;
Expand Down Expand Up @@ -65,6 +66,9 @@ switch (requestedScene) {
case "flightDemo":
scene = await createFlightDemoScene(engine);
break;
case "neutronStar":
scene = await createNeutronStarScene(engine);
break;
default:
scene = await createAutomaticLandingScene(engine);
}
Expand Down
Loading

0 comments on commit 86a031b

Please sign in to comment.