Atmospheric light scattering implementation

I'm trying to implement an atmospheric scattering in OpenGL. I'm using this "paper" as tutorial.

However I have some difficulties to understand certain points and to figure out some constants.

Basically I've to implement these formulas:

Firstly i don't know if s is the distance from eye to the dome or the distance from eye to the light source (here sun) position. Same for the angle theta I can't figure out if it's the angle from ground to sun or to the dome position the eye is looking at.

Secondly in this slide:

...it tells me the blue color of the sky will appears. I know it's cause of Rayleigh scattering but there is something i can't understand. All the calculation in the formulas above give me a scalar: so how a white light of the sun wich is basically a vec3(1,1,1), will become blue when I multiply it by scalars, it will only get in gray scale because I will have for result for example vec3(0.8,0.8,0.8). I mean , if some different sky color appears, I must multiply the sun light with a vec3 to change the RGB value differently.

Now I encountered some difficulties to implement my shader. Here is the code for the sky shader:

#version 330
in vec3 vpoint;
in vec2 vtexcoord;

out vec2 uv;
out vec3 atmos;

uniform mat4 M;
uniform mat4 V;
uniform mat4 P;

mat4 MVP = P*V*M;

//uniform vec3 lpos;
vec3 lpos = vec3(100,0,0);

uniform vec3 cpos;

vec3 br = vec3(5.5e-6, 13.0e-6, 22.4e-6);
vec3 bm = vec3(21e-6);

float g = -0.75f;

vec3 Esun = vec3(2000,2000,2000);

vec3 Br(float theta) {
    return 3/(16*3.14) * br * (1+cos(theta)*cos(theta));
}

vec3 Bm(float theta) {
    return 1/(4*3.14) * bm * ((1 - g)*(1 - g))/(pow(1+g*g-2*g*cos(theta),3/2));
}

vec3 atmospheric(float theta, float s) {
    return (Br(theta)*Bm(theta))/(br+bm) * Esun * (1- exp( -(br+bm)*s ));
}

void main() {
    gl_Position = MVP * vec4(vpoint, 1.0);
    uv = vtexcoord;

    vec3 domePos = vec3(M*vec4(vpoint,1.0));

    vec3 ldir = lpos - domePos;

    float s = length(domePos-cpos);
    float theta = acos(dot(normalize(ldir-domePos),normalize(domePos-
    cpos)*vec3(1,1,0)));

    atmos = atmospheric(theta,s)*1000000*5;
}

I don't get what I'm expected, here is what I get:

I only have the blue, and no reddish sunset, yet the sun is low and according to the different tutorials I have seen, i should see some reddish color appear when the sun goes low.

Warning I'm not an expert on this field, take this with a grain of salt.

This pretty much says it all.

s is the distance between the vertex/pixel and the camera.
θ is the angle between the sun and the line of sight.

In order to compute θ you need to know the "yellow line" and the "line of sight".
The latter is ordinary shader math; the former is just a way to express how high on the sky the sun is. You can model it as a ray from the sun to a point on the ground.

All the formula above gives you vectors.
L₀ is a vector.
E_sun is also a vector.

The slides basically says that the physic concepts like Radiance and Irradiance (E_sun) are continuous on the spectrum and one should use a Spectral Power Distribution to describe lights and colors.
A faster approach however is to do the math only on three points of the spectrum, the one for the R, G and B wavelengths.
In practice this says that E_sun is a vector describing the irradiance of the sun for the three RGB wavelength.

The blue of the sky comes from the parameter β_R which depends on θ which depends on the "line of sight" which depends on the altitude of the fragment of the sky being coloured.