How to correctly project equirectangular to azimuthal in HLSL?

Question

I'm developing a small Astroids clone for a PufferSphere device in Unity URP. I'm using a cubemap to equirectangular shader for the camera, to capture the 3d sphere ingame.

How would I correctly transform the equirectangular projection to azimuthal? The desired result is explained in this thread here, but as a GL shader: How to do a shader to convert to azimuthal_equidistant

But I'm running into issues, that either the shader does nothing or the result has diamond shaped artefacts.

The shader I'm using(I marked the section I added with //equirectengular to azimuthal START):

// Upgrade NOTE: replaced 'mul(UNITY_MATRIX_MVP,*)' with 'UnityObjectToClipPos(*)'

Shader "Conversion/CubemapToEquirectangularTEST" {

    Properties{
          _MainTex("Cubemap (RGB)", CUBE) = "" {}
    }

        Subshader{
            Pass {
                ZTest Always Cull Off ZWrite Off
                Fog { Mode off }

                CGPROGRAM
                    #pragma vertex vert
                    #pragma fragment frag
                    #pragma fragmentoption ARB_precision_hint_fastest
        //#pragma fragmentoption ARB_precision_hint_nicest
        #include "UnityCG.cginc"

        #define PI    3.141592653589793
        #define TWOPI 6.283185307179587

        struct v2f {
            float4 pos : POSITION;
            float2 uv : TEXCOORD0;
        };

        samplerCUBE _MainTex;

        v2f vert(appdata_img v)
        {
            v2f o;
            o.pos = UnityObjectToClipPos(v.vertex);
            o.uv = v.texcoord.xy * float2(TWOPI, PI);
            return o;
        }

        fixed4 frag(v2f i) : COLOR
        {
            float theta = i.uv.y;
            float phi = i.uv.x;
            float3 unit = float3(0,0,0);

            unit.x = sin(phi) * sin(theta) * -1;
            unit.y = cos(theta) * -1;
            unit.z = cos(phi) * sin(theta) * -1;


            //equirectengular to azimuthal START
            float a,d;
            d=length(i.pos);
            
            if (d<1.0)      // inside projected sphere surface
            {
                a=atan2(i.pos.x,i.pos.y);
                if (a<0.0) a+=TWOPI;
                if (a>TWOPI) a-=TWOPI;
                unit.x=a/TWOPI;
                unit.y=d;
            }
            //equirectengular to azimuthal END
            return texCUBE(_MainTex, unit);
        }
    ENDCG
}
    }
        Fallback Off
}

Current output:

Desired output: