Artifacts in atmospheric scattering shader

I am trying to implement the atmospheric scattering shader found here: https://lightshaderdevlog.wordpress.com/

I have implemented the algorithm in GLSL and it looks good, however for some reason black artifacts render: https://i.sstatic.net/OZ6t1.jpg

It seems that the artifacts are only present when looking at specific camera angles and in those specific vertices.

Here is my code:

Vertex Shader

 #version 330 core
layout (location = 0) in vec3 aPos;
layout (location = 1) in vec2 aTexCoord;
layout (location = 2) in vec3 aNormal;

// Vertex
uniform vec3 star_pos;      // point light source position (Sun) [GCS]
uniform vec3 planet_pos;    // planet center position [GCS]
uniform float overglow;     // cos(angle) of terminator propagation\

uniform mat4 modelView;
uniform mat4 projection;
uniform float og_farPlaneDistance;
uniform float u_logarithmicDepthConstant;

out vec3 lpos;          // point light source position (Sun) [camera space]
out vec3 ppos;          // planet center position [camera space]
out vec3 pixel_pos;     // fragment position [camera space]
out vec3 pixel_nor;     // fragment surface normal
out vec2 pixel_txy;     // fragment surface texture coord
out float angleIncidence;
out vec4 colAtmosphere;
out vec3 lightDir;

const float PI = 3.14159265f;
const float transitionWidth = 0.1f;
const float fresnelExponent = 20f;


vec4 modelToClipCoordinates(vec4 position, mat4 modelViewPerspectiveMatrix, float depthConstant, float farPlaneDistance){
    vec4 clip = modelViewPerspectiveMatrix * position;

    clip.z = ((2.0 * log(depthConstant * clip.z + 1.0) / log(depthConstant * farPlaneDistance + 1.0)) - 1.0) * clip.w;
    return clip;
}

void main()
    {


    lpos=(modelView*vec4(star_pos,1.0)).xyz;
    ppos=(modelView*vec4(planet_pos,1.0)).xyz;

    lightDir =normalize(lpos-ppos);

    pixel_pos=vec3(modelView*vec4(aPos,1.0));
    pixel_nor=mat3(transpose(inverse(modelView))) * aNormal;
    pixel_txy=aTexCoord;

    vec3 viewDir = normalize(-pixel_pos);


    angleIncidence = acos(dot(lightDir, pixel_nor)) / PI;

    float shadeFactor = 0.1 * (1 - angleIncidence) + 0.9 * (1 - (clamp(angleIncidence, 0.5, 0.5 + transitionWidth) - 0.5) / transitionWidth);

    float angleToViewer = sin(acos(dot(pixel_nor, viewDir)));

    float perspectiveFactor = 0.3 + 0.2 * pow(angleToViewer, fresnelExponent) + 0.5 * pow(angleToViewer, fresnelExponent * 20);

    colAtmosphere = vec4(perspectiveFactor*shadeFactor);

    gl_Position = modelToClipCoordinates(vec4(aPos, 1.0), projection * modelView, u_logarithmicDepthConstant, og_farPlaneDistance);
    }

Fragment Shader

#version 330 core
out vec4 FragColor;

// Fragment
uniform vec3 star_pos;      // point light source position (Sun) [GCS]
uniform vec3 planet_pos;    // planet center position [GCS]
//uniform vec3 planet_r;        // planet radius
uniform float overglow;     // cos(angle) of terminator propagation

uniform sampler2D diffuseTex;
uniform sampler2D txratm;
in vec3 lpos;           // point light source position (Sun) [camera space]
in vec3 ppos;           // planet center position [camera space]
in vec3 pixel_pos;      // fragment position [camera space]
in vec3 pixel_nor;      // fragment surface normal
in vec2 pixel_txy;      // fragment surface texture coord
in float angleIncidence;
in vec4 colAtmosphere;
in vec3 lightDir;

