Make neutron stars matter jets volumetric

src/shaders/matterjet.glsl

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

src/ts/playground.ts

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