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