feat: preparation for self shadowing

didymos_hera.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "id": "5fcbc291-92ea-4233-a7ad-f335293b327a",
    "metadata": {},
    "outputs": [],
@@ -23,7 +23,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "id": "3b5cc6c6-cd81-47b6-8593-6a5aaf2ca473",
    "metadata": {},
    "outputs": [],
@@ -35,7 +35,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "id": "aa8ad5dd-c343-49eb-98b0-f2b007c8223b",
    "metadata": {},
    "outputs": [],

draft.ipynb

@@ -77,6 +77,40 @@
     "    def build() -> Asteroid:\n",
     "        return Asteroid(self.__spice_info, self.__mesh, self.__plot)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "id": "0b0b6eb3-0a99-4ca8-ab51-c676ed9e66b9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ET_ = world.ephemeris_time\n",
+    "FRAME_PRIMARY = didymos.frame\n",
+    "FRAME_SECONDARY = dimorphos.frame\n",
+    "NAME_PRIMARY = didymos.name\n",
+    "NAME_SECONDARY = dimorphos.name\n",
+    "ABCORR_NONE = \"NONE\"\n",
+    "ABCORR = \"LT+S\"\n",
+    "CENTERS_PRIMARY = didymos.center_facets()\n",
+    "CENTERS_SECONDARY = dimorphos.center_facets()\n",
+    "DSK_HANDLE_PRIMARY = didymos.handle_dsk()\n",
+    "DSK_HANDLE_SECONDARY = dimorphos.handle_dsk()\n",
+    "\n",
+    "\n",
+    "ORIENTATION_PRIMARY_FROM_FRAME_SECONDARY = spice.pxform(FRAME_PRIMARY, FRAME_SECONDARY, ET_)\n",
+    "POSITION_PRIMARY_FROM_FRAME_SECONDARY, _ = spice.spkpos(NAME_PRIMARY, ET_, FRAME_SECONDARY, ABCORR_NONE, NAME_SECONDARY)\n",
+    "POSITION_SUN_FROM_FRAME_SECONDARY, _ = spice.spkpos(\"SUN\", ET_, FRAME_SECONDARY, ABCORR, NAME_SECONDARY)\n",
+    "LIGHTS_PRIMARY_CURRENT = mesh.compute_light(CENTERS_PRIMARY, DSK_HANDLE_SECONDARY, POSITION_PRIMARY_FROM_FRAME_SECONDARY, ORIENTATION_PRIMARY_FROM_FRAME_SECONDARY, POSITION_SUN_FROM_FRAME_SECONDARY)\n",
+    "\n",
+    "ORIENTATION_SECONDARY_FROM_FRAME_PRIMARY = spice.pxform(FRAME_SECONDARY, FRAME_PRIMARY, ET_)\n",
+    "POSITION_SECONDARY_FROM_FRAME_PRIMARY, _ = spice.spkpos(NAME_SECONDARY, ET_, FRAME_PRIMARY, ABCORR_NONE, NAME_PRIMARY)\n",
+    "POSITION_SUN_FROM_FRAME_PRIMARY, _ = spice.spkpos(\"SUN\", ET_, FRAME_PRIMARY, ABCORR, NAME_PRIMARY)\n",
+    "LIGHTS_SECONDARY_CURRENT = mesh.compute_light(CENTERS_SECONDARY, DSK_HANDLE_PRIMARY, POSITION_SECONDARY_FROM_FRAME_PRIMARY, ORIENTATION_SECONDARY_FROM_FRAME_PRIMARY, POSITION_SUN_FROM_FRAME_PRIMARY)\n",
+    "\n",
+    "didymos.mesh[SCALAR] = LIGHTS_PRIMARY_CURRENT\n",
+    "dimorphos.mesh[SCALAR] = LIGHTS_SECONDARY_CURRENT"
+   ]
   }
  ],
  "metadata": {

sptool.py

+from __future__ import annotations
+
 from pathlib import Path
 
 import numpy
@@ -10,6 +12,10 @@ import tool
 import mesh
 
 
+def state(target: str, ephemeris_time: float, frame: str, abcorr: str, origin: str) -> tool.ARRAY_STATE:
+    return numpy.stack(spice.spkezr(target, ephemeris_time, frame, abcorr, origin)[0], axis=0)
+
+
 class Settings:
     def __init__(self, off_screen: bool = False) -> None:
         self.off_screen = off_screen
@@ -37,28 +43,63 @@ class Camera:
     def boundaries(self) -> tool.ARRAY_MATRIX4:
         _, _, _, _, BOUNDARIES = spice.getfvn(self.name, self.NUMBER_BOUNDS)
         return BOUNDARIES
-    
-    
+        
 
-class Spacecraft:
-    def __init__(self, name: str, frame: str, camera: Camera | None = None) -> None:
+class Body:
+    def __init__(self, world: "World", name: str, frame: str) -> None:
+        self.world = world  
         self.name = name
         self.frame = frame
+    
+    def __eq__(self, other: Body) -> bool:
+        return self.name.lower() == other.name.lower()
+
+    def state(self, frame: str | None = None, origin: str | None = None, ephemeris_time: float | None = None, abcorr: str = "NONE") -> tool.ARRAY_STATE:
+        if frame is None:
+            frame = self.world.frame
+        if origin is None:
+            origin = self.world.origin
+        if ephemeris_time is None:
+            ephemeris_time = self.world.ephemeris_time
+        return state(self.name, ephemeris_time, frame, abcorr, origin)
+    
+    def state_from(self, other: Body, frame: str | None = None, ephemeris_time: float | None = None, abcorr: str = "NONE") -> tool.ARRAY_STATE:
+        if frame is None:
+            frame = other.frame
+        if ephemeris_time is None:
+            ephemeris_time = self.world.ephemeris_time
+        return state(self.name, ephemeris_time, frame, abcorr, other.name)
+
+    def state_sun(self, frame: str | None = None, ephemeris_time: float | None = None, abcorr: str = "NONE") -> tool.ARRAY_STATE:
+        if frame is None:
+            frame = self.frame
+        if ephemeris_time is None:
+            ephemeris_time = self.world.ephemeris_time
+        return state("SUN", ephemeris_time, frame, abcorr, self.name)
+    
+    def orientation(self) -> tool.ARRAY_MATRIX:
+        return spice.pxform(self.frame, self.world.frame, self.world.ephemeris_time)
+
+    def orientation_from(self, other: Body) -> tool.ARRAY_MATRIX:
+        return spice.pxform(self.frame, other.frame, self.world.ephemeris_time)
+
+    
+class Spacecraft(Body):
+    def __init__(self, world: "World", name: str, frame: str, camera: Camera | None = None) -> None:
+        super().__init__(world, name, frame)
         self.camera: Camera | None = camera
         
     def set_camera(self, name: str, fov: (float, float), resolution: (float, float), vector_up: tool.ARRAY_VECTOR = numpy.array([0.0, 0.0, 1.0])) -> None:
         self.camera = Camera(name, fov, resolution, vector_up)
 
         
-class Asteroid:
-    def __init__(self, name: str, rotation_period: float, frame: str | None = None) -> None:
-        self.name = name
+class Asteroid(Body):
+    def __init__(self, world: "World", name: str, rotation_period: float, frame: str | None = None) -> None:
         self.rotation_period = rotation_period
         
         if frame is None:
-            self.frame = spice.cnmfrm(name)[1]
-        else:
-            self.frame = frame
+            frame = spice.cnmfrm(name)[1]
+        super().__init__(world, name, frame)
         
         # self.mesh: pyvista.DataSet | None = None
         # will be defined after the following method.
@@ -109,35 +150,67 @@ class Asteroid:
     def center_facets(self) -> pyvista.pyvista_ndarray[("N", 3), numpy.float64]:
         return self.mesh.cell_centers().points
 
-    def compute_visibility(self, OTHER: "Asteroid", POSITION_OBSERVER: tool.ARRAY_VECTOR, EPHEMERIS_TIME: float) -> tool.ARRAY_LIST:
-        """Spice distances are in kilometers. 1.0 means visible, 0.0 means not."""
+    def compute_visibility(self, other: Asteroid, position_observer: tool.ARRAY_VECTOR) -> tool.ARRAY_LIST:
+        """Spice distances are in kilometers. 1.0 means visible, 0.0 means not.
+        `position_observer` should be defined with respect to the other asteroid in it's fixed frame."""
         # TODO: parallelize facet loop.
         
-        FRAME_REF = "J2000ECLIP"
-        ABCORR_NONE = "NONE"
         ABCORR = "LT+S"
+        lights = numpy.ones(self.center_facets().shape[0])
         
-        shadows = numpy.ones(self.center_facets().shape[0])
-        
-        POSITION_FROM_OTHER = spice.spkpos(self.name, EPHEMERIS_TIME, OTHER.frame, ABCORR_NONE, OTHER.name)[0]
-        if tool.angle_between(POSITION_FROM_OTHER, POSITION_OBSERVER) < numpy.pi / 2.0:
-            return shadows
+        STATE_FROM_OTHER_FRAME = self.state_from(other)
+        if tool.angle_between(STATE_FROM_OTHER_FRAME[:3], position_observer) < numpy.pi / 2.0:
+            return lights
 
-        HANDLE_DSK_OTHER = OTHER.handle_dsk()
-        DLADSC_DSK_OTHER = OTHER.dladsc_dsk()
-        LENGTH_RAY = tool.magnitude(POSITION_FROM_OTHER) * 2.0
-        ORIENTATION_FROM_OTHER = spice.pxform(self.frame, OTHER.frame, EPHEMERIS_TIME)
+        HANDLE_DSK_OTHER = other.handle_dsk()
+        DLADSC_DSK_OTHER = other.dladsc_dsk()
+        LENGTH_RAY = tool.magnitude(STATE_FROM_OTHER_FRAME[:3]) * 2.0
+        ORIENTATION_FROM_OTHER = self.orientation_from(other)
 
         for INDEX, CENTER in enumerate(self.center_facets()):
-            CENTER_FROM_OTHER = ORIENTATION_FROM_OTHER.dot(CENTER) + POSITION_FROM_OTHER
-            RAY = tool.unit_vector(CENTER_FROM_OTHER - POSITION_OBSERVER) * LENGTH_RAY
+            CENTER_FROM_OTHER = ORIENTATION_FROM_OTHER.dot(CENTER) + STATE_FROM_OTHER_FRAME[:3]
+            RAY = tool.unit_vector(CENTER_FROM_OTHER - position_observer) * LENGTH_RAY
             VERTEX = CENTER_FROM_OTHER - RAY
             _, _, FOUND = spice.dskx02(HANDLE_DSK_OTHER, DLADSC_DSK_OTHER, VERTEX, RAY)
             
             if FOUND:
-                shadows[INDEX] = 0.0
+                lights[INDEX] = 0.0
                 
-        return shadows 
+        return lights
+    
+    def compute_self_visibility(self, position_observer: tool.ARRAY_VECTOR) -> tool.ARRAY_LIST:
+        """Spice distances are in kilometers. 1.0 means visible, 0.0 means not.
+        `position_observer` should be defined with respect to the asteroid in it's fixed frame."""
+        # TODO: parallelize facet loop.
+        
+        ABCORR = "LT+S"
+        lights = numpy.ones(self.center_facets().shape[0])
+
+        HANDLE_DSK = self.handle_dsk()
+        DLADSC_DSK = self.dladsc_dsk()
+
+        RADII = spice.bodvrd(self.name, "RADII", 3)[1]
+        RADIUS_MAX = max(RADII)
+        DISTANCE_RAY = RADIUS_MAX * 3.0
+        FACTOR_DISTANCE_RAY = 4.0
+
+        for INDEX, CENTER in enumerate(self.center_facets()):
+            VERTEX = CENTER + tool.unit_vector(position_observer) * DISTANCE_RAY
+            RAY = (CENTER - VERTEX) * FACTOR_DISTANCE_RAY
+            PLID, _, FOUND = spice.dskx02(HANDLE_DSK, DLADSC_DSK, VERTEX, RAY)
+            
+            if FOUND and PLID != INDEX + 1:
+                lights[INDEX] = 0.0
+                
+        return lights
+        
+    def compute_shadows(self, other: Asteroid) -> tool.ARRAY_LIST:
+        STATE_SUN_FROM_OTHER_FRAME = other.state_sun(abcorr="LT+S")
+        return self.compute_visibility(other, STATE_SUN_FROM_OTHER_FRAME[:3])
+
+    def compute_self_shadows(self) -> tool.ARRAY_LIST:
+        STATE_SUN = self.state_sun(abcorr="LT+S")
+        return self.compute_self_visibility(STATE_SUN[:3])
 
     
 class World:
@@ -221,12 +294,12 @@ class World:
         raise ValueError(f"Spacecraft {name} does not exist. {self.string_list_spacecrafts()}")
     
     def add_asteroid(self, name: str, rotation_period: float, frame: str | None = None) -> Asteroid:
-        asteroid = Asteroid(name, rotation_period, frame)
+        asteroid = Asteroid(self, name, rotation_period, frame)
         self.asteroids.append(asteroid)
         return self.asteroids[-1]
         
     def add_spacecraft(self, name: str, frame: str) -> Spacecraft:
-        spacecraft = Spacecraft(name, frame)
+        spacecraft = Spacecraft(self, name, frame)
         self.spacecrafts.append(spacecraft)
         return self.spacecrafts[-1]