eclipsefinder.cpp
上传用户:center1979
上传日期:2022-07-26
资源大小:50633k
文件大小:8k
源码类别:

OpenGL

开发平台:

Visual C++

  1. // eclipsefinder.cpp by Christophe Teyssier <chris@teyssier.org>
  2. // adapted form wineclipses.cpp by Kendrix <kendrix@wanadoo.fr>
  3. // 
  4. // Copyright (C) 2001, Chris Laurel <claurel@shatters.net>
  5. //
  6. // Compute Solar Eclipses for our Solar System planets
  7. //
  8. // This program is free software; you can redistribute it and/or
  9. // modify it under the terms of the GNU General Public License
  10. // as published by the Free Software Foundation; either version 2
  11. // of the License, or (at your option) any later version.
  12. #include <cstring>
  13. #include <sstream>
  14. #include <algorithm>
  15. #include <set>
  16. #include <cassert>
  17. #include "eclipsefinder.h"
  18. #include "celmath/mathlib.h"
  19. #include "celmath/ray.h"
  20. #include "celmath/distance.h"
  21. using namespace std;
  22. Eclipse::Eclipse(int Y, int M, int D) :
  23.     body(NULL)
  24. {
  25.     date = new astro::Date(Y, M, D);
  26. }
  27. Eclipse::Eclipse(double JD) :
  28.     body(NULL)
  29. {
  30.     date = new astro::Date(JD);
  31. }
  32. // TODO: share this constant and function with render.cpp
  33. static const float MinRelativeOccluderRadius = 0.005f;
  34. bool EclipseFinder::testEclipse(const Body& receiver, const Body& caster,
  35.                                 double now) const
  36. {
  37.     // Ignore situations where the shadow casting body is much smaller than
  38.     // the receiver, as these shadows aren't likely to be relevant.  Also,
  39.     // ignore eclipses where the caster is not an ellipsoid, since we can't
  40.     // generate correct shadows in this case.
  41.     if (caster.getRadius() >= receiver.getRadius() * MinRelativeOccluderRadius &&
  42.         caster.isEllipsoid())
  43.     {
  44.         // All of the eclipse related code assumes that both the caster
  45.         // and receiver are spherical.  Irregular receivers will work more
  46.         // or less correctly, but casters that are sufficiently non-spherical
  47.         // will produce obviously incorrect shadows.  Another assumption we
  48.         // make is that the distance between the caster and receiver is much
  49.         // less than the distance between the sun and the receiver.  This
  50.         // approximation works everywhere in the solar system, and likely
  51.         // works for any orbitally stable pair of objects orbiting a star.
  52.         Point3d posReceiver = receiver.getAstrocentricPosition(now);
  53.         Point3d posCaster = caster.getAstrocentricPosition(now);
  54.         const Star* sun = receiver.getSystem()->getStar();
  55.         assert(sun != NULL);
  56.         double distToSun = posReceiver.distanceFromOrigin();
  57.         float appSunRadius = (float) (sun->getRadius() / distToSun);
  58.         Vec3d dir = posCaster - posReceiver;
  59.         double distToCaster = dir.length() - receiver.getRadius();
  60.         float appOccluderRadius = (float) (caster.getRadius() / distToCaster);
  61.         // The shadow radius is the radius of the occluder plus some additional
  62.         // amount that depends upon the apparent radius of the sun.  For
  63.         // a sun that's distant/small and effectively a point, the shadow
  64.         // radius will be the same as the radius of the occluder.
  65.         float shadowRadius = (1 + appSunRadius / appOccluderRadius) *
  66.             caster.getRadius();
  67.         // Test whether a shadow is cast on the receiver.  We want to know
  68.         // if the receiver lies within the shadow volume of the caster.  Since
  69.         // we're assuming that everything is a sphere and the sun is far
  70.         // away relative to the caster, the shadow volume is a
  71.         // cylinder capped at one end.  Testing for the intersection of a
  72.         // singly capped cylinder is as simple as checking the distance
  73.         // from the center of the receiver to the axis of the shadow cylinder.
  74.         // If the distance is less than the sum of the caster's and receiver's
  75.         // radii, then we have an eclipse.
  76.         float R = receiver.getRadius() + shadowRadius;
  77.         double dist = distance(posReceiver,
  78.                                Ray3d(posCaster, posCaster - Point3d(0, 0, 0)));
  79.         if (dist < R)
  80.         {
  81.             // Ignore "eclipses" where the caster and receiver have
  82.             // intersecting bounding spheres.
  83.             if (distToCaster > caster.getRadius())
  84.                 return true;
  85.         }
  86.     }
  87.     return false;
  88. }
  89. double EclipseFinder::findEclipseSpan(const Body& receiver, const Body& caster,
  90.                        double now, double dt) const
  91. {
  92.     double t = now;
  93.     while (testEclipse(receiver, caster, t))
  94.         t += dt;
  95.     return t;
  96. }
  97. int EclipseFinder::CalculateEclipses()
  98. {
  99.     Simulation* sim = appCore->getSimulation();
  100.     
  101.     Eclipse* eclipse;
  102.     double* JDback = NULL;
  103.     
  104.     int nIDplanetetofindon = 0;
  105.     int nSattelites = 0;
  106.     const SolarSystem* sys = sim->getNearestSolarSystem();
  107.     
  108.     toProcess = false;
  109.     
  110.     if ((!sys))
  111.     {
  112.         eclipse = new Eclipse(0.);
  113.         eclipse->planete = "None";
  114.         Eclipses_.insert(Eclipses_.end(), *eclipse);
  115.         delete eclipse;
  116.         return 1;
  117.     }
  118.     else if (sys->getStar()->getCatalogNumber() != 0)
  119.     {
  120.         eclipse = new Eclipse(0.);
  121.         eclipse->planete = "None";
  122.         Eclipses_.insert(Eclipses_.end(), *eclipse);
  123.         delete eclipse;
  124.         return 1;
  125.     }
  126.     PlanetarySystem* system = sys->getPlanets();
  127.     int nbPlanets = system->getSystemSize();
  128.     for (int i = 0; i < nbPlanets; ++i)
  129.     {
  130.         Body* planete = system->getBody(i);
  131.         if (planete != NULL)
  132.             if (strPlaneteToFindOn == planete->getName())
  133.             {
  134.                 nIDplanetetofindon = i;
  135.                 PlanetarySystem* satellites = planete->getSatellites();
  136.                 if (satellites)
  137.                 {
  138.                     nSattelites = satellites->getSystemSize();
  139.                     JDback = new double[nSattelites];
  140.                     memset(JDback, 0, nSattelites*sizeof(double));
  141.                 }
  142.                 break;
  143.             }
  144.     }
  145.     Body* planete = system->getBody(nIDplanetetofindon);
  146.     while (JDfrom < JDto)
  147.     {
  148.         PlanetarySystem* satellites = planete->getSatellites();
  149.         if (satellites)
  150.         {
  151.             for (int j = 0; j < nSattelites; ++j)
  152.             {
  153.                 Body* caster = NULL;
  154.                 Body* receiver = NULL;
  155.                 bool test = false;
  156.                 if (satellites->getBody(j)->getClassification() != Body::Spacecraft)
  157.                 {
  158.                     if (type == Eclipse::Solar)
  159.                     {
  160.                         caster = satellites->getBody(j);
  161.                         receiver = planete;
  162.                     }
  163.                     else
  164.                     {
  165.                         caster = planete;
  166.                         receiver = satellites->getBody(j);
  167.                     }
  168.                     test = testEclipse(*receiver, *caster, JDfrom);
  169.                 }
  170.                 if (test && JDfrom - JDback[j] > 1)
  171.                 {
  172.                     JDback[j] = JDfrom;
  173.                     eclipse = new Eclipse(JDfrom);
  174.                     eclipse->startTime = findEclipseSpan(*receiver, *caster,
  175.                                                          JDfrom,
  176.                                                          -1.0 / (24.0 * 60.0));
  177.                     eclipse->endTime   = findEclipseSpan(*receiver, *caster,
  178.                                                          JDfrom,
  179.                                                          1.0 / (24.0 * 60.0));
  180.                     eclipse->body = receiver;
  181.                     eclipse->planete = planete->getName();
  182.                     eclipse->sattelite = satellites->getBody(j)->getName();
  183.                     Eclipses_.insert(Eclipses_.end(), *eclipse);
  184.                     delete eclipse;
  185.                 }
  186.             }
  187.         }
  188.         JDfrom += 1.0 / 24.0;
  189.     }
  190.     if (JDback)
  191.         delete JDback;
  192.     if (Eclipses_.empty())
  193.     {
  194.         eclipse = new Eclipse(0.);
  195.         eclipse->planete = "None";
  196.         Eclipses_.insert(Eclipses_.end(), *eclipse);
  197.         delete eclipse;
  198.     }
  199.     return 0;
  200. }