void main()
    {
    float li;
    vec3 c,lt_dir,c0;
    vec4 c1;
    lt_dir=lightDir; // vector from fragment to point light source
    li=dot(pixel_nor,lt_dir)+overglow;
    if (li<0.0) {
        li=0.0;
    }
    if (li>1.0) {
        li=1.0;
    }

    vec2 gradientLevel = vec2(angleIncidence, 0.5);
    c1 = colAtmosphere * texture(txratm, gradientLevel) * 1.4;

    c0=texture(diffuseTex,pixel_txy).rgb * li;

    c = c1.a * c1.rgb + (vec3(1.0, 1.0, 1.0) - c1.a) * c0;;


    FragColor=vec4(c,1.0);
//  gl_FragDepth=0.0;
    }

Solution

I have solved my issue. The problem was in the vertex shader's acos functions. For some reason, the dot product between the normal and the light and view directions was giving a value higher than one(Which I am not sure why it happened, given that all of these values should be normalized).

This was simply fixed by checking the value returned by the dot product, and selecting the minimum between it and the value of 1.

Here is the updated vertex shader code:

#version 330 core
layout (location = 0) in vec3 aPos;
layout (location = 1) in vec2 aTexCoord;
layout (location = 2) in vec3 aNormal;

// Vertex
uniform vec3 star_pos;      // point light source position (Sun) [GCS]
uniform vec3 planet_pos;    // planet center position [GCS]
uniform float overglow;     // cos(angle) of terminator propagation\

uniform mat4 modelView;
uniform mat4 projection;
uniform float og_farPlaneDistance;
uniform float u_logarithmicDepthConstant;

out vec3 lpos;          // point light source position (Sun) [camera space]
out vec3 ppos;          // planet center position [camera space]
out vec3 pixel_pos;     // fragment position [camera space]
out vec3 pixel_nor;     // fragment surface normal
out vec2 pixel_txy;     // fragment surface texture coord
out float angleIncidence;
out vec4 colAtmosphere; //color of the atmosphere
out vec3 lightDir; //direction of light in camera space

const float PI = 3.14159265;
const float transitionWidth = 0.1; //How prominent the atmosphere is
const float fresnelExponent = 20;


vec4 modelToClipCoordinates(vec4 position, mat4 modelViewPerspectiveMatrix, float depthConstant, float farPlaneDistance){
    vec4 clip = modelViewPerspectiveMatrix * position;

    clip.z = ((2.0 * log(depthConstant * clip.z + 1.0) / log(depthConstant * farPlaneDistance + 1.0)) - 1.0) * clip.w;
    return clip;
}

void main()
    {


    lpos=(modelView*vec4(star_pos,1.0)).xyz;
    ppos=(modelView*vec4(planet_pos,1.0)).xyz;

    lightDir =normalize(lpos-ppos);

    pixel_pos=vec3(modelView*vec4(aPos,1.0));
    pixel_nor=normalize(mat3(transpose(inverse(modelView))) * aNormal);
    pixel_txy=aTexCoord;

    vec3 viewDir = normalize(-pixel_pos);

    float dotProd = dot(lightDir, pixel_nor);
    dotProd = min(dotProd,1.0);

    angleIncidence = acos(dotProd) / PI;

    float shadeFactor = 0.1 * (1 - angleIncidence) + 0.9 * (1 - (clamp(angleIncidence, 0.5, 0.5 + transitionWidth) - 0.5) / transitionWidth);

    float dotProd2 = dot(pixel_nor, viewDir);
    dotProd2 = min(dotProd2,1.0);

    float angleToViewer = sin(acos(dotProd2));

    float perspectiveFactor = 0.3 + 0.2 * pow(angleToViewer, fresnelExponent) + 0.5 * pow(angleToViewer, fresnelExponent * 20);

    colAtmosphere = vec4(perspectiveFactor*shadeFactor);

    gl_Position = modelToClipCoordinates(vec4(aPos, 1.0), projection * modelView, u_logarithmicDepthConstant, og_farPlaneDistance);
    }