Fixed raytracer problems due to numeric scale.

I had lots of problems with using AU as my scale units.
By changing to 10000*km units, the vector equation solver
works correctly again, and I actually get an image of Jupiter
and its moons. However, it does not match the test photo:
Ganymede does not appear close to the planet, nor is there
a shadow of it on the planet. I will have to debug that
separately.

demo/c/raytrace/debug.cpp

@@ -0,0 +1,88 @@
+/*
+    debug.cpp
+
+    Copyright (C) 2013 by Don Cross  -  http://cosinekitty.com/raytrace
+
+    This software is provided 'as-is', without any express or implied
+    warranty. In no event will the author be held liable for any damages
+    arising from the use of this software.
+
+    Permission is granted to anyone to use this software for any purpose,
+    including commercial applications, and to alter it and redistribute it
+    freely, subject to the following restrictions:
+
+    1. The origin of this software must not be misrepresented; you must not
+       claim that you wrote the original software. If you use this software
+       in a product, an acknowledgment in the product documentation would be
+       appreciated but is not required.
+
+    2. Altered source versions must be plainly marked as such, and must not be
+       misrepresented as being the original software.
+
+    3. This notice may not be removed or altered from any source
+       distribution.
+
+    -------------------------------------------------------------------------
+
+    Implements output operators and other stuff 
+    useful for debugging the ray tracer.
+*/
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include "imager.h"
+
+namespace Imager
+{
+    void Indent(std::ostream& output, int depth)
+    {
+        const int numSpaces = 4 * depth;
+        for (int i=0; i < numSpaces; ++i)
+        {
+            output << ' ';
+        }
+    }
+
+    std::ostream& operator<< (std::ostream &output, const Color& color)
+    {
+        using namespace std;
+
+        output.precision(3);
+        output << scientific;
+        output << "(R " << color.red << ", G " << color.green << ", B " << color.blue << ")";
+        return output;
+    }
+
+    std::ostream& operator<< (std::ostream& output, const Intersection& intersection)
+    {
+        using namespace std;
+
+        output << "{ ";
+        if (intersection.solid)
+        {
+            output << intersection.solid->GetTag();
+        }
+        if (intersection.tag)
+        {
+            output << ":" << intersection.tag;
+        }
+        output.precision(3);
+        output << fixed;
+        output << " d=" << intersection.distanceSquared;
+        output << ", p=" << intersection.point;
+        output << ", n=" << intersection.surfaceNormal;
+        output << " }";
+        return output;
+    }
+
+    std::ostream& operator<< (std::ostream& output, const Vector& vector)
+    {
+        using namespace std;
+
+        output.precision(3);
+        output << fixed;
+        output << "(" << vector.x << ", " << vector.y << ", " << vector.z << ")";
+        return output;
+    }
+}

demo/c/raytrace/main.cpp

@@ -88,12 +88,16 @@ int JupiterImage(const char *filename, int width, astro_time_t time)
 
     Scene scene(Color(0.0, 0.0, 0.0));
 
-    const double equ_radius_au = JUPITER_EQUATORIAL_RADIUS_KM / KM_PER_AU;
-    const double pol_radius_au = JUPITER_POLAR_RADIUS_KM / KM_PER_AU;
-    Spheroid *planet = new Spheroid(equ_radius_au, equ_radius_au, pol_radius_au);
+    const double km_scale = 10000.0;   // makes Jupiter's radius about 7 units.
+    const double au_scale = km_scale / KM_PER_AU;
+
+    const double equ_radius = JUPITER_EQUATORIAL_RADIUS_KM / km_scale;
+    const double pol_radius = JUPITER_POLAR_RADIUS_KM / km_scale;
+    Spheroid *planet = new Spheroid(equ_radius, equ_radius, pol_radius);
     scene.AddSolidObject(planet);
     planet->SetFullMatte(Color(1.0, 1.0, 0.7));
-    planet->Move(Vector(jupiter.x, jupiter.y, jupiter.z));
+    planet->Move(Vector(jupiter.x, jupiter.y, jupiter.z) / au_scale);
+    planet->SetTag("Jupiter");
 
     // Reorient Jupiter's rotation axis to match the calculated orientation.
     planet->RotateY(90.0 - axis.dec);       // ?? not sure about direction
@@ -102,51 +106,54 @@ int JupiterImage(const char *filename, int width, astro_time_t time)
     // Add Jupiter's moons to the scene.
 
     // Io
-    const double io_radius = IO_RADIUS_KM / KM_PER_AU;
+    const double io_radius = IO_RADIUS_KM / km_scale;
     Spheroid *io = new Spheroid(io_radius, io_radius, io_radius);
     scene.AddSolidObject(io);
     io->SetFullMatte(Color(1.0, 1.0, 1.0));     // ??? Actual moon colors ???
     io->Move(Vector(
         jm.moon[JM_IO].x + jupiter.x,
         jm.moon[JM_IO].y + jupiter.y,
-        jm.moon[JM_IO].z + jupiter.z)
+        jm.moon[JM_IO].z + jupiter.z) / au_scale
     );
+    io->SetTag("Io");
 
     // Europa
