Instantly share code, notes, and snippets.

Embed
What would you like to do?
GLSL Astronomy compute shader
// 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);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment