2 #include "CreatureCreatureCollision.hpp"
5 #include "Resource.hpp"
7 #include "Simulation.hpp"
10 #include "TileType.hpp"
12 #include "../app/Assets.hpp"
13 #include "../creature/Composition.hpp"
14 #include "../creature/Creature.hpp"
15 #include "../graphics/Viewport.hpp"
16 #include "../math/const.hpp"
17 #include "../math/geometry.hpp"
18 #include "../math/OctaveNoise.hpp"
19 #include "../math/SimplexNoise.hpp"
24 #include <glm/gtc/matrix_transform.hpp>
25 #include <glm/gtx/euler_angles.hpp>
26 #include <glm/gtx/io.hpp>
27 #include <glm/gtx/transform.hpp>
49 , inverse_orbital(1.0)
59 void Body::SetSimulation(Simulation &s) noexcept {
61 for (auto child : children) {
62 child->SetSimulation(s);
66 void Body::SetParent(Body &p) {
71 parent->AddChild(*this);
74 void Body::UnsetParent() {
75 if (!HasParent()) return;
76 parent->RemoveChild(*this);
80 void Body::AddChild(Body &c) {
81 children.push_back(&c);
82 c.SetSimulation(*sim);
85 void Body::RemoveChild(Body &c) {
86 auto entry = std::find(children.begin(), children.end(), &c);
87 if (entry != children.end()) {
88 children.erase(entry);
92 double Body::Inertia() const noexcept {
93 // assume solid sphere for now
94 return (2.0/5.0) * Mass() * pow(Radius(), 2);
97 double Body::GravitationalParameter() const noexcept {
101 double Body::OrbitalPeriod() const noexcept {
103 return PI * 2.0 * sqrt(pow(orbit.SemiMajorAxis(), 3) / (G * (parent->Mass() + Mass())));
109 double Body::RotationalPeriod() const noexcept {
110 if (std::abs(angular) < std::numeric_limits<double>::epsilon()) {
111 return std::numeric_limits<double>::infinity();
113 return PI * 2.0 * Inertia() / angular;
117 double Body::DayLength() const noexcept {
119 return RotationalPeriod();
121 double year = OrbitalPeriod();
122 double sidereal = RotationalPeriod();
123 double grade = (angular < 0.0 ? -1.0 : 1.0) * (std::abs(axis_tilt.x) > PI * 0.5 ? -1.0 : 1.0);
124 return std::abs((year * sidereal) / ( year + (grade * sidereal)));
127 double Body::SphereOfInfluence() const noexcept {
129 return orbit.SemiMajorAxis() * std::pow(Mass() / Parent().Mass(), 2.0 / 5.0);
131 return std::numeric_limits<double>::infinity();
135 glm::dmat4 Body::ToUniverse() const noexcept {
137 const Body *b = this;
138 while (b->HasParent()) {
139 m = b->ToParent() * m;
145 glm::dmat4 Body::FromUniverse() const noexcept {
147 const Body *b = this;
148 while (b->HasParent()) {
149 m *= b->FromParent();
156 std::vector<creature::Creature *> ccache;
157 std::vector<CreatureCreatureCollision> collisions;
160 void Body::Tick(double dt) {
161 rotation += dt * AngularMomentum() / Inertia();
163 ccache = Creatures();
164 for (creature::Creature *c : ccache) {
167 // first remove creatures so they don't collide
168 for (auto c = Creatures().begin(); c != Creatures().end();) {
169 if ((*c)->Removable()) {
171 c = Creatures().erase(c);
179 void Body::Cache() noexcept {
182 orbit.Matrix(PI * 2.0 * (GetSimulation().Time() / OrbitalPeriod()))
183 * glm::eulerAngleXY(axis_tilt.x, axis_tilt.y);
185 glm::eulerAngleYX(-axis_tilt.y, -axis_tilt.x)
186 * orbit.InverseMatrix(PI * 2.0 * (GetSimulation().Time() / OrbitalPeriod()));
188 orbital = glm::eulerAngleXY(axis_tilt.x, axis_tilt.y);
189 inverse_orbital = glm::eulerAngleYX(-axis_tilt.y, -axis_tilt.x);
191 local = glm::eulerAngleY(rotation);
192 inverse_local = glm::eulerAngleY(-rotation);
195 void Body::CheckCollision() noexcept {
196 if (Creatures().size() < 2) return;
198 auto end = Creatures().end();
199 for (auto i = Creatures().begin(); i != end; ++i) {
200 math::AABB i_box((*i)->CollisionBounds());
201 glm::dmat4 i_mat((*i)->CollisionTransform());
202 for (auto j = (i + 1); j != end; ++j) {
203 glm::dvec3 diff((*i)->GetSituation().Position() - (*j)->GetSituation().Position());
204 double max_dist = ((*i)->Size() + (*j)->Size()) * 1.74;
205 if (glm::length2(diff) > max_dist * max_dist) continue;
206 math::AABB j_box((*j)->CollisionBounds());
207 glm::dmat4 j_mat((*j)->CollisionTransform());
210 if (Intersect(i_box, i_mat, j_box, j_mat, normal, depth)) {
211 collisions.push_back({ **i, **j, normal, depth });
215 for (auto &c : collisions) {
216 c.A().OnCollide(c.B());
217 c.B().OnCollide(c.A());
218 c.A().GetSituation().Move(c.Normal() * (c.Depth() * -0.5));
219 c.B().GetSituation().Move(c.Normal() * (c.Depth() * 0.5));
220 c.A().GetSituation().Accelerate(c.Normal() * -glm::dot(c.Normal(), c.AVel()));
221 c.B().GetSituation().Accelerate(c.Normal() * -glm::dot(c.Normal(), c.BVel()));
225 void Body::AddCreature(creature::Creature *c) {
226 creatures.push_back(c);
229 void Body::RemoveCreature(creature::Creature *c) {
230 auto entry = std::find(creatures.begin(), creatures.end(), c);
231 if (entry != creatures.end()) {
232 creatures.erase(entry);
237 CreatureCreatureCollision::~CreatureCreatureCollision() {
240 const glm::dvec3 &CreatureCreatureCollision::APos() const noexcept {
241 return a->GetSituation().Position();
244 const glm::dvec3 &CreatureCreatureCollision::AVel() const noexcept {
245 return a->GetSituation().Velocity();
248 const glm::dvec3 &CreatureCreatureCollision::BPos() const noexcept {
249 return b->GetSituation().Position();
252 const glm::dvec3 &CreatureCreatureCollision::BVel() const noexcept {
253 return b->GetSituation().Velocity();
269 double Orbit::SemiMajorAxis() const noexcept {
273 Orbit &Orbit::SemiMajorAxis(double s) noexcept {
278 double Orbit::Eccentricity() const noexcept {
282 Orbit &Orbit::Eccentricity(double e) noexcept {
287 double Orbit::Inclination() const noexcept {
291 Orbit &Orbit::Inclination(double i) noexcept {
296 double Orbit::LongitudeAscending() const noexcept {
300 Orbit &Orbit::LongitudeAscending(double l) noexcept {
305 double Orbit::ArgumentPeriapsis() const noexcept {
309 Orbit &Orbit::ArgumentPeriapsis(double a) noexcept {
314 double Orbit::MeanAnomaly() const noexcept {
318 Orbit &Orbit::MeanAnomaly(double m) noexcept {
325 double mean2eccentric(double M, double e) {
326 double E = M; // eccentric anomaly, solve M = E - e sin E
327 // limit to 100 steps to prevent deadlocks in impossible situations
328 for (int i = 0; i < 100; ++i) {
329 double dE = (E - e * sin(E) - M) / (1 - e * cos(E));
331 if (std::abs(dE) < 1.0e-6) break;
338 glm::dmat4 Orbit::Matrix(double t) const noexcept {
340 double E = mean2eccentric(M, ecc);
342 // coordinates in orbital plane, P=x, Q=-z
343 double P = sma * (cos(E) - ecc);
344 double Q = sma * sin(E) * sqrt(1 - (ecc * ecc));
346 return glm::yawPitchRoll(asc, inc, arg) * glm::translate(glm::dvec3(P, 0.0, -Q));
349 glm::dmat4 Orbit::InverseMatrix(double t) const noexcept {
351 double E = mean2eccentric(M, ecc);
352 double P = sma * (cos(E) - ecc);
353 double Q = sma * sin(E) * sqrt(1 - (ecc * ecc));
354 return glm::translate(glm::dvec3(-P, 0.0, Q)) * glm::transpose(glm::yawPitchRoll(asc, inc, arg));
358 Planet::Planet(int sidelength)
360 , sidelength(sidelength)
361 , tiles(TilesTotal())
363 Radius(double(sidelength) / 2.0);
370 /// map p onto cube, s gives the surface, u and v the position in [-1,1]
371 void cubemap(const glm::dvec3 &p, int &s, double &u, double &v) noexcept {
372 const glm::dvec3 p_abs(glm::abs(p));
373 const glm::bvec3 p_pos(glm::greaterThan(p, glm::dvec3(0.0)));
374 double max_axis = 0.0;
376 if (p_pos.x && p_abs.x >= p_abs.y && p_abs.x >= p_abs.z) {
382 if (!p_pos.x && p_abs.x >= p_abs.y && p_abs.x >= p_abs.z) {
388 if (p_pos.y && p_abs.y >= p_abs.x && p_abs.y >= p_abs.z) {
394 if (!p_pos.y && p_abs.y >= p_abs.x && p_abs.y >= p_abs.z) {
400 if (p_pos.z && p_abs.z >= p_abs.x && p_abs.z >= p_abs.y) {
406 if (!p_pos.z && p_abs.z >= p_abs.x && p_abs.z >= p_abs.y) {
415 /// get p from cube, s being surface, u and v the position in [-1,1],
416 /// gives a vector from the center to the surface
417 glm::dvec3 cubeunmap(int s, double u, double v) {
420 case 0: return glm::dvec3(u, v, 1.0); // +Z
421 case 1: return glm::dvec3(1.0, u, v); // +X
422 case 2: return glm::dvec3(v, 1.0, u); // +Y
423 case 3: return glm::dvec3(-u, -v, -1.0); // -Z
424 case 4: return glm::dvec3(-1.0, -u, -v); // -X
425 case 5: return glm::dvec3(-v, -1.0, -u); // -Y
430 Tile &Planet::TileAt(const glm::dvec3 &p) noexcept {
434 cubemap(p, srf, u, v);
435 int x = glm::clamp(int(u * Radius() + Radius()), 0, sidelength - 1);
436 int y = glm::clamp(int(v * Radius() + Radius()), 0, sidelength - 1);
437 return TileAt(srf, x, y);
440 const Tile &Planet::TileAt(const glm::dvec3 &p) const noexcept {
444 cubemap(p, srf, u, v);
445 int x = glm::clamp(int(u * Radius() + Radius()), 0, sidelength - 1);
446 int y = glm::clamp(int(v * Radius() + Radius()), 0, sidelength - 1);
447 return TileAt(srf, x, y);
450 const TileType &Planet::TileTypeAt(const glm::dvec3 &p) const noexcept {
451 return GetSimulation().TileTypes()[TileAt(p).type];
454 Tile &Planet::TileAt(int surface, int x, int y) noexcept {
455 return tiles[IndexOf(surface, x, y)];
458 const Tile &Planet::TileAt(int surface, int x, int y) const noexcept {
459 return tiles[IndexOf(surface, x, y)];
462 const TileType &Planet::TypeAt(int srf, int x, int y) const noexcept {
463 return GetSimulation().TileTypes()[TileAt(srf, x, y).type];
466 glm::dvec3 Planet::TileCenter(int srf, int x, int y, double e) const noexcept {
467 double u = (double(x) - Radius() + 0.5) / Radius();
468 double v = (double(y) - Radius() + 0.5) / Radius();
469 return glm::normalize(cubeunmap(srf, u, v)) * (Radius() + e);
472 void Planet::BuildVAO() {
473 vao.reset(new graphics::SimpleVAO<Attributes, unsigned int>);
475 vao->BindAttributes();
476 vao->EnableAttribute(0);
477 vao->EnableAttribute(1);
478 vao->EnableAttribute(2);
479 vao->EnableAttribute(3);
480 vao->EnableAttribute(4);
481 vao->EnableAttribute(5);
482 vao->AttributePointer<glm::vec3>(0, false, offsetof(Attributes, position));
483 vao->AttributePointer<glm::vec3>(1, false, offsetof(Attributes, normal));
484 vao->AttributePointer<glm::vec3>(2, false, offsetof(Attributes, tex_coord));
485 vao->AttributePointer<float>(3, false, offsetof(Attributes, shiny));
486 vao->AttributePointer<float>(4, false, offsetof(Attributes, glossy));
487 vao->AttributePointer<float>(5, false, offsetof(Attributes, metallic));
488 vao->ReserveAttributes(TilesTotal() * 4, GL_STATIC_DRAW);
490 auto attrib = vao->MapAttributes(GL_WRITE_ONLY);
491 float offset = Radius();
494 // up +Z +X +Y -Z -X -Y
496 for (int index = 0, surface = 0; surface < 6; ++surface) {
497 for (int y = 0; y < sidelength; ++y) {
498 for (int x = 0; x < sidelength; ++x, ++index) {
500 pos[0][(surface + 0) % 3] = float(x + 0) - offset;
501 pos[0][(surface + 1) % 3] = float(y + 0) - offset;
502 pos[0][(surface + 2) % 3] = offset;
503 pos[1][(surface + 0) % 3] = float(x + 0) - offset;
504 pos[1][(surface + 1) % 3] = float(y + 1) - offset;
505 pos[1][(surface + 2) % 3] = offset;
506 pos[2][(surface + 0) % 3] = float(x + 1) - offset;
507 pos[2][(surface + 1) % 3] = float(y + 0) - offset;
508 pos[2][(surface + 2) % 3] = offset;
509 pos[3][(surface + 0) % 3] = float(x + 1) - offset;
510 pos[3][(surface + 1) % 3] = float(y + 1) - offset;
511 pos[3][(surface + 2) % 3] = offset;
513 const TileType &t = TypeAt(surface, x, y);
514 float tex = t.texture;
515 const float tex_v_begin = surface < 3 ? 1.0f : 0.0f;
516 const float tex_v_end = surface < 3 ? 0.0f : 1.0f;
518 attrib[4 * index + 0].position = glm::normalize(pos[0]) * (surface < 3 ? offset : -offset);
519 attrib[4 * index + 0].normal = pos[0];
520 attrib[4 * index + 0].tex_coord[0] = 0.0f;
521 attrib[4 * index + 0].tex_coord[1] = tex_v_begin;
522 attrib[4 * index + 0].tex_coord[2] = tex;
523 attrib[4 * index + 0].shiny = t.shiny;
524 attrib[4 * index + 0].glossy = t.glossy;
525 attrib[4 * index + 0].metallic = t.metallic;
527 attrib[4 * index + 1].position = glm::normalize(pos[1]) * (surface < 3 ? offset : -offset);
528 attrib[4 * index + 1].normal = pos[1];
529 attrib[4 * index + 1].tex_coord[0] = 0.0f;
530 attrib[4 * index + 1].tex_coord[1] = tex_v_end;
531 attrib[4 * index + 1].tex_coord[2] = tex;
532 attrib[4 * index + 1].shiny = t.shiny;
533 attrib[4 * index + 1].glossy = t.glossy;
534 attrib[4 * index + 1].metallic = t.metallic;
536 attrib[4 * index + 2].position = glm::normalize(pos[2]) * (surface < 3 ? offset : -offset);
537 attrib[4 * index + 2].normal = pos[2];
538 attrib[4 * index + 2].tex_coord[0] = 1.0f;
539 attrib[4 * index + 2].tex_coord[1] = tex_v_begin;
540 attrib[4 * index + 2].tex_coord[2] = tex;
541 attrib[4 * index + 2].shiny = t.shiny;
542 attrib[4 * index + 2].glossy = t.glossy;
543 attrib[4 * index + 2].metallic = t.metallic;
545 attrib[4 * index + 3].position = glm::normalize(pos[3]) * (surface < 3 ? offset : -offset);
546 attrib[4 * index + 3].normal = pos[3];
547 attrib[4 * index + 3].tex_coord[0] = 1.0f;
548 attrib[4 * index + 3].tex_coord[1] = tex_v_end;
549 attrib[4 * index + 3].tex_coord[2] = tex;
550 attrib[4 * index + 3].shiny = t.shiny;
551 attrib[4 * index + 3].glossy = t.glossy;
552 attrib[4 * index + 3].metallic = t.metallic;
558 vao->ReserveElements(TilesTotal() * 6, GL_STATIC_DRAW);
560 auto element = vao->MapElements(GL_WRITE_ONLY);
562 for (int surface = 0; surface < 3; ++surface) {
563 for (int y = 0; y < sidelength; ++y) {
564 for (int x = 0; x < sidelength; ++x, ++index) {
565 element[6 * index + 0] = 4 * index + 0;
566 element[6 * index + 1] = 4 * index + 2;
567 element[6 * index + 2] = 4 * index + 1;
568 element[6 * index + 3] = 4 * index + 1;
569 element[6 * index + 4] = 4 * index + 2;
570 element[6 * index + 5] = 4 * index + 3;
574 for (int surface = 3; surface < 6; ++surface) {
575 for (int y = 0; y < sidelength; ++y) {
576 for (int x = 0; x < sidelength; ++x, ++index) {
577 element[6 * index + 0] = 4 * index + 0;
578 element[6 * index + 1] = 4 * index + 1;
579 element[6 * index + 2] = 4 * index + 2;
580 element[6 * index + 3] = 4 * index + 2;
581 element[6 * index + 4] = 4 * index + 1;
582 element[6 * index + 5] = 4 * index + 3;
590 void Planet::Draw(app::Assets &assets, graphics::Viewport &viewport) {
594 vao->DrawTriangles(TilesTotal() * 6);
598 void GenerateEarthlike(const Set<TileType> &tiles, Planet &p) noexcept {
599 math::SimplexNoise elevation_gen(0);
600 math::SimplexNoise variation_gen(45623752346);
602 const int ice = tiles["ice"].id;
603 const int ocean = tiles["ocean"].id;
604 const int water = tiles["water"].id;
605 const int sand = tiles["sand"].id;
606 const int grass = tiles["grass"].id;
607 const int tundra = tiles["tundra"].id;
608 const int taiga = tiles["taiga"].id;
609 const int desert = tiles["desert"].id;
610 const int mntn = tiles["mountain"].id;
611 const int algae = tiles["algae"].id;
612 const int forest = tiles["forest"].id;
613 const int jungle = tiles["jungle"].id;
614 const int rock = tiles["rock"].id;
615 const int wheat = tiles["wheat"].id;
617 constexpr double ocean_thresh = -0.2;
618 constexpr double water_thresh = 0.0;
619 constexpr double beach_thresh = 0.05;
620 constexpr double highland_thresh = 0.4;
621 constexpr double mountain_thresh = 0.5;
623 const glm::dvec3 axis(0.0, 1.0, 0.0);
624 const double cap_thresh = std::abs(std::cos(p.AxialTilt().x));
625 const double equ_thresh = std::abs(std::sin(p.AxialTilt().x)) / 2.0;
626 const double fzone_start = equ_thresh - (equ_thresh - cap_thresh) / 3.0;
627 const double fzone_end = cap_thresh + (equ_thresh - cap_thresh) / 3.0;
629 for (int surface = 0; surface <= 5; ++surface) {
630 for (int y = 0; y < p.SideLength(); ++y) {
631 for (int x = 0; x < p.SideLength(); ++x) {
632 glm::dvec3 to_tile = p.TileCenter(surface, x, y);
633 double near_axis = std::abs(glm::dot(glm::normalize(to_tile), axis));
634 if (near_axis > cap_thresh) {
635 p.TileAt(surface, x, y).type = ice;
638 float elevation = math::OctaveNoise(
640 glm::vec3(to_tile / p.Radius()),
643 5 / p.Radius(), // frequency
647 float variation = math::OctaveNoise(
649 glm::vec3(to_tile / p.Radius()),
652 16 / p.Radius(), // frequency
656 if (elevation < ocean_thresh) {
657 p.TileAt(surface, x, y).type = ocean;
658 } else if (elevation < water_thresh) {
659 if (variation > 0.3) {
660 p.TileAt(surface, x, y).type = algae;
662 p.TileAt(surface, x, y).type = water;
664 } else if (elevation < beach_thresh) {
665 p.TileAt(surface, x, y).type = sand;
666 } else if (elevation < highland_thresh) {
667 if (near_axis < equ_thresh) {
668 if (variation > 0.6) {
669 p.TileAt(surface, x, y).type = grass;
670 } else if (variation > 0.2) {
671 p.TileAt(surface, x, y).type = sand;
673 p.TileAt(surface, x, y).type = desert;
675 } else if (near_axis < fzone_start) {
676 if (variation > 0.4) {
677 p.TileAt(surface, x, y).type = forest;
678 } else if (variation < -0.5) {
679 p.TileAt(surface, x, y).type = jungle;
680 } else if (variation > -0.02 && variation < 0.02) {
681 p.TileAt(surface, x, y).type = wheat;
683 p.TileAt(surface, x, y).type = grass;
685 } else if (near_axis < fzone_end) {
686 p.TileAt(surface, x, y).type = tundra;
688 p.TileAt(surface, x, y).type = taiga;
690 } else if (elevation < mountain_thresh) {
691 if (variation > 0.3) {
692 p.TileAt(surface, x, y).type = mntn;
694 p.TileAt(surface, x, y).type = rock;
697 p.TileAt(surface, x, y).type = mntn;
705 void GenerateTest(const Set<TileType> &tiles, Planet &p) noexcept {
706 for (int surface = 0; surface <= 5; ++surface) {
707 for (int y = 0; y < p.SideLength(); ++y) {
708 for (int x = 0; x < p.SideLength(); ++x) {
709 if (x == p.SideLength() / 2 && y == p.SideLength() / 2) {
710 p.TileAt(surface, x, y).type = surface;
712 p.TileAt(surface, x, y).type = (x == p.SideLength()/2) + (y == p.SideLength()/2) + 6;
731 std::vector<TileType::Yield>::const_iterator TileType::FindResource(int r) const {
732 auto yield = resources.cbegin();
733 for (; yield != resources.cend(); ++yield) {
734 if (yield->resource == r) {
741 std::vector<TileType::Yield>::const_iterator TileType::FindBestResource(const creature::Composition &comp) const {
742 auto best = resources.cend();
743 double best_value = 0.0;
744 for (auto yield = resources.cbegin(); yield != resources.cend(); ++yield) {
745 double value = comp.Get(yield->resource);
746 if (value > best_value) {