-    const double eu_radius = EUROPA_RADIUS_KM / KM_PER_AU;
+    const double eu_radius = EUROPA_RADIUS_KM / km_scale;
     Spheroid *europa = new Spheroid(eu_radius, eu_radius, eu_radius);
     scene.AddSolidObject(europa);
     europa->SetFullMatte(Color(1.0, 1.0, 1.0));     // ??? Actual moon colors ???
     europa->Move(Vector(
         jm.moon[JM_EUROPA].x + jupiter.x,
         jm.moon[JM_EUROPA].y + jupiter.y,
-        jm.moon[JM_EUROPA].z + jupiter.z)
+        jm.moon[JM_EUROPA].z + jupiter.z) / au_scale
     );
+    europa->SetTag("Europa");
 
     // Ganymede
-    const double gan_radius = GANYMEDE_RADIUS_KM / KM_PER_AU;
+    const double gan_radius = GANYMEDE_RADIUS_KM / km_scale;
     Spheroid *ganymede = new Spheroid(gan_radius, gan_radius, gan_radius);
     scene.AddSolidObject(ganymede);
     ganymede->SetFullMatte(Color(1.0, 1.0, 1.0));     // ??? Actual moon colors ???
     ganymede->Move(Vector(
         jm.moon[JM_GANYMEDE].x + jupiter.x,
         jm.moon[JM_GANYMEDE].y + jupiter.y,
-        jm.moon[JM_GANYMEDE].z + jupiter.z)
+        jm.moon[JM_GANYMEDE].z + jupiter.z) / au_scale
     );
 
     // Callisto
-    const double cal_radius = CALLISTO_RADIUS_KM / KM_PER_AU;
+    const double cal_radius = CALLISTO_RADIUS_KM / km_scale;
     Spheroid *callisto = new Spheroid(cal_radius, cal_radius, cal_radius);
     scene.AddSolidObject(callisto);
     callisto->SetFullMatte(Color(1.0, 1.0, 1.0));     // ??? Actual moon colors ???
     callisto->Move(Vector(
         jm.moon[JM_CALLISTO].x + jupiter.x,
         jm.moon[JM_CALLISTO].y + jupiter.y,
-        jm.moon[JM_CALLISTO].z + jupiter.z)
+        jm.moon[JM_CALLISTO].z + jupiter.z) / au_scale
     );
+    callisto->SetTag("Callisto");
 
     // Add the Sun as the point light source.
-    scene.AddLightSource(LightSource(Vector(sun.x, sun.y, sun.z), Color(1.0, 1.0, 1.0)));
+    scene.AddLightSource(LightSource(Vector(sun.x, sun.y, sun.z) / au_scale, Color(1.0, 1.0, 1.0)));
 
     // Aim the camera at the center of Jupiter.
     // Start with an identity matrix, which leaves the camera pointing in the -z direction,
@@ -178,7 +185,7 @@ int JupiterImage(const char *filename, int width, astro_time_t time)
 
     // Magnify the image as appropriate for Jupiter's closest approach (opposition).
     // Render the image.
-    const double zoom = 200.0;
+    const double zoom = 300.0;
     scene.SaveImage(filename, (size_t)width, (size_t)width, zoom, 4);
 
     return 0;

demo/c/raytrace/run

@@ -0,0 +1,2 @@
+#!/bin/bash
+./raytrace jupiter.png 1000 Jupiter 2021-12-11T23:12:00Z

demo/c/raytrace/spheroid.cpp

@@ -24,7 +24,7 @@
 
     -------------------------------------------------------------------------
     Implements class Spheroid.
-    A spheroid is like a sphere, only it may have three different 
+    A spheroid is like a sphere, only it may have three different
     diameters in the x-, y-, and z-directions.
 */
 
@@ -34,8 +34,8 @@
 namespace Imager
 {
     void Spheroid::ObjectSpace_AppendAllIntersections(
-        const Vector& vantage, 
-        const Vector& direction, 
+        const Vector& vantage,
+        const Vector& direction,
         IntersectionList& intersectionList) const
     {
         double u[2];
@@ -55,9 +55,9 @@ namespace Imager
                 intersection.distanceSquared = displacement.MagnitudeSquared();
                 intersection.point = vantage + displacement;
 
-                // The surface normal vector was calculated by expressing the spheroid as a 
+                // The surface normal vector was calculated by expressing the spheroid as a
                 // function z(x,y) = sqrt(1 - (x/a)^2 - (y/b)^2),
-                // taking partial derivatives dz/dx = (c*c*x)/(a*a*z), dz/dy = (c*c*y)/(b*b*z), 
+                // taking partial derivatives dz/dx = (c*c*x)/(a*a*z), dz/dy = (c*c*y)/(b*b*z),
                 // and using these to calculate the vectors <1, 0, dz/dx> and <0, 1, dy,dz>.
                 // The normalized cross product of these two vectors yields the surface normal vector.
                 const double x = intersection.point.x;
@@ -74,7 +74,7 @@ namespace Imager
                     }
                     else
                     {
-                        // The equation devolves to an ellipse on the xy plane : 
+                        // The equation devolves to an ellipse on the xy plane :
                         // (x^2)/(a^2) + (y^2)/(b^2) = 1.
                         intersection.surfaceNormal = Vector(-1.0, -(a2*y)/(b2*x), 0.0).UnitVector();
                     }