Commit 88c581f9 authored by Dario Seyb's avatar Dario Seyb

started on atmosphere precomputations

parent b51ecfc7
......@@ -34,3 +34,4 @@
# Other
/imgui.ini
*.trace
transmittance.png
This diff is collapsed.
This diff is collapsed.
......@@ -19,6 +19,8 @@
#include <glm/glm.hpp>
#include <glm/gtx/transform.hpp>
#include <engine/utils/stb_image_write.h>
bool OrbitalSimulationSystem::startup() {
RESOLVE_DEPENDENCY(m_events);
RESOLVE_DEPENDENCY(m_scene);
......@@ -224,7 +226,7 @@ void OrbitalSimulationSystem::simulateOrbitals(float dt, float totalTime) {
p.unpack<Planet, Transform>(planet, transform);
// dt here is interpreted as a DAY NUMBER. TODO: Conversion/scale
pos = applyKepler(planet, dt);
pos = applyKepler(planet, dt );
if (planet->name != "Moon") {
minimapItems[p.id()].component<Transform>()->position = glm::dvec3(pos.x, pos.z, pos.y -70 );
......@@ -236,8 +238,123 @@ void OrbitalSimulationSystem::simulateOrbitals(float dt, float totalTime) {
}
}
// Mass is in solar masses!
double getAtmosphereDistance(glm::dvec3 o, glm::dvec3 l, double r) {
double a = -glm::dot(l, o);
double d = a * a - glm::dot(o, o) + r * r;
if (d < 0) {
return 0;
}
double dsqrt = glm::sqrt(d);
// We are inside the atmosphere, return the distance to the far intersection
return a + dsqrt;
}
double getPlanetDistance(glm::dvec3 o, glm::dvec3 l, double r) {
double a = -glm::dot(l, o);
double d = a * a - glm::dot(o, o) + r * r;
if (d < 0) {
return -1;
}
double dsqrt = glm::sqrt(d);
// We are outside the planet, return the distance to the near intersection
return a - dsqrt;
}
const int TRANSMITTANCE_INTEGRAL_SAMPLES = 25;
double limit(double r, double mu, double Rg, double RL) {
double dout = -r * mu + sqrt(r * r * (mu * mu - 1.0) + RL * RL);
double delta2 = r * r * (mu * mu - 1.0) + Rg * Rg;
if (delta2 >= 0.0) {
double din = -r * mu - sqrt(delta2);
if (din >= 0.0) {
dout = glm::min(dout, din);
}
}
return dout;
}
double opticalDepth(double referenceHeight, double height, double viewAngle, double Rg, double RL) {
double result = 0.0;
double dx = limit(height, viewAngle, Rg, RL + 1) / double(TRANSMITTANCE_INTEGRAL_SAMPLES);
double xi = 0.0;
double yi = exp(-(height - Rg) / referenceHeight);
for (int i = 1; i <= TRANSMITTANCE_INTEGRAL_SAMPLES; ++i) {
double xj = double(i) * dx;
double yj = exp(-(glm::sqrt(height * height + xj * xj + 2.0 * xj * height * viewAngle) - Rg) / referenceHeight);
result += (yi + yj) / 2.0 * dx;
xi = xj;
yi = yj;
}
return viewAngle < -sqrt(1.0 - (Rg / height) * (Rg / height)) ? 1e9 : result;
}
glm::dvec3 computeTransmittance(glm::dvec3 constBetaRayleigh, double constBetaMie, double planetRadius, double atmosphereRadius, double averageDensityHeight, double viewAngle, double height ) {
viewAngle = glm::acos(viewAngle);
glm::dvec3 pos = { 0, height, 0 };
glm::dvec3 dir = glm::normalize(glm::rotateX(glm::dvec3(0, 1, 0), viewAngle));
glm::dvec3 result = { 0, 0, 0 };
int stepCount = 500;
double atmosphereDistance = getAtmosphereDistance(pos, dir, atmosphereRadius);
double planetDistance = getPlanetDistance(pos, dir, planetRadius);
double totalDistance = (atmosphereDistance < planetDistance || planetDistance < 0) ? atmosphereDistance : planetDistance;
double stepSize = totalDistance / stepCount;
double opticalDepthRayleigh = 0;
double opticalDepthMie = 0;
// Integrate
for (int i = 0; i < stepCount; i++) {
height = glm::length(pos);
double altitude = (height - planetRadius);
opticalDepthRayleigh += glm::exp(-altitude / 8000.0) * stepSize;
opticalDepthMie += glm::exp(-altitude / 1200.0) * stepSize;
pos += dir * stepSize;
}
result = constBetaRayleigh * opticalDepthRayleigh + glm::dvec3(constBetaMie/0.9) * opticalDepthMie;
result = glm::exp(-result);
return result;
}
void precomputeAtmosphereLookupTexture(double molecularNumberDensity, double airIOR, glm::dvec3 rgbWavelengths, double planetRadius, double atmosphereRadius, double averageDensityHeight) {
const double wavelengthIndependentFactor = (8.0 * glm::pow3(M_PI) * glm::pow2(glm::pow3(airIOR) - 1.0));
const glm::dvec3 constBetaRayleigh = wavelengthIndependentFactor / (3 * molecularNumberDensity * glm::pow4(rgbWavelengths));
const double constBetaMie = glm::pow(210.0, -5.0); //wavelengthIndependentFactor / (3 * molecularNumberDensity * glm::pow(10.0, 36.0)); //
const int WIDTH = 512;
const int HEIGHT = 512;
auto transmittance = new GLubyte[WIDTH * HEIGHT * 4];
auto pixel = transmittance;
for (int y = 0; y < HEIGHT; y++) {
double height = (1.0 - double(y) / HEIGHT) * (atmosphereRadius - planetRadius) + planetRadius;
for (int x = 0; x < WIDTH; x++) {
double viewAngle = -0.15 + (double(x + 1) / WIDTH) * (1.0 + 0.15); //M_PI * (1.0 - ); // 98.0/180
auto val = glm::vec4(computeTransmittance(constBetaRayleigh, constBetaMie, planetRadius, atmosphereRadius, averageDensityHeight - planetRadius, viewAngle, height), 1.0);
pixel[0] = (GLubyte)(val.r * 255);
pixel[1] = (GLubyte)(val.g * 255);
pixel[2] = (GLubyte)(val.b * 255);
pixel[3] = (GLubyte)(val.a * 255);
pixel += 4;
}
}
stbi_write_png("transmittance.png", WIDTH, HEIGHT, 4, transmittance, 0);
}
// Mass is in solar masses!
Entity OrbitalSimulationSystem::addPlanet(Transform::Handle sun, std::string n, double m, double r, double e, double a, double i, double N, double P, double T) {
// Mass, Radius, Eccentricity, Semimajor axis, Inclination, Ascending Node, Arg. of Periapsis, time at perihelion
......@@ -261,6 +378,10 @@ Entity OrbitalSimulationSystem::addPlanet(Transform::Handle sun, std::string n,
planetEntity.assign<Drawable>(defaultGeom, earthMat, 2, m_renderer->getRenderPassId("Main"_sh));
// Add an atmosphere to the planet
/*precomputeAtmosphereLookupTexture(2.55 * glm::pow(10.0, -11.0 ), // premultiply by nm^4 => 25 - 4*9
1.000206102, glm::dvec3{ 680, 550, 440 }, // in nm
6360000, 6420000, planetTransform->scale.x*1.02);*/
auto atmosphere = m_scene->create();
auto atmosphereTransform = atmosphere.assign<Transform>();
atmosphereTransform->parent = planetTransform;
......
......@@ -393,8 +393,7 @@ void AtmosphereTestScene::loadMainSceneResources() {
.attributeLocations(vaoSkybox->getAttributeLocations())
.fragmentDataLocations(m_renderer->getGBufferLocations()));
auto skyboxTexture =
Texture2DFileManager::the()->get(Texture2DCreator("skybox.png"));
auto skyboxTexture = Texture2DFileManager::the()->get(Texture2DCreator("skybox.png"));
skyboxGeometry = {vaoSkybox};
skyboxGeometry.vao->setMode(GL_PATCHES);
......
#define STB_IMAGE_IMPLEMENTATION
#include <engine/utils/stb_image.h>
\ No newline at end of file
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include <engine/utils/stb_image_write.h>
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment