Hi everyone. I’m trying to write a sky shader based on Nishita’s model that allows you to specify altitude as one of the sky parameters. I’m following this article here: http://www.scratchapixel.com/lessons/3d-advanced-lessons/simulating-the-colors-of-the-sky/atmospheric-scattering/

So far I’ve had some luck, and my shader compiles, but I can’t seem to set the sun direction properly - whatever vector I specify the sun always seems to end up below or at the horizon. I’m clearly doing something wrong with some vector maths somewhere, but I can’t track it down. Can anyone see what I’m doing wrong?

Current shader code:

```
#include <stdosl.h>
shader nishita_sky (
float Height = 1,
float Intensity = 20,
vector InSunDirection = vector(0,0,1),
vector InDirection = P,
output color Rayleigh = 0.0,
output color Mie = 0.0,
output color Sky = 0.0)
{
/* Set out my important variables */
float Re = 6360e3; // Radius of the earth
float Ra = 6420e3; // Radius of the atmosphere
float Hr = 7994; // Rayleigh scale height
float Hm = 1200; // Mie scale height
float g = 0.76; // Anisotropy term for Mie scattering.
vector BetaR = vector(5.8e-6,13.0e-6,22.4e-6);
vector BetaM = vector(21e-5);
vector SumR = vector(0);
vector SumM = vector(0);
/* Compute intersection point of camera ray with atmosphere*/
// Five possible situations:
// 1. No intersection, return black.
// 2. Camera inside planet. Return black.
// 3. Camera inside atmosphere. Return correct d value.
// 4. Camera inside atmoshere, facing earth, return d for collision with earth.
// 5. Camera outside atmosphere, can see through to other side. Return two d values one for each entry and exit.
// 6. Camera outside atmosphere, can't see through to other side. Return two d values, one for entry, one for earth.
// For this we need d, the distance along the camera ray we need to travel to reach the edge of the atmosphere.
vector Direction = normalize (InDirection);
vector SunDirection = normalize (InSunDirection);
// Our start and end points for our integral, these should be overwritten below.
vector Pu = vector(0,0,0);
vector Pv = vector(0,0,0);
if (Height > 0) { // Checks height is above ground.
vector Pc = vector (0, 0, Re + Height);
Pu = Pc;
float b = 2 * dot (Direction, Pc);
float bsquared = pow(b,2);
float fourac = 4 * dot(Direction, Direction) * (1 - pow(Ra,2));
if (bsquared - fourac > 0) { //Checks that our view direction intersects with the atmosphere.
float d1 = (-b - pow(bsquared - fourac, 0.5)) / 2;
float d2 = (-b + pow(bsquared - fourac, 0.5)) / 2;
// Check solutions for our outcomes:
if (Pc[2] > Ra) { //Checks if camera is outside atmosphere.
//We've already ruled out no solutions, so we expect two, they should both be positive.
if (d1 > 0 && d2 > 0) {
float d_entry = min(d1,d2);
float d_exit = max(d1,d2);
vector Pe = Pc + d_entry * Direction; // Pe is the entry point to the atmosphere.
Pu = Pe; // Overwrite Pu with Pe, its now our start point for the integral.
// Does the view pass through earth? If so find Pr.
float b_earth = 2 * dot (Direction, Pc);
float bsquared_earth = 2 * b_earth;
float fourac_earth = 4 * dot (Direction,Direction) * (1 - pow(Re, 2));
if (bsquared_earth - fourac_earth > 0) { // Check if ray hits earth.
float d1_earth = (-b_earth - pow(bsquared_earth - fourac_earth, 0.5)) / 2;
float d2_earth = (-b_earth + pow(bsquared_earth - fourac_earth, 0.5)) / 2;
if (d1_earth > 0 && d2_earth > 0) { // Solutions aren't negative.
float d_entry = min(d1_earth,d2_earth);
vector Pr = Pc + d_entry * Direction;
Pv = Pr;
}
} else { // View doesn't hit earth, exits atmosphere.
vector Pa = Pc + d_exit * Direction;
Pv = Pa;
}
}
} else {
// Our camera is inside the atmosphere.
// We expect one negative solution, one positive, we only want the positive.
float d_exit = max(d1,d2);
vector Pa = Pc + d_exit * Direction;
Pv = Pa;
}
}
}
/* Now compute Rayleigh and Mie Phases, these are constant terms, not part of the integral.*/
float mu = dot(Direction, SunDirection);
float phaseR = (3/(16 * M_PI)) * (1 + mu*mu);
float phaseM = (3/(8 * M_PI)) * ((1 - g*g) * (1 + mu*mu) / ((2 + g*g) * pow(1 + g*g - 2 * g * mu, 1.5)));
/* Create samples along camera ray. */
int numSamples = 16;
int numSamplesL = 8;
float segmentLength = distance (Pu, Pv) / numSamples;
float opticalDepthR = 0;
float opticalDepthM = 0;
// Now begin taking samples for the integtral.
for (int i = 0; i < numSamples; i++) {
vector samplePosition = Pu + segmentLength * Direction * (i + 0.5);
float sampleHeight = length(samplePosition) - Re;
/* Get optical depth for light at this sample: */
float Hr_sample = exp(-sampleHeight/Hr) * segmentLength;
float Hm_sample = exp(-sampleHeight/Hm) * segmentLength;
opticalDepthR += Hr_sample;
opticalDepthM += Hm_sample;
/* Find light sample termination point Ps for this sample: */
// There are 2 possibilities here (we now know we are in the atmosphere because only lines within the atmosphere are sampled):
// 1. The sun is visible from this position.
// 2. It isn't, because its behind the earth.
float opticalDepthLR = 0;
float opticalDepthLM = 0;
// Check that the vector in the direction of the sun from the sample position doesnt intersect the earth.
float b_s = 2 * dot (SunDirection,samplePosition);
float bsquared_s = pow(b_s, 2);
float fourac_s_earth = 4 * dot(SunDirection,SunDirection) * (1 - pow(Re,2));
int canseesun = 0;
if (bsquared_s - fourac_s_earth > 0) { // Check if our vector towards the sun goes through the earth. If so dont bother sampling.
// There are Solutions! But if they're both negative we can carry on.
float d1_s_earth = (-b_s - pow(bsquared_s - fourac_s_earth, 0.5)) / 2;
float d2_s_earth = (-b_s + pow(bsquared_s - fourac_s_earth, 0.5)) / 2;
if (d1_s_earth < 0 && d2_s_earth < 0) {
//Both of our earth intersection vectors are behind the camera. So we can see the sun.
canseesun = 1;
}
} else { //Doesn't hit earth. Can see sun.
canseesun = 1;
}
if (canseesun) {
// Calculate distance to Ps and do samples.
float fourac_s = 4 * dot(SunDirection, SunDirection) * (1 - pow(Ra,2)); // This really shouldn't be negative at any time, that would mean we we inside a sphere but somehow the sun wasn't outside it...
float d1_s = (-b_s - pow(bsquared_s - fourac_s, 0.5)) / 2;
float d2_s = (-b_s + pow(bsquared_s - fourac_s, 0.5)) / 2;
// One of the two above should be negative, pick the positive(negative?) one.
float d_s = 0;
if (d1_s > 0) {
d_s = d1_s;
} else {
d_s = d2_s;
}
vector Ps = samplePosition + d_s * SunDirection;
/* Create samples along sun ray*/
float segmentLengthL = d_s / numSamplesL;
for (int j = 0; j < numSamplesL; j++) {
vector samplePositionL = samplePosition + segmentLengthL * SunDirection * (j + 0.5);
/* Test if light ray is below ground, and discard if it is. */
float sampleHeightL = length(samplePositionL) - Re;
if (sampleHeightL > 0) {
// ignore sampleheights < 0 they're inside the planet...
/* Get optical depth for light at this sample */
float Hlr_sample = exp(-sampleHeightL/Hr)*segmentLengthL;
float Hlm_sample = exp(-sampleHeightL/Hm)*segmentLengthL;
opticalDepthLR += Hlr_sample;
opticalDepthLM += Hlm_sample;
}
}
}
/* once we've done all our samples we can calculate our attenuation for the light */
vector tauR = BetaR * (opticalDepthR + opticalDepthLR);
vector tauM = BetaM * (opticalDepthM + opticalDepthLM);
vector attnR = vector (exp(-tauR[0]), exp(-tauR[1]), exp(-tauR[2]));
vector attnM = vector (exp(-tauM));
SumR += Hr_sample * attnR;
SumM += Hm_sample * attnM;
}
//printf("%.2f
", Direction);
Rayleigh = color (SumR * phaseR * BetaR) * Intensity;
Mie = color (SumM * phaseM * BetaM) * Intensity;
Sky = Rayleigh + Mie;
}
```

Here’s what the output looks like currently (you’re looking at the sky thorough a fisheye lens here):