VR Schwarzschild black hole shader (works with SPS and SPS-I)

BH.shader

Shader "Misha/Schwartzschild"
{
    Properties
    {
        [Header(Black Hole)]
        _EventHorizonRadius("Event Horizon Radius", Range(0.0, 0.5)) = 0.075

        [Header(Accretion)]
        [HDR] _AccretionColor("Color", Color) = (1,1,1,1)
        _AccretionDensity("Accretion Density", Range(0, 1000)) = 200.0
        _AccretionThickness("Accretion Thickness", Range(0, 0.2)) = 0.05
        _AccretionRadius("Accretion Radius", Range(0, 1)) = 0.5
        _AccretionInnerRadius("Accretion Inner Radius", Range(0, 1)) = 0.1
        _Opacity("Opacity", Range(0, 100)) = 0.5

        [Header(Quality)]
        _StepSize("Step Size", Range(0, 0.02)) = 0.016
        _MaxSteps("Max Steps", int) = 128
        _DepthThickness("Depth Occlusion Thickness", Range(0, 0.1)) = 0.057
        [Toggle] _UseDepth("Use Depth Occlusion", int) = 1
        _MaxRadius("Max Ray Trace Distance", Range(0.5, 10)) = 0.5
    }

    SubShader
    {
        Tags { "Queue" = "Transparent" }

        GrabPass
        {
            "_BackgroundTexture"
        }

        Pass
        {
            ZWrite Off
            ZTest Always
            Cull Front
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            #include "UnityCG.cginc"

            float3 GetCameraPosition() { return _WorldSpaceCameraPos; }

            struct v2f
            {
                float3 worldPos : TEXCOORD0;
                float2 screenPos : TEXCOORD1;
                float4 pos : SV_POSITION;
            #if defined(UNITY_VERTEX_OUTPUT_STEREO)
                UNITY_VERTEX_OUTPUT_STEREO
            #endif
            };

            struct appdata
            {
                float4 vertex : POSITION;
            #if defined(UNITY_VERTEX_INPUT_INSTANCE_ID)
                UNITY_VERTEX_INPUT_INSTANCE_ID
            #endif
            };

            UNITY_DECLARE_SCREENSPACE_TEXTURE(_BackgroundTexture);
            UNITY_DECLARE_DEPTH_TEXTURE(_CameraDepthTexture);
            
            float4 _AccretionColor;
            float _Opacity;
            float _AccretionDensity;
            float _AccretionThickness;
            float _AccretionRadius;
            float _AccretionInnerRadius;
            float _StepSize;
            float _DepthThickness;
            int _MaxSteps;
            float _BHStrength;
            float _MaxRadius;
            float _EventHorizonRadius;
            int _UseDepth;

            float3 PositionToScreen(float3 pos)
            {
                float4 clipPos = UnityObjectToClipPos(float4(pos, 1));
                return float3(ComputeGrabScreenPos(clipPos).xy / clipPos.w, clipPos.z / clipPos.w);
            }

            float4 SampleBackground(float2 uv)
            {
                return UNITY_SAMPLE_SCREENSPACE_TEXTURE(_BackgroundTexture, uv);
            }

            float4 SampleDepth(float2 uv)
            {   
                uv.y = _ProjectionParams.x * .5 + .5 - uv.y * _ProjectionParams.x;
                return SAMPLE_DEPTH_TEXTURE(_CameraDepthTexture, uv);
            }

            v2f vert(appdata v) {
                v2f o;
                UNITY_SETUP_INSTANCE_ID( v );
                UNITY_INITIALIZE_OUTPUT( v2f, o );
                UNITY_INITIALIZE_VERTEX_OUTPUT_STEREO( o );
                o.worldPos = mul(unity_ObjectToWorld, v.vertex).xyz;
                o.pos = UnityObjectToClipPos(v.vertex);
                o.screenPos = ComputeGrabScreenPos(o.pos) / o.pos.w;
                return o;
            }

            float2 SphereIntersection(float3 ro, float3 rd, float3 sphPos, float sphRad)
            {
                float3 oc = ro - sphPos;
                float b = dot(oc, rd);
                float c = dot(oc, oc) - sphRad * sphRad;
                float h = b * b - c;
                if (h < 0.0) return float2(-1.0, -1.0);
                h = sqrt(h);
                return float2(-b - h, -b + h);
            }

            float BHStrength(float radius)
            {
                return pow(radius / 0.75, 3.0);
            }

            //the black hole "force"
            float3 force(float dist, float3 p)
            {
                return -1.5 * BHStrength(_EventHorizonRadius) * p * pow(dist, -5.0) * smoothstep(1.0, 0.9, dist); //smoothstep to make space flat outside the region
            }

            float sqr(float x) { return x * x; }

            float accretionDisk(float dist, float z)
            {
                //just uniform exponential disk
                float h = _AccretionThickness;
                float R = _AccretionRadius;
                float r = _AccretionInnerRadius;
                float heigh_weight = exp(- sqr(z / h));
                float rad_weight = exp(- dist / R) * smoothstep(0.05, 1.0, 1.0 - dist / R) * smoothstep(r, r*1.2, dist);
                return _AccretionDensity * heigh_weight * rad_weight;
            }

            float InterleavedGradientNoise(float2 pixelPos)
            {
                return frac(52.9829189f * frac(0.06711056f*float(pixelPos.x) + 0.00583715f*float(pixelPos.y)));  
            }

            half4 frag(v2f i, bool isFront : SV_IsFrontFace) : SV_Target
            {
                UNITY_SETUP_STEREO_EYE_INDEX_POST_VERTEX(i);
                //this blog post is used as reference
                //https://rantonels.github.io/starless/ 

                float3 cam = GetCameraPosition();
                float3 rd = i.worldPos - cam;
                rd = normalize(mul(unity_WorldToObject, float4(rd, 0.0)).xyz);
                cam = mul(unity_WorldToObject, float4(cam, 1.0)).xyz;

                float2 t = SphereIntersection(cam, rd, float3(0, 0, 0), 0.5); //boundary of flat space
                float td = max(0.0, t.x);
                cam += td * rd;

                if(t.x >= 0.0) //if outside the ray marching region - check if region is behind the object
                {
                    float3 screenPos = PositionToScreen(cam);
                    float this_depth = SampleDepth(screenPos.xy);
                    if(screenPos.z < this_depth)
                    {
                        discard;
                    }
                }
            
                float noise = InterleavedGradientNoise(int2(i.pos.xy));
                float3 rv = rd;
                float3 ro = cam;
                ro += noise * _StepSize * rv; 

                float accretion_rho = 0.0;
                float accretion_color = 0.0;
                bool hit_horizon = false;
                bool stopped = false;
                float steps = 0.0;
                float worldToObjectScale = length(mul(unity_WorldToObject, float4(1, 1, 1, 0)).xyz);
                [loop] 
                for (int j = 0; j < _MaxSteps; j++)
                {
                 
                    float dt = _StepSize;// / length(rv); //will make integration step size constant
                    float dist = length(ro);
                    if(dist < 0.5) //is inside curved space region
                    {
                        float rho = accretionDisk(dist, ro.y);
                        accretion_rho += _Opacity*rho * dt;
                        float opacity = exp(-accretion_rho);
                        accretion_color += rho * dt * opacity;
                        //integrate velocity
                        rv += force(dist, ro) * dt;
                    }
                    //integrate position
                    ro += rv * dt;

                    if (dist < _EventHorizonRadius) { //if has hit the event horizon
                        steps = float(j);
                        hit_horizon = true; 
                        break;
                    }

                    if(dist > _MaxRadius || accretion_rho > 32.0) //if outside the region or occluded
                    {
                        steps = float(j);
                        break;
                    }

                    if(_UseDepth) //if depth occlusion is enabled
                    { 
                        float3 screenPos = PositionToScreen(ro);
                        float depth_scene = worldToObjectScale*LinearEyeDepth(SampleDepth(screenPos.xy));
                        float depth_ro = worldToObjectScale*LinearEyeDepth(screenPos.z);
                        float diff = depth_scene - depth_ro;
                        if((diff < 0.0) && (diff > -_DepthThickness)) //if has hit an object in the scene
                        {
                            steps = float(j);
                            stopped = true;
                            break;
                        }
                    }
                }

                float3 screenPos = PositionToScreen(ro);
                float depth = SampleDepth(screenPos.xy);
                
                bool behind = screenPos.z < depth;

                float4 bgcolor;
                if(_UseDepth && behind && !stopped) //if background sample occluded - fallback to no distortion
                {
                    screenPos.xy = i.screenPos;
                }

                float4 background = SampleBackground(screenPos.xy) * (1.0 - float(hit_horizon));
                float4 accretionColor = _AccretionColor * accretion_color;
                float opacity = 1.0 - exp(-accretion_rho);
                return float4(lerp(background.rgb, accretionColor.rgb, opacity), 1.0);
            }
            ENDCG
        }
    }
}