src/world/world.cpp

   1 #include "Body.hpp"
   2 #include "Orbit.hpp"
   3 #include "Planet.hpp"
   4 #include "Simulation.hpp"
   5 #include "Sun.hpp"
   6 #include "Tile.hpp"
   7
   8 #include "../const.hpp"
   9 #include "../app/Assets.hpp"
  10 #include "../graphics/Viewport.hpp"
  11 #include "../rand/OctaveNoise.hpp"
  12 #include "../rand/SimplexNoise.hpp"
  13
  14 #include <algorithm>
  15 #include <cmath>
  16 #include <iostream>
  17 #include <glm/gtc/matrix_transform.hpp>
  18 #include <glm/gtx/euler_angles.hpp>
  19 #include <glm/gtx/io.hpp>
  20 #include <glm/gtx/transform.hpp>
  21
  22 using blobs::G;
  23 using blobs::PI_2p0;
  24
  25 using std::sin;
  26 using std::cos;
  27 using std::pow;
  28 using std::sqrt;
  29
  30
  31 namespace blobs {
  32 namespace world {
  33
  34 Body::Body()
  35 : sim(nullptr)
  36 , parent(nullptr)
  37 , children()
  38 , mass(1.0)
  39 , radius(1.0)
  40 , orbit()
  41 , surface_tilt(0.0, 0.0)
  42 , axis_tilt(0.0, 0.0)
  43 , rotation(0.0)
  44 , angular(0.0)
  45 , orbital(1.0)
  46 , inverse_orbital(1.0)
  47 , local(1.0)
  48 , inverse_local(1.0) {
  49 }
  50
  51 Body::~Body() {
  52 }
  53
  54 void Body::SetSimulation(Simulation &s) noexcept {
  55         sim = &s;
  56         for (auto child : children) {
  57                 child->SetSimulation(s);
  58         }
  59 }
  60
  61 void Body::SetParent(Body &p) {
  62         if (HasParent()) {
  63                 UnsetParent();
  64         }
  65         parent = &p;
  66         parent->AddChild(*this);
  67 }
  68
  69 void Body::UnsetParent() {
  70         if (!HasParent()) return;
  71         parent->RemoveChild(*this);
  72         parent = nullptr;
  73 }
  74
  75 void Body::AddChild(Body &c) {
  76         children.push_back(&c);
  77         c.SetSimulation(*sim);
  78 }
  79
  80 void Body::RemoveChild(Body &c) {
  81         auto entry = std::find(children.begin(), children.end(), &c);
  82         if (entry != children.end()) {
  83                 children.erase(entry);
  84         }
  85 }
  86
  87 double Body::Inertia() const noexcept {
  88         // assume solid sphere for now
  89         return (2.0/5.0) * Mass() * pow(Radius(), 2);
  90 }
  91
  92 double Body::GravitationalParameter() const noexcept {
  93         return G * Mass();
  94 }
  95
  96 double Body::OrbitalPeriod() const noexcept {
  97         if (parent) {
  98                 return PI_2p0 * sqrt(pow(orbit.SemiMajorAxis(), 3) / (G * (parent->Mass() + Mass())));
  99         } else {
 100                 return 0.0;
 101         }
 102 }
 103
 104 double Body::RotationalPeriod() const noexcept {
 105         if (std::abs(angular) < std::numeric_limits<double>::epsilon()) {
 106                 return std::numeric_limits<double>::infinity();
 107         } else {
 108                 return PI_2p0 * Inertia() / angular;
 109         }
 110 }
 111
 112 glm::dmat4 Body::ToUniverse() const noexcept {
 113         glm::dmat4 m(1.0);
 114         const Body *b = this;
 115         while (b->HasParent()) {
 116                 m = b->ToParent() * m;
 117                 b = &b->Parent();
 118         }
 119         return m;
 120 }
 121
 122 glm::dmat4 Body::FromUniverse() const noexcept {
 123         glm::dmat4 m(1.0);
 124         const Body *b = this;
 125         while (b->HasParent()) {
 126                 m *= b->FromParent();
 127                 b = &b->Parent();
 128         }
 129         return m;
 130 }
 131
 132 void Body::Cache() noexcept {
 133         if (parent) {
 134                 orbital =
 135                         orbit.Matrix(PI_2p0 * (GetSimulation().Time() / OrbitalPeriod()))
 136                         * glm::eulerAngleXY(axis_tilt.x, axis_tilt.y);
 137                 inverse_orbital =
 138                         glm::eulerAngleYX(-axis_tilt.y, -axis_tilt.x)
 139                         * orbit.InverseMatrix(PI_2p0 * (GetSimulation().Time() / OrbitalPeriod()));
 140         } else {
 141                 orbital = glm::eulerAngleXY(axis_tilt.x, axis_tilt.y);
 142                 inverse_orbital = glm::eulerAngleYX(-axis_tilt.y, -axis_tilt.x);
 143         }
 144         local =
 145                 glm::eulerAngleY(rotation)
 146                 * glm::eulerAngleXY(surface_tilt.x, surface_tilt.y);
 147         inverse_local =
 148                 glm::eulerAngleYX(-surface_tilt.y, -surface_tilt.x)
 149                 * glm::eulerAngleY(-rotation);
 150 }
 151
 152
 153 Orbit::Orbit()
 154 : sma(1.0)
 155 , ecc(0.0)
 156 , inc(0.0)
 157 , asc(0.0)
 158 , arg(0.0)
 159 , mna(0.0) {
 160 }
 161
 162 Orbit::~Orbit() {
 163 }
 164
 165 double Orbit::SemiMajorAxis() const noexcept {
 166         return sma;
 167 }
 168
 169 Orbit &Orbit::SemiMajorAxis(double s) noexcept {
 170         sma = s;
 171         return *this;
 172 }
 173
 174 double Orbit::Eccentricity() const noexcept {
 175         return ecc;
 176 }
 177
 178 Orbit &Orbit::Eccentricity(double e) noexcept {
 179         ecc = e;
 180         return *this;
 181 }
 182
 183 double Orbit::Inclination() const noexcept {
 184         return inc;
 185 }
 186
 187 Orbit &Orbit::Inclination(double i) noexcept {
 188         inc = i;
 189         return *this;
 190 }
 191
 192 double Orbit::LongitudeAscending() const noexcept {
 193         return asc;
 194 }
 195
 196 Orbit &Orbit::LongitudeAscending(double l) noexcept {
 197         asc = l;
 198         return *this;
 199 }
 200
 201 double Orbit::ArgumentPeriapsis() const noexcept {
 202         return arg;
 203 }
 204
 205 Orbit &Orbit::ArgumentPeriapsis(double a) noexcept {
 206         arg = a;
 207         return *this;
 208 }
 209
 210 double Orbit::MeanAnomaly() const noexcept {
 211         return mna;
 212 }
 213
 214 Orbit &Orbit::MeanAnomaly(double m) noexcept {
 215         mna = m;
 216         return *this;
 217 }
 218
 219 namespace {
 220
 221 double mean2eccentric(double M, double e) {
 222         double E = M; // eccentric anomaly, solve M = E - e sin E
 223         // limit to 100 steps to prevent deadlocks in impossible situations
 224         for (int i = 0; i < 100; ++i) {
 225                 double dE = (E - e * sin(E) - M) / (1 - e * cos(E));
 226                 E -= dE;
 227                 if (abs(dE) < 1.0e-6) break;
 228         }
 229         return E;
 230 }
 231
 232 }
 233
 234 glm::dmat4 Orbit::Matrix(double t) const noexcept {
 235         double M = mna + t;
 236         double E = mean2eccentric(M, ecc);
 237
 238         // coordinates in orbital plane, P=x, Q=-z
 239         double P = sma * (cos(E) - ecc);
 240         double Q = sma * sin(E) * sqrt(1 - (ecc * ecc));
 241
 242         return glm::yawPitchRoll(asc, inc, arg) * glm::translate(glm::dvec3(P, 0.0, -Q));
 243 }
 244
 245 glm::dmat4 Orbit::InverseMatrix(double t) const noexcept {
 246         double M = mna + t;
 247         double E = mean2eccentric(M, ecc);
 248         double P = sma * (cos(E) - ecc);
 249         double Q = sma * sin(E) * sqrt(1 - (ecc * ecc));
 250         return glm::translate(glm::dvec3(-P, 0.0, Q)) * glm::transpose(glm::yawPitchRoll(asc, inc, arg));
 251 }
 252
 253
 254 Planet::Planet(int sidelength)
 255 : Body()
 256 , sidelength(sidelength)
 257 , tiles(new Tile[TilesTotal()])
 258 , vao() {
 259         Radius(double(sidelength) / 2.0);
 260 }
 261
 262 Planet::~Planet() {
 263 }
 264
 265 glm::dvec3 Planet::TileCenter(int surface, int x, int y) const noexcept {
 266         glm::dvec3 center(0.0f);
 267         center[(surface + 0) % 3] = x + 0.5 - Radius();
 268         center[(surface + 1) % 3] = y + 0.5 - Radius();
 269         center[(surface + 2) % 3] = surface < 3 ? Radius() : -Radius();
 270         return center;
 271 }
 272
 273 void Planet::BuildVAOs() {
 274         vao.Bind();
 275         vao.BindAttributes();
 276         vao.EnableAttribute(0);
 277         vao.EnableAttribute(1);
 278         vao.AttributePointer<glm::vec3>(0, false, offsetof(Attributes, position));
 279         vao.AttributePointer<glm::vec3>(1, false, offsetof(Attributes, tex_coord));
 280         vao.ReserveAttributes(TilesTotal() * 4, GL_STATIC_DRAW);
 281         {
 282                 auto attrib = vao.MapAttributes(GL_WRITE_ONLY);
 283                 float offset = Radius();
 284
 285                 // srf  0  1  2  3  4  5
 286                 //  up +Z +X +Y -Z -X -Y
 287
 288                 for (int index = 0, surface = 0; surface < 6; ++surface) {
 289                         for (int y = 0; y < sidelength; ++y) {
 290                                 for (int x = 0; x < sidelength; ++x, ++index) {
 291                                         float tex = TileAt(surface, x, y).type;
 292                                         const float tex_u_begin = surface < 3 ? 1.0f : 0.0f;
 293                                         const float tex_u_end = surface < 3 ? 0.0f : 1.0f;
 294                                         attrib[4 * index + 0].position[(surface + 0) % 3] = x + 0 - offset;
 295                                         attrib[4 * index + 0].position[(surface + 1) % 3] = y + 0 - offset;
 296                                         attrib[4 * index + 0].position[(surface + 2) % 3] = surface < 3 ? offset : -offset;
 297                                         attrib[4 * index + 0].tex_coord[0] = tex_u_begin;
 298                                         attrib[4 * index + 0].tex_coord[1] = 1.0f;
 299                                         attrib[4 * index + 0].tex_coord[2] = tex;
 300
 301                                         attrib[4 * index + 1].position[(surface + 0) % 3] = x + 0 - offset;
 302                                         attrib[4 * index + 1].position[(surface + 1) % 3] = y + 1 - offset;
 303                                         attrib[4 * index + 1].position[(surface + 2) % 3] = surface < 3 ? offset : -offset;
 304                                         attrib[4 * index + 1].tex_coord[0] = tex_u_end;
 305                                         attrib[4 * index + 1].tex_coord[1] = 1.0f;
 306                                         attrib[4 * index + 1].tex_coord[2] = tex;
 307
 308                                         attrib[4 * index + 2].position[(surface + 0) % 3] = x + 1 - offset;
 309                                         attrib[4 * index + 2].position[(surface + 1) % 3] = y + 0 - offset;
 310                                         attrib[4 * index + 2].position[(surface + 2) % 3] = surface < 3 ? offset : -offset;
 311                                         attrib[4 * index + 2].tex_coord[0] = tex_u_begin;
 312                                         attrib[4 * index + 2].tex_coord[1] = 0.0f;
 313                                         attrib[4 * index + 2].tex_coord[2] = tex;
 314
 315                                         attrib[4 * index + 3].position[(surface + 0) % 3] = x + 1 - offset;
 316                                         attrib[4 * index + 3].position[(surface + 1) % 3] = y + 1 - offset;
 317                                         attrib[4 * index + 3].position[(surface + 2) % 3] = surface < 3 ? offset : -offset;
 318                                         attrib[4 * index + 3].tex_coord[0] = tex_u_end;
 319                                         attrib[4 * index + 3].tex_coord[1] = 0.0f;
 320                                         attrib[4 * index + 3].tex_coord[2] = tex;
 321                                 }
 322                         }
 323                 }
 324         }
 325         vao.BindElements();
 326         vao.ReserveElements(TilesTotal() * 6, GL_STATIC_DRAW);
 327         {
 328                 auto element = vao.MapElements(GL_WRITE_ONLY);
 329                 int index = 0;
 330                 for (int surface = 0; surface < 3; ++surface) {
 331                         for (int y = 0; y < sidelength; ++y) {
 332                                 for (int x = 0; x < sidelength; ++x, ++index) {
 333                                         element[6 * index + 0] = 4 * index + 0;
 334                                         element[6 * index + 1] = 4 * index + 2;
 335                                         element[6 * index + 2] = 4 * index + 1;
 336                                         element[6 * index + 3] = 4 * index + 1;
 337                                         element[6 * index + 4] = 4 * index + 2;
 338                                         element[6 * index + 5] = 4 * index + 3;
 339                                 }
 340                         }
 341                 }
 342                 for (int surface = 3; surface < 6; ++surface) {
 343                         for (int y = 0; y < sidelength; ++y) {
 344                                 for (int x = 0; x < sidelength; ++x, ++index) {
 345                                         element[6 * index + 0] = 4 * index + 0;
 346                                         element[6 * index + 1] = 4 * index + 1;
 347                                         element[6 * index + 2] = 4 * index + 2;
 348                                         element[6 * index + 3] = 4 * index + 2;
 349                                         element[6 * index + 4] = 4 * index + 1;
 350                                         element[6 * index + 5] = 4 * index + 3;
 351                                 }
 352                         }
 353                 }
 354         }
 355         vao.Unbind();
 356 }
 357
 358 void Planet::Draw(app::Assets &assets, graphics::Viewport &viewport) {
 359         vao.Bind();
 360         const glm::mat4 &MV = assets.shaders.planet_surface.MV();
 361         assets.shaders.planet_surface.SetNormal(glm::vec3(MV * glm::vec4(0.0f, 0.0f, 1.0f, 0.0f)));
 362         vao.DrawTriangles(TilesPerSurface() * 6, TilesPerSurface() * 6 * 0);
 363         assets.shaders.planet_surface.SetNormal(glm::vec3(MV * glm::vec4(1.0f, 0.0f, 0.0f, 0.0f)));
 364         vao.DrawTriangles(TilesPerSurface() * 6, TilesPerSurface() * 6 * 1);
 365         assets.shaders.planet_surface.SetNormal(glm::vec3(MV * glm::vec4(0.0f, 1.0f, 0.0f, 0.0f)));
 366         vao.DrawTriangles(TilesPerSurface() * 6, TilesPerSurface() * 6 * 2);
 367         assets.shaders.planet_surface.SetNormal(glm::vec3(MV * glm::vec4(0.0f, 0.0f, -1.0f, 0.0f)));
 368         vao.DrawTriangles(TilesPerSurface() * 6, TilesPerSurface() * 6 * 3);
 369         assets.shaders.planet_surface.SetNormal(glm::vec3(MV * glm::vec4(-1.0f, 0.0f, 0.0f, 0.0f)));
 370         vao.DrawTriangles(TilesPerSurface() * 6, TilesPerSurface() * 6 * 4);
 371         assets.shaders.planet_surface.SetNormal(glm::vec3(MV * glm::vec4(0.0f, -1.0f, 0.0f, 0.0f)));
 372         vao.DrawTriangles(TilesPerSurface() * 6, TilesPerSurface() * 6 * 5);
 373 }
 374
 375
 376 void GenerateEarthlike(Planet &p) noexcept {
 377         rand::SimplexNoise elevation_gen(0);
 378
 379         constexpr int ice = 0;
 380         constexpr int grass = 3;
 381         constexpr int water = 4;
 382         constexpr int sand = 5;
 383         constexpr int rock = 8;
 384
 385         constexpr double water_thresh = 0.0;
 386         constexpr double beach_thresh = 0.1;
 387         constexpr double mountain_thresh = 0.5;
 388
 389         const glm::dvec3 axis(glm::dvec4(0.0, 1.0, 0.0, 0.0) * glm::eulerAngleXY(p.SurfaceTilt().x, p.SurfaceTilt().y));
 390         const double cap_thresh = std::cos(p.AxialTilt().x);
 391
 392         for (int surface = 0; surface <= 5; ++surface) {
 393                 for (int y = 0; y < p.SideLength(); ++y) {
 394                         for (int x = 0; x < p.SideLength(); ++x) {
 395                                 glm::dvec3 to_tile = p.TileCenter(surface, x, y);
 396                                 double near_axis = std::abs(glm::dot(glm::normalize(to_tile), axis));
 397                                 if (near_axis > cap_thresh) {
 398                                         p.TileAt(surface, x, y).type = ice;
 399                                         continue;
 400                                 }
 401                                 float elevation = rand::OctaveNoise(
 402                                         elevation_gen,
 403                                         to_tile / p.Radius(),
 404                                         3,   // octaves
 405                                         0.5, // persistence
 406                                         2 / p.Radius(), // frequency
 407                                         2,   // amplitude
 408                                         2    // growth
 409                                 );
 410                                 if (elevation < water_thresh) {
 411                                         p.TileAt(surface, x, y).type = water;
 412                                 } else if (elevation < beach_thresh) {
 413                                         p.TileAt(surface, x, y).type = sand;
 414                                 } else if (elevation < mountain_thresh) {
 415                                         p.TileAt(surface, x, y).type = grass;
 416                                 } else {
 417                                         p.TileAt(surface, x, y).type = rock;
 418                                 }
 419                         }
 420                 }
 421         }
 422         p.BuildVAOs();
 423 }
 424
 425 void GenerateTest(Planet &p) noexcept {
 426         for (int surface = 0; surface <= 5; ++surface) {
 427                 for (int y = 0; y < p.SideLength(); ++y) {
 428                         for (int x = 0; x < p.SideLength(); ++x) {
 429                                 if (x == p.SideLength() / 2 && y == p.SideLength() / 2) {
 430                                         p.TileAt(surface, x, y).type = surface;
 431                                 } else {
 432                                         p.TileAt(surface, x, y).type = (x == p.SideLength()/2) + (y == p.SideLength()/2) + 6;
 433                                 }
 434                         }
 435                 }
 436         }
 437         p.BuildVAOs();
 438 }
 439
 440
 441 Sun::Sun()
 442 : Body() {
 443 }
 444
 445 Sun::~Sun() {
 446 }
 447
 448 }
 449 }