GLSL Astronomy compute shader

astronomy.glsl

// Copyright (C) 2018, Benjamin "BeRo" Rosseaux (benjamin@rosseaux.de) - License: CC0
// Hint: It's not super accurate, but it should be good enough for games and demos with sky rendering.
#version 430
layout(local_size_x = 1, local_size_y = 1, local_size_z = 1) in;  

uniform vec3 tc; // Time constant

layout(std430) buffer ssboAstronomy {
   vec3 sunPosition;
   float sunDistance;
   vec3 moonPosition;
   float moonDistance;
   float moonPhase;
   mat3 celestialTransform;
};                                      

const float HalfPI = 1.5707963267948966, // 1.570796326794896619231,
            PI = 3.141592653589793, // 3.141592653589793238463,
            PI2 = 6.283185307179586; // 6.283185307179586476925            

// The AMD windows GPU drivers do not seem to like constant double precision values here, 
// so no const prefix in this case.
double astronomicalUnitsToKiloMeters = 149597871.0;
                                     
const vec2 sinCosOffsets = vec2(0.0, HalfPI);
  
double calculateJulianDate(int year, int month, int day){
  year -= int(month <= 2);
  return double((2 - (year / 100)) + (year / 400)) +
         floor(double(year) * 365.25) +
         floor(30.6001 * double(month + 1 + (int(month <= 2) * 12))) +
         double(day) +
         1720994.5;
}
        
double getElapsedJulianDays(double julianDate){
  return julianDate - 2451545.0;
}

float julianDateToGreenwichMeanSiderealTime(double julianDate){
	double theta = ((floor(julianDate - 0.5) + 0.5) - 2451545.0) / 36525.0;
	return float(mod((6.697374558 + (theta * (2400.051336 + (theta * 0.000025862)))) + 
                   ((fract(julianDate - 0.5) * 24.0) * 1.002737909), 
                   24.0));
}                  

float localMeanSiderealTimeInRadians(double julianDate,
                                     float longitude){
	return radians((julianDateToGreenwichMeanSiderealTime(julianDate) * 15.0) + longitude);
}

vec3 azimuthAltitudeToDirection(const in vec2 azimuthAltitude){
  vec4 sinCosAzimuthAltitude = sin(azimuthAltitude.xxyy + sinCosOffsets.xyxy);  
  return vec3((sinCosAzimuthAltitude.xy * sinCosAzimuthAltitude.w) * vec2(-1.0, 1.0),
              sinCosAzimuthAltitude.z).xzy;
}

vec2 convertEclipticToEquatorialRad(double julianDate,
                                    vec2 longitudeLatitude){
  double elapsedJulianDays = getElapsedJulianDays(julianDate),
         elapsedJulianCenturies = elapsedJulianDays / 36525.0,
         squaredElapsedJulianCenturies = elapsedJulianCenturies * elapsedJulianCenturies;
  float earthAxialTilt = radians(float(mod((23.0 + (26.0 / 60.0) + (21.448 / 3600.0)) - 
                                           (((46.8150 * elapsedJulianCenturies) + 
                                            (0.00059 * squaredElapsedJulianCenturies) - 
                                            (0.001813 * squaredElapsedJulianCenturies * elapsedJulianCenturies)) / 3600.0), 360.0)));
  vec4 sinCosLongitudeLatitude = sin(longitudeLatitude.xxyy + sinCosOffsets.xyxy); 
  vec2 sinCosEarthAxialTilt = sin(vec2(earthAxialTilt) + sinCosOffsets); 
  vec3 vector = vec3(sinCosLongitudeLatitude.y * sinCosLongitudeLatitude.w,
                     (sinCosEarthAxialTilt.yx * sinCosLongitudeLatitude.x * sinCosLongitudeLatitude.w) +
                     (vec2(-sinCosEarthAxialTilt.x, sinCosEarthAxialTilt.y) * sinCosLongitudeLatitude.z)
                    );                
  return vec2(atan(vector.y, vector.x),           // right ascension
              atan(vector.z, length(vector.xy))); // declination
}

vec2 convertEquatorialToHorizontal(double julianDate,
                                   vec2 longitudeLatitude,
                                   vec2 rightAscensionDeclination){
  vec2 sinCosLatitude = sin(vec2(radians(longitudeLatitude.y)) + sinCosOffsets),
       sinCosHourAngle = sin(vec2(localMeanSiderealTimeInRadians(julianDate, longitudeLatitude.x) - rightAscensionDeclination.x) + sinCosOffsets), 
       sinCosDeclination = sin(vec2(rightAscensionDeclination.y) + sinCosOffsets); 
  return vec2(mod(mod(atan(-sinCosHourAngle.x, (tan(rightAscensionDeclination.y) * sinCosLatitude.y) - (sinCosHourAngle.y * sinCosLatitude.x)), PI2) + PI2, PI2), // azimuth
	            asin((sinCosLatitude.y * sinCosHourAngle.y * sinCosDeclination.y) + (sinCosDeclination.x * sinCosLatitude.x))); // altitude
}

vec2 getSunPosition(double julianDate, 
                    vec2 longitudeLatitude){
  double elapsedJulianDays = getElapsedJulianDays(julianDate);
  float meanAnomaly = radians(float(mod(357.5291 + (0.98560028 * elapsedJulianDays), 360.0))),
        equationOfCenter = radians((1.9148 * sin(meanAnomaly)) + 
                                   (0.02 * sin(2.0 * meanAnomaly)) + 
                                   (0.0003 * sin(3.0 * meanAnomaly))),
        perihelionOfEarth = radians(102.9372),
        eclipticLongitude = meanAnomaly + equationOfCenter + perihelionOfEarth + float(PI);
  return convertEquatorialToHorizontal(julianDate, 
                                       longitudeLatitude, 
                                       convertEclipticToEquatorialRad(julianDate,
                                                                      vec2(float(eclipticLongitude), 0.0)));
}

double getSunDistance(double julianDate){
  double elapsedJulianDays = getElapsedJulianDays(julianDate);
  float meanAnomaly = radians(float(mod(357.5291 + (0.98560028 * elapsedJulianDays), 360.0)));
  return ((1.00014 - (0.01671 * cos(meanAnomaly))) - (0.00014 * cos(2.0 * meanAnomaly))) * astronomicalUnitsToKiloMeters;
}

vec2 getMoonPosition(double julianDate, 
                     vec2 longitudeLatitude){
  double elapsedJulianDays = getElapsedJulianDays(julianDate);
  float eclipticLongitude = radians(float(mod((218.316 + (13.176396 * elapsedJulianDays)), 360.0))),
        meanAnomaly = radians(float(mod((134.963 + (13.064993 * elapsedJulianDays)), 360.0))),
        meanDistance = radians(float(mod((93.272 + (13.229350 * elapsedJulianDays)), 360.0))),
        longitude = eclipticLongitude + radians(sin(meanAnomaly) * 6.289), 
        latitude = radians(sin(meanDistance) * 5.128);
  return convertEquatorialToHorizontal(julianDate, 
                                       longitudeLatitude, 
                                       convertEclipticToEquatorialRad(julianDate,
                                                                      vec2(longitude, latitude)));
}

double getMoonDistance(double julianDate){
  return 385001.0 - (20905.0 * cos(radians(float(mod((134.963 + (13.064993 * getElapsedJulianDays(julianDate))), 360.0)))));
}

float getMoonPhase(double julianDate){
  return float(2.0 - abs(2.0 - (4.0 * abs(fract((julianDate - 2454488.0665) / 29.531026)))));
}

mat3 getCelestialTransform(double julianDate, 
                           vec2 longitudeLatitude){
  vec4 sinCosAxisZX = sin(vec2(-(HalfPI + localMeanSiderealTimeInRadians(julianDate, longitudeLatitude.x)),
                               radians(longitudeLatitude.y) - HalfPI).xxyy + 
                          sinCosOffsets.xyxy);
  return mat3(sinCosAxisZX.y, sinCosAxisZX.x, 0.0, -sinCosAxisZX.x, sinCosAxisZX.y, 0.0, 0.0, 0.0, 1.0) *
         mat3(1.0, 0.0, 0.0, 0.0, sinCosAxisZX.w, sinCosAxisZX.z, 0.0, -sinCosAxisZX.z, sinCosAxisZX.w);
}

void main(){
  const vec2 observerLongitudeLatitude = vec2(51.18539, 6.44172); // Latitude and longitude coordinates for Mönchengladbach, Germany
  double julianDate = calculateJulianDate(2018, 7, 17) + double(tc.y / (3.0 * 256.0)); 
  sunPosition = azimuthAltitudeToDirection(getSunPosition(julianDate, observerLongitudeLatitude)) * (sunDistance = float(getSunDistance(julianDate)));
  moonPosition = azimuthAltitudeToDirection(getMoonPosition(julianDate, observerLongitudeLatitude)) * (moonDistance = float(getMoonDistance(julianDate)));
  moonPhase = getMoonPhase(julianDate);
  celestialTransform = getCelestialTransform(julianDate, observerLongitudeLatitude);
}