From f4e9c3b19d844c033e5757043ff892567ae8cb6c Mon Sep 17 00:00:00 2001
From: sochowski <aidansochowski64@gmail.com>
Date: Thu, 2 Dec 2021 18:06:24 -0600
Subject: [PATCH 01/35] Adding temp PR for spherical rendering

---
 examples/amr_volume_rendering.py              | 21 ++++++++-
 yt_idv/scene_components/base_component.py     |  2 +-
 yt_idv/shaders/ray_tracing.frag.glsl          |  4 ++
 yt_idv/shaders/shaderlist.yaml                |  6 ++-
 .../shaders/spherical_grid_position.vert.glsl | 45 +++++++++++++++++++
 5 files changed, 74 insertions(+), 4 deletions(-)
 create mode 100644 yt_idv/shaders/spherical_grid_position.vert.glsl

diff --git a/examples/amr_volume_rendering.py b/examples/amr_volume_rendering.py
index 17fd8eae..15b4352c 100644
--- a/examples/amr_volume_rendering.py
+++ b/examples/amr_volume_rendering.py
@@ -1,9 +1,26 @@
+import numpy as np
 import yt
 
 import yt_idv
 
-ds = yt.load_sample("IsolatedGalaxy")
+# Spherical Test (to line 20)
+fake_data = {"density": np.random.random((256, 256, 256))}
+ds = yt.load_uniform_grid(
+    fake_data,
+    [256, 256, 256],
+    bbox=np.array([[0.0, 1.0], [0.0, np.pi], [0.0, 2 * np.pi]]),
+    nprocs=4096,
+    geometry="spherical",
+)
 
 rc = yt_idv.render_context(height=800, width=800, gui=True)
-sg = rc.add_scene(ds, "density", no_ghost=True)
+dd = ds.all_data()
+dd.max_level = 1
+sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
 rc.run()
+
+# Cartesian Test (to line 25)
+# ds = yt.load_sample("IsolatedGalaxy")
+# rc = yt_idv.render_context(height=800, width=800, gui=True)
+# sg = rc.add_scene(ds, "density", no_ghost=True)
+# rc.run()
diff --git a/yt_idv/scene_components/base_component.py b/yt_idv/scene_components/base_component.py
index 19d0f8f9..6b2659f3 100644
--- a/yt_idv/scene_components/base_component.py
+++ b/yt_idv/scene_components/base_component.py
@@ -56,7 +56,7 @@ class SceneComponent(traitlets.HasTraits):
     # These attributes are
     cmap_min = traitlets.CFloat(None, allow_none=True)
     cmap_max = traitlets.CFloat(None, allow_none=True)
-    cmap_log = traitlets.Bool(True)
+    cmap_log = traitlets.Bool(False)
     scale = traitlets.CFloat(1.0)
 
     @traitlets.observe("display_bounds")
diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index db3c15e3..b5da51b2 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -32,6 +32,10 @@ vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
 
 void main()
 {
+    // Draws the block outline
+    // output_color = vec4(1.0);
+    // return;
+
     // Obtain screen coordinates
     // https://www.opengl.org/wiki/Compute_eye_space_from_window_space#From_gl_FragCoord
     vec3 ray_position = v_model.xyz;
diff --git a/yt_idv/shaders/shaderlist.yaml b/yt_idv/shaders/shaderlist.yaml
index ee936b42..5967a4d8 100644
--- a/yt_idv/shaders/shaderlist.yaml
+++ b/yt_idv/shaders/shaderlist.yaml
@@ -134,6 +134,9 @@ shader_definitions:
     grid_position:
       info: Pass some grid positions on through
       source: grid_position.vert.glsl
+    spherical_grid_position:
+      info: Pass some spherical grid positions on through and translate them to cartesian 
+      source: spherical_grid_position.vert.glsl
     draw_lines:
       info: A line drawing vertex shader
       source: draw_lines.vert.glsl
@@ -158,7 +161,8 @@ component_shaders:
     default_value: max_intensity
     max_intensity:
       description: Maximum Intensity
-      first_vertex: grid_position
+      first_vertex: spherical_grid_position
+      # first_vertex: grid_position
       first_geometry: grid_expand
       first_fragment: max_intensity
       second_vertex: passthrough
diff --git a/yt_idv/shaders/spherical_grid_position.vert.glsl b/yt_idv/shaders/spherical_grid_position.vert.glsl
new file mode 100644
index 00000000..b297dfcb
--- /dev/null
+++ b/yt_idv/shaders/spherical_grid_position.vert.glsl
@@ -0,0 +1,45 @@
+in vec4 model_vertex; // The location of the vertex in model space
+in vec3 in_dx;
+in vec3 in_left_edge;
+in vec3 in_right_edge;
+flat out vec4 vv_model;
+flat out mat4 vinverse_proj;
+flat out mat4 vinverse_mvm;
+flat out mat4 vinverse_pmvm;
+flat out vec3 vdx;
+flat out vec3 vleft_edge;
+flat out vec3 vright_edge;
+
+vec4 spherical_to_cartesian4(vec4 v) {
+    int theta = 2;
+    int phi = 1;
+    return vec4(
+        v[0] * sin(v[phi]) * cos(v[theta]),
+        v[0] * sin(v[phi]) * sin(v[theta]),
+        v[0] * cos(v[phi]),
+        v[3]
+    );
+}
+
+vec3 spherical_to_cartesian3(vec3 v) {
+    int theta = 2;
+    int phi = 1;
+    return vec3(
+        v[0] * sin(v[phi]) * cos(v[theta]),
+        v[0] * sin(v[phi]) * sin(v[theta]),
+        v[0] * cos(v[phi])
+    );
+}
+
+void main()
+{
+    vv_model = spherical_to_cartesian4(model_vertex);
+    vinverse_proj = inverse(projection);
+    // inverse model-view-matrix
+    vinverse_mvm = inverse(modelview);
+    vinverse_pmvm = inverse(projection * modelview);
+    gl_Position = projection * modelview * vv_model;
+    vdx = in_dx;
+    vleft_edge = spherical_to_cartesian3(in_left_edge);
+    vright_edge = spherical_to_cartesian3(in_right_edge);
+}

From ab78cabe969e455f70d3605038d424ea88846f94 Mon Sep 17 00:00:00 2001
From: sochowski <aidansochowski64@gmail.com>
Date: Wed, 8 Dec 2021 11:17:51 -0600
Subject: [PATCH 02/35] Merged grid_position and spherical_grid_position to
 switch vertex shader functionality based on dataset geometry

---
 yt_idv/scene_components/base_component.py |  3 +++
 yt_idv/shaders/grid_position.vert.glsl    | 26 +++++++++++++++++++----
 yt_idv/shaders/known_uniforms.inc.glsl    |  5 ++++-
 yt_idv/shaders/shaderlist.yaml            |  6 +-----
 4 files changed, 30 insertions(+), 10 deletions(-)

diff --git a/yt_idv/scene_components/base_component.py b/yt_idv/scene_components/base_component.py
index 6b2659f3..a2bde541 100644
--- a/yt_idv/scene_components/base_component.py
+++ b/yt_idv/scene_components/base_component.py
@@ -231,6 +231,9 @@ def run_program(self, scene):
             return
         with self.fb.bind(True):
             with self.program1.enable() as p:
+                p._set_uniform(
+                    "is_spherical", self.data.data_source.ds.geometry == "spherical"
+                )
                 scene.camera._set_uniforms(scene, p)
                 self._set_uniforms(scene, p)
                 p._set_uniform("iso_tolerance", float(self.iso_tolerance))
diff --git a/yt_idv/shaders/grid_position.vert.glsl b/yt_idv/shaders/grid_position.vert.glsl
index 7aa95c7a..bebd3b69 100644
--- a/yt_idv/shaders/grid_position.vert.glsl
+++ b/yt_idv/shaders/grid_position.vert.glsl
@@ -10,16 +10,34 @@ flat out vec3 vdx;
 flat out vec3 vleft_edge;
 flat out vec3 vright_edge;
 
+vec3 transform_vec3(vec3 v) {
+    if (is_spherical) {
+        int theta = 2;
+        int phi = 1;
+        return vec3(
+            v[0] * sin(v[phi]) * cos(v[theta]),
+            v[0] * sin(v[phi]) * sin(v[theta]),
+            v[0] * cos(v[phi])
+        );
+    } else {
+        return v;
+    }
+}
+
+vec4 transform_vec4(vec4 v) {
+    vec3 v3 = transform_vec3(vec3(v));
+    return vec4(v3[0], v3[1], v3[2], v[3]);
+}
+
 void main()
 {
-    vv_model = model_vertex;
+    vv_model = transform_vec4(model_vertex);
     vinverse_proj = inverse(projection);
     // inverse model-view-matrix
     vinverse_mvm = inverse(modelview);
     vinverse_pmvm = inverse(projection * modelview);
     gl_Position = projection * modelview * model_vertex;
     vdx = vec3(in_dx);
-    vleft_edge = vec3(in_left_edge);
-    vright_edge = vec3(in_right_edge);
-
+    vleft_edge = transform_vec3(in_left_edge);
+    vright_edge = transform_vec3(in_right_edge);
 }
diff --git a/yt_idv/shaders/known_uniforms.inc.glsl b/yt_idv/shaders/known_uniforms.inc.glsl
index c1d1ed83..f7bd8e7d 100644
--- a/yt_idv/shaders/known_uniforms.inc.glsl
+++ b/yt_idv/shaders/known_uniforms.inc.glsl
@@ -56,4 +56,7 @@ uniform int iso_num_layers;
 uniform float iso_layers[32];
 uniform bool iso_log;
 uniform float iso_min;
-uniform float iso_max;
\ No newline at end of file
+uniform float iso_max;
+
+// spherical coordinates
+uniform bool is_spherical;
\ No newline at end of file
diff --git a/yt_idv/shaders/shaderlist.yaml b/yt_idv/shaders/shaderlist.yaml
index 5967a4d8..ee936b42 100644
--- a/yt_idv/shaders/shaderlist.yaml
+++ b/yt_idv/shaders/shaderlist.yaml
@@ -134,9 +134,6 @@ shader_definitions:
     grid_position:
       info: Pass some grid positions on through
       source: grid_position.vert.glsl
-    spherical_grid_position:
-      info: Pass some spherical grid positions on through and translate them to cartesian 
-      source: spherical_grid_position.vert.glsl
     draw_lines:
       info: A line drawing vertex shader
       source: draw_lines.vert.glsl
@@ -161,8 +158,7 @@ component_shaders:
     default_value: max_intensity
     max_intensity:
       description: Maximum Intensity
-      first_vertex: spherical_grid_position
-      # first_vertex: grid_position
+      first_vertex: grid_position
       first_geometry: grid_expand
       first_fragment: max_intensity
       second_vertex: passthrough

From af7663dfb296a02fb7c5d31f24ae5038616904ba Mon Sep 17 00:00:00 2001
From: sochowski <aidansochowski64@gmail.com>
Date: Mon, 20 Dec 2021 10:53:08 -0600
Subject: [PATCH 03/35] Temporary PR for spherical rendering 2

---
 yt_idv/shaders/grid_position.vert.glsl | 8 +++++---
 yt_idv/shaders/ray_tracing.frag.glsl   | 2 +-
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/yt_idv/shaders/grid_position.vert.glsl b/yt_idv/shaders/grid_position.vert.glsl
index bebd3b69..06bfbb13 100644
--- a/yt_idv/shaders/grid_position.vert.glsl
+++ b/yt_idv/shaders/grid_position.vert.glsl
@@ -36,8 +36,10 @@ void main()
     // inverse model-view-matrix
     vinverse_mvm = inverse(modelview);
     vinverse_pmvm = inverse(projection * modelview);
-    gl_Position = projection * modelview * model_vertex;
+    gl_Position = projection * modelview * transform_vec4(model_vertex);
     vdx = vec3(in_dx);
-    vleft_edge = transform_vec3(in_left_edge);
-    vright_edge = transform_vec3(in_right_edge);
+    vec3 temp_left_edge = transform_vec3(in_left_edge);
+    vec3 temp_right_edge = transform_vec3(in_right_edge);
+    vleft_edge = min(temp_left_edge, temp_right_edge);
+    vright_edge = max(temp_left_edge, temp_right_edge);
 }
diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index b5da51b2..80d006b3 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -65,7 +65,7 @@ void main()
     temp_t = min(tmax.xx, tmax.yz);
     float t1 = min(temp_t.x, temp_t.y);
     t0 = max(t0, 0.0);
-    if (t1 <= t0) discard;
+    // if (t1 <= t0) discard;
 
     // Some more discussion of this here:
     //  http://prideout.net/blog/?p=64

From ba3ed4812e67d7364025bfa0fc59236ca8a16033 Mon Sep 17 00:00:00 2001
From: Matthew Turk <matthewturk@gmail.com>
Date: Mon, 20 Dec 2021 15:32:31 -0600
Subject: [PATCH 04/35] Update to do transformations in geometry shader

---
 examples/spherical_subset_volume_rendering.py | 29 +++++++++++++++++++
 yt_idv/shaders/grid_expand.geom.glsl          | 16 +++++++++-
 yt_idv/shaders/grid_position.vert.glsl        | 10 +++----
 3 files changed, 48 insertions(+), 7 deletions(-)
 create mode 100644 examples/spherical_subset_volume_rendering.py

diff --git a/examples/spherical_subset_volume_rendering.py b/examples/spherical_subset_volume_rendering.py
new file mode 100644
index 00000000..f8975253
--- /dev/null
+++ b/examples/spherical_subset_volume_rendering.py
@@ -0,0 +1,29 @@
+import numpy as np
+import yt
+
+import yt_idv
+
+# Spherical Test (to line 20)
+
+NDIM = 32
+
+
+bbox = np.array([[0.0, 0.5], [0.0, np.pi / 8], [0.0, 2 * np.pi / 8]])
+
+fake_data = {"density": np.random.random((NDIM, NDIM, NDIM))}
+ds = yt.load_uniform_grid(
+    fake_data, [NDIM, NDIM, NDIM], bbox=bbox, geometry="spherical",
+)
+
+rc = yt_idv.render_context(height=800, width=800, gui=True)
+dd = ds.all_data()
+dd.max_level = 1
+sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
+sg.camera.focus = [0.0, 0.0, 0.0]
+rc.run()
+
+# Cartesian Test (to line 25)
+# ds = yt.load_sample("IsolatedGalaxy")
+# rc = yt_idv.render_context(height=800, width=800, gui=True)
+# sg = rc.add_scene(ds, "density", no_ghost=True)
+# rc.run()
diff --git a/yt_idv/shaders/grid_expand.geom.glsl b/yt_idv/shaders/grid_expand.geom.glsl
index edf0955e..9b63be5c 100644
--- a/yt_idv/shaders/grid_expand.geom.glsl
+++ b/yt_idv/shaders/grid_expand.geom.glsl
@@ -36,6 +36,20 @@ uniform vec3 arrangement[8] = vec3[](
 
 uniform int aindex[14] = int[](6, 7, 4, 5, 1, 7, 3, 6, 2, 4, 0, 1, 2, 3);
 
+vec3 transform_vec3(vec3 v) {
+    if (is_spherical) {
+        int theta = 2;
+        int phi = 1;
+        return vec3(
+            v[0] * sin(v[phi]) * cos(v[theta]),
+            v[0] * sin(v[phi]) * sin(v[theta]),
+            v[0] * cos(v[phi])
+        );
+    } else {
+        return v;
+    }
+}
+
 void main() {
 
     vec4 center = gl_in[0].gl_Position;
@@ -45,7 +59,7 @@ void main() {
     vec4 newPos;
 
     for (int i = 0; i < 14; i++) {
-        newPos = vec4(vleft_edge[0] + width * arrangement[aindex[i]], 1.0);
+        newPos = vec4(transform_vec3(vleft_edge[0] + width * arrangement[aindex[i]]), 1.0);
         gl_Position = projection * modelview * newPos;
         left_edge = vleft_edge[0];
         right_edge = vright_edge[0];
diff --git a/yt_idv/shaders/grid_position.vert.glsl b/yt_idv/shaders/grid_position.vert.glsl
index 06bfbb13..8c18267c 100644
--- a/yt_idv/shaders/grid_position.vert.glsl
+++ b/yt_idv/shaders/grid_position.vert.glsl
@@ -31,15 +31,13 @@ vec4 transform_vec4(vec4 v) {
 
 void main()
 {
-    vv_model = transform_vec4(model_vertex);
+    vv_model = model_vertex;
     vinverse_proj = inverse(projection);
     // inverse model-view-matrix
     vinverse_mvm = inverse(modelview);
     vinverse_pmvm = inverse(projection * modelview);
-    gl_Position = projection * modelview * transform_vec4(model_vertex);
+    gl_Position = projection * modelview * model_vertex;
     vdx = vec3(in_dx);
-    vec3 temp_left_edge = transform_vec3(in_left_edge);
-    vec3 temp_right_edge = transform_vec3(in_right_edge);
-    vleft_edge = min(temp_left_edge, temp_right_edge);
-    vright_edge = max(temp_left_edge, temp_right_edge);
+    vleft_edge = vec3(in_left_edge);
+    vright_edge = vec3(in_right_edge);
 }

From ca6dbfb762f5fb17be16a24850cde54905bbe32f Mon Sep 17 00:00:00 2001
From: Matthew Turk <matthewturk@gmail.com>
Date: Mon, 20 Dec 2021 15:37:28 -0600
Subject: [PATCH 05/35] add constant color shader

---
 examples/spherical_subset_volume_rendering.py |   5 +-
 yt_idv/shaders/shaderlist.yaml                | 109 ++++++++++--------
 2 files changed, 65 insertions(+), 49 deletions(-)

diff --git a/examples/spherical_subset_volume_rendering.py b/examples/spherical_subset_volume_rendering.py
index f8975253..4a939afa 100644
--- a/examples/spherical_subset_volume_rendering.py
+++ b/examples/spherical_subset_volume_rendering.py
@@ -7,8 +7,9 @@
 
 NDIM = 32
 
-
-bbox = np.array([[0.0, 0.5], [0.0, np.pi / 8], [0.0, 2 * np.pi / 8]])
+bbox = np.array(
+    [[0.0, 0.5], [np.pi / 8, 2 * np.pi / 8], [2 * np.pi / 8, 3 * np.pi / 8]]
+)
 
 fake_data = {"density": np.random.random((NDIM, NDIM, NDIM))}
 ds = yt.load_uniform_grid(
diff --git a/yt_idv/shaders/shaderlist.yaml b/yt_idv/shaders/shaderlist.yaml
index ee936b42..53f2d4a8 100644
--- a/yt_idv/shaders/shaderlist.yaml
+++ b/yt_idv/shaders/shaderlist.yaml
@@ -4,124 +4,132 @@ shader_definitions:
       info: A constant value applied
       source: constant.frag.glsl
       blend_func:
-      - one
-      - one
+        - one
+        - one
       blend_equation: func add
     apply_colormap:
-      info: A second pass fragment shader used to apply a colormap to the result of
+      info:
+        A second pass fragment shader used to apply a colormap to the result of
         the first pass rendering
       source: apply_colormap.frag.glsl
       blend_func:
-      - src alpha
-      - dst alpha
+        - src alpha
+        - dst alpha
       blend_equation: func add
     expand_1d:
       info: This expands a 1D texture along the y dimension
       source: expand_1d.frag.glsl
       blend_func:
-      - one
-      - zero
+        - one
+        - zero
       blend_equation: func add
     draw_blocks:
-      info: A first pass fragment shader that performs ray casting using transfer function.
+      info:
+        A first pass fragment shader that performs ray casting using transfer function.
         See :ref:`volume-rendering-method` for more details.
       source: block_outline.frag.glsl
       blend_func:
-      - src alpha
-      - one minus src alpha
+        - src alpha
+        - one minus src alpha
       blend_equation: func add
     isocontour:
       info: A first pass fragment shader that renders isocontour layers.
       source: isocontour.frag.glsl
       blend_func:
-      - src alpha
-      - dst alpha
+        - src alpha
+        - dst alpha
       blend_equation: func add
     max_intensity:
-      info: A first pass fragment shader that computes Maximum Intensity Projection
+      info:
+        A first pass fragment shader that computes Maximum Intensity Projection
         of the data. See :ref:`projection-types` for more information.
       source:
-      - ray_tracing.frag.glsl
-      - max_intensity.frag.glsl
+        - ray_tracing.frag.glsl
+        - max_intensity.frag.glsl
       blend_func:
-      - one
-      - one
+        - one
+        - one
       blend_equation: max
     mesh:
       info: A vertex shader used for unstructured mesh rendering.
       source: mesh.frag.glsl
       depth_test: less
       blend_func:
-      - one
-      - zero
+        - one
+        - zero
       blend_equation: func add
     noop:
-      info: A second pass fragment shader that performs no operation.  Usually used
+      info:
+        A second pass fragment shader that performs no operation.  Usually used
         if the first pass already took care of applying proper color to the data
       source: noop.frag.glsl
     passthrough:
-      info: A first pass fragment shader that performs no operation. Used for debug
+      info:
+        A first pass fragment shader that performs no operation. Used for debug
         puproses. It's distinct from NoOpFragmentShader, because of the number of uniforms
       source: passthrough.frag.glsl
       blend_func:
-      - src alpha
-      - dst alpha
+        - src alpha
+        - dst alpha
       blend_equation: func add
     draw_lines:
       info: A line drawing fragment shader
       source: draw_lines.frag.glsl
       blend_func:
-      - one
-      - zero
+        - one
+        - zero
       blend_equation: func add
     projection:
-      info: A first pass fragment shader that performs unweighted integration of the
+      info:
+        A first pass fragment shader that performs unweighted integration of the
         data along the line of sight. See :ref:`projection-types` for more information.
       source:
-      - ray_tracing.frag.glsl
-      - projection.frag.glsl
+        - ray_tracing.frag.glsl
+        - projection.frag.glsl
       blend_func:
-      - one
-      - one
+        - one
+        - one
       blend_equation: func add
     text_overlay:
       info: A simple text overlay shader
       source: textoverlay.frag.glsl
       blend_func:
-      - src alpha
-      - one minus src alpha
+        - src alpha
+        - one minus src alpha
       blend_equation: func add
     transfer_function:
-      info: A first pass fragment shader that performs ray casting using transfer function.
+      info:
+        A first pass fragment shader that performs ray casting using transfer function.
         See :ref:`volume-rendering-method` for more details.
       source:
-      - ray_tracing.frag.glsl
-      - transfer_function.frag.glsl
+        - ray_tracing.frag.glsl
+        - transfer_function.frag.glsl
       blend_func_separate:
-      - one minus dst alpha
-      - one
-      - one minus dst alpha
-      - one
+        - one minus dst alpha
+        - one
+        - one minus dst alpha
+        - one
       blend_equation_separate:
-      - func add
-      - func add
+        - func add
+        - func add
     sph_kernel:
       info: Sample pre-integrated SPH kernel
       source: sph_kernel.frag.glsl
       blend_func:
-      - one
-      - one
+        - one
+        - one
       blend_equation: func add
     field_value:
       info: Use field values as input
       source: field_value.frag.glsl
       blend_func:
-      - one
-      - one
+        - one
+        - one
       blend_equation: func add
   vertex:
     default:
-      info: A first pass vertex shader that tranlates the location of vertices from
+      info:
+        A first pass vertex shader that tranlates the location of vertices from
         the world coordinates to the viewing plane coordinates
       source: default.vert.glsl
     mesh:
@@ -184,6 +192,13 @@ component_shaders:
       first_fragment: isocontour
       second_vertex: passthrough
       second_fragment: apply_colormap
+    constant:
+      description: Constant Color
+      first_vertex: grid_position
+      first_geometry: grid_expand
+      first_fragment: constant
+      second_vertex: passthrough
+      second_fragment: apply_colormap
   octree_block_rendering:
     default_value: max_intensity
     max_intensity:
@@ -277,4 +292,4 @@ component_shaders:
       first_vertex: particle
       first_fragment: field_value
       second_vertex: passthrough
-      second_fragment: apply_colormap
\ No newline at end of file
+      second_fragment: apply_colormap

From d840d3f80e6f5afb743f2eb9afce1f576645c731 Mon Sep 17 00:00:00 2001
From: sochowski <aidansochowski64@gmail.com>
Date: Thu, 20 Jan 2022 21:51:20 -0600
Subject: [PATCH 06/35] Organized spherical ray-tracing

---
 yt_idv/shaders/ray_tracing.frag.glsl | 44 ++++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 80d006b3..e0220986 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -15,6 +15,20 @@ bool within_bb(vec3 pos)
     return all(left) && all(right);
 }
 
+vec3 transform_vec3(vec3 v) {
+    if (is_spherical) {
+        int theta = 2;
+        int phi = 1;
+        return vec3(
+            v[0] * sin(v[phi]) * cos(v[theta]),
+            v[0] * sin(v[phi]) * sin(v[theta]),
+            v[0] * cos(v[phi])
+        );
+    } else {
+        return v;
+    }
+}
+
 vec3 get_offset_texture_position(sampler3D tex, vec3 tex_curr_pos)
 {
     ivec3 texsize = textureSize(tex, 0); // lod (mipmap level) always 0?
@@ -30,6 +44,31 @@ vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
 //   void (vec3 tex_curr_pos, inout vec4 curr_color, float tdelta, float t,
 //         vec3 direction);
 
+void sph_main()
+{
+    vec3 ray_position = v_model.xyz;
+    vec3 p0 = camera_pos.xyz;
+    vec4 curr_color = vec4(0.0);
+
+    vec3 temp_left_edge = transform_vec3(left_edge);
+    vec3 temp_right_edge = transform_vec3(right_edge);
+    vec3 c_left_edge = min(temp_left_edge, temp_right_edge);
+    vec3 c_right_edge = max(temp_left_edge, temp_right_edge);
+
+    vec3 range = (c_right_edge + dx/2.0) - (c_left_edge - dx/2.0);
+    vec3 nzones = range / dx;
+    vec3 ndx = 1.0/nzones;
+
+    vec3 tex_curr_pos = (ray_position - c_left_edge) / range;
+    tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
+    bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
+    if (sampled) {
+        output_color = curr_color;
+    } else {
+        output_color = vec4(0.0);
+    }
+}
+
 void main()
 {
     // Draws the block outline
@@ -38,6 +77,11 @@ void main()
 
     // Obtain screen coordinates
     // https://www.opengl.org/wiki/Compute_eye_space_from_window_space#From_gl_FragCoord
+    if (is_spherical) {
+        sph_main();
+        return;
+    }
+
     vec3 ray_position = v_model.xyz;
 
     // Five samples

From 4cd83b35f1cfc18465672d6289661a08d466a977 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Tue, 2 Aug 2022 20:11:55 +0000
Subject: [PATCH 07/35] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 examples/spherical_subset_volume_rendering.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/examples/spherical_subset_volume_rendering.py b/examples/spherical_subset_volume_rendering.py
index 4a939afa..682bf724 100644
--- a/examples/spherical_subset_volume_rendering.py
+++ b/examples/spherical_subset_volume_rendering.py
@@ -13,7 +13,10 @@
 
 fake_data = {"density": np.random.random((NDIM, NDIM, NDIM))}
 ds = yt.load_uniform_grid(
-    fake_data, [NDIM, NDIM, NDIM], bbox=bbox, geometry="spherical",
+    fake_data,
+    [NDIM, NDIM, NDIM],
+    bbox=bbox,
+    geometry="spherical",
 )
 
 rc = yt_idv.render_context(height=800, width=800, gui=True)

From 726aa82c844d565fdbc7bee4b5ea5b248346ec9c Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Tue, 2 Aug 2022 15:15:25 -0500
Subject: [PATCH 08/35] spherical example to its own, resetn cartesian

---
 examples/amr_spherical_volume_rendering.py | 19 +++++++++++++++++++
 examples/amr_volume_rendering.py           | 21 ++-------------------
 2 files changed, 21 insertions(+), 19 deletions(-)
 create mode 100644 examples/amr_spherical_volume_rendering.py

diff --git a/examples/amr_spherical_volume_rendering.py b/examples/amr_spherical_volume_rendering.py
new file mode 100644
index 00000000..c84ac5d1
--- /dev/null
+++ b/examples/amr_spherical_volume_rendering.py
@@ -0,0 +1,19 @@
+import numpy as np
+import yt
+
+import yt_idv
+
+fake_data = {"density": np.random.random((256, 256, 256))}
+ds = yt.load_uniform_grid(
+    fake_data,
+    [256, 256, 256],
+    bbox=np.array([[0.0, 1.0], [0.0, np.pi], [0.0, 2 * np.pi]]),
+    nprocs=4096,
+    geometry="spherical",
+)
+
+rc = yt_idv.render_context(height=800, width=800, gui=True)
+dd = ds.all_data()
+dd.max_level = 1
+sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
+rc.run()
diff --git a/examples/amr_volume_rendering.py b/examples/amr_volume_rendering.py
index 15b4352c..17fd8eae 100644
--- a/examples/amr_volume_rendering.py
+++ b/examples/amr_volume_rendering.py
@@ -1,26 +1,9 @@
-import numpy as np
 import yt
 
 import yt_idv
 
-# Spherical Test (to line 20)
-fake_data = {"density": np.random.random((256, 256, 256))}
-ds = yt.load_uniform_grid(
-    fake_data,
-    [256, 256, 256],
-    bbox=np.array([[0.0, 1.0], [0.0, np.pi], [0.0, 2 * np.pi]]),
-    nprocs=4096,
-    geometry="spherical",
-)
+ds = yt.load_sample("IsolatedGalaxy")
 
 rc = yt_idv.render_context(height=800, width=800, gui=True)
-dd = ds.all_data()
-dd.max_level = 1
-sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
+sg = rc.add_scene(ds, "density", no_ghost=True)
 rc.run()
-
-# Cartesian Test (to line 25)
-# ds = yt.load_sample("IsolatedGalaxy")
-# rc = yt_idv.render_context(height=800, width=800, gui=True)
-# sg = rc.add_scene(ds, "density", no_ghost=True)
-# rc.run()

From 3cf8186b799137034b95ef5d182f338c8f6b06ec Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Thu, 1 Dec 2022 14:47:17 -0600
Subject: [PATCH 09/35] experimenting...

---
 yt_idv/opengl_support.py                  |  4 +-
 yt_idv/scene_components/base_component.py |  6 +++
 yt_idv/scene_data/block_collection.py     |  2 -
 yt_idv/shaders/grid_expand.geom.glsl      | 47 +++++++++++++----------
 yt_idv/shaders/grid_position.vert.glsl    | 34 ++++++++--------
 yt_idv/shaders/known_uniforms.inc.glsl    |  3 ++
 yt_idv/shaders/ray_tracing.frag.glsl      | 40 +++++++++----------
 7 files changed, 75 insertions(+), 61 deletions(-)

diff --git a/yt_idv/opengl_support.py b/yt_idv/opengl_support.py
index 316579d6..1d13eaf2 100644
--- a/yt_idv/opengl_support.py
+++ b/yt_idv/opengl_support.py
@@ -353,7 +353,9 @@ def bind(self, program=None):
                 _ = GL.glEnableVertexAttribArray(loc)
         _ = GL.glBindBuffer(GL.GL_ARRAY_BUFFER, self.id)
         if loc >= 0:
-            GL.glVertexAttribPointer(loc, self.each, self.opengl_type, False, 0, None)
+            each_n = self.each
+            print((loc, each_n))
+            GL.glVertexAttribPointer(loc, each_n, self.opengl_type, False, 0, None)
         yield
         if loc >= 0:
             GL.glDisableVertexAttribArray(loc)
diff --git a/yt_idv/scene_components/base_component.py b/yt_idv/scene_components/base_component.py
index 210073df..4b9817a0 100644
--- a/yt_idv/scene_components/base_component.py
+++ b/yt_idv/scene_components/base_component.py
@@ -235,6 +235,12 @@ def run_program(self, scene):
                 p._set_uniform(
                     "is_spherical", self.data.data_source.ds.geometry == "spherical"
                 )
+                if self.data.data_source.ds.geometry == "spherical":
+                    axis_id = self.data.data_source.ds.coordinates.axis_id
+                    p._set_uniform("id_theta", axis_id["theta"])
+                    p._set_uniform("id_r", axis_id["r"])
+                    p._set_uniform("id_phi", axis_id["phi"])
+
                 scene.camera._set_uniforms(scene, p)
                 self._set_uniforms(scene, p)
                 p._set_uniform("iso_tolerance", float(self.iso_tolerance))
diff --git a/yt_idv/scene_data/block_collection.py b/yt_idv/scene_data/block_collection.py
index 4c82be28..49f0c681 100644
--- a/yt_idv/scene_data/block_collection.py
+++ b/yt_idv/scene_data/block_collection.py
@@ -62,8 +62,6 @@ def add_data(self, field, no_ghost=False):
             vert.append([1.0, 1.0, 1.0, 1.0])
             dds = (block.RightEdge - block.LeftEdge) / block.source_mask.shape
             dx.append(dds.tolist())
-            le.append(block.LeftEdge.tolist())
-            re.append(block.RightEdge.tolist())
         for (g, node, (sl, _dims, _gi)) in self.data_source.tiles.slice_traverse():
             block = node.data
             self.blocks_by_grid[g.id - g._id_offset].append((id(block), i))
diff --git a/yt_idv/shaders/grid_expand.geom.glsl b/yt_idv/shaders/grid_expand.geom.glsl
index 9b63be5c..d7fa2bab 100644
--- a/yt_idv/shaders/grid_expand.geom.glsl
+++ b/yt_idv/shaders/grid_expand.geom.glsl
@@ -1,9 +1,9 @@
 layout ( points ) in;
 layout ( triangle_strip, max_vertices = 14 ) out;
 
-flat in vec3 vdx[];
-flat in vec3 vleft_edge[];
-flat in vec3 vright_edge[];
+flat in vec3 vdx[];  // spherical
+flat in vec3 vleft_edge[]; // spherical
+flat in vec3 vright_edge[]; // spherical
 
 flat out vec3 dx;
 flat out vec3 left_edge;
@@ -13,10 +13,10 @@ flat out mat4 inverse_mvm;
 flat out mat4 inverse_pmvm;
 out vec4 v_model;
 
-flat in mat4 vinverse_proj[];
-flat in mat4 vinverse_mvm[];
-flat in mat4 vinverse_pmvm[];
-flat in vec4 vv_model[];
+flat in mat4 vinverse_proj[];  // cartesian
+flat in mat4 vinverse_mvm[]; // cartesian
+flat in mat4 vinverse_pmvm[]; // cartesian
+flat in vec4 vv_model[];  // spherical
 
 flat out ivec3 texture_offset;
 
@@ -38,12 +38,10 @@ uniform int aindex[14] = int[](6, 7, 4, 5, 1, 7, 3, 6, 2, 4, 0, 1, 2, 3);
 
 vec3 transform_vec3(vec3 v) {
     if (is_spherical) {
-        int theta = 2;
-        int phi = 1;
         return vec3(
-            v[0] * sin(v[phi]) * cos(v[theta]),
-            v[0] * sin(v[phi]) * sin(v[theta]),
-            v[0] * cos(v[phi])
+            v[id_r] * sin(v[id_phi]) * cos(v[id_theta]),
+            v[id_r] * sin(v[id_phi]) * sin(v[id_theta]),
+            v[id_r] * cos(v[id_phi])
         );
     } else {
         return v;
@@ -52,24 +50,31 @@ vec3 transform_vec3(vec3 v) {
 
 void main() {
 
-    vec4 center = gl_in[0].gl_Position;
+    vec4 center = gl_in[0].gl_Position;  // cartesian
 
-    vec3 width = vright_edge[0] - vleft_edge[0];
+    vec3 width = vright_edge[0] - vleft_edge[0];  // spherical
 
     vec4 newPos;
+    vec3 current_pos;
 
     for (int i = 0; i < 14; i++) {
-        newPos = vec4(transform_vec3(vleft_edge[0] + width * arrangement[aindex[i]]), 1.0);
-        gl_Position = projection * modelview * newPos;
-        left_edge = vleft_edge[0];
-        right_edge = vright_edge[0];
+        // walks through each vertex of the triangle strip, emits it
+        current_pos = vec3(vleft_edge[0] + width * arrangement[aindex[i]]);
+        newPos = vec4(transform_vec3(current_pos), 1.0);
+        // need to emit gl_positionin cartesian, but pass native coords out in v_model
+        gl_Position = projection * modelview * newPos;  // cartesian
+        left_edge = vleft_edge[0];  // spherical
+        right_edge = vright_edge[0]; // spherical
         inverse_pmvm = vinverse_pmvm[0];
         inverse_proj = vinverse_proj[0];
         inverse_mvm = vinverse_mvm[0];
-        dx = vdx[0];
-        v_model = newPos;
+        dx = vdx[0];  // spherical
+        v_model = vec4(current_pos, 1.0);  // now spherical!
         texture_offset = ivec3(0);
         EmitVertex();
+        // EmitVertex emits the current values of output variables to the
+        // current output primitive on the first (and possibly only) primitive
+        // stream
+        // https://registry.khronos.org/OpenGL-Refpages/gl4/html/EmitVertex.xhtml
     }
-
 }
diff --git a/yt_idv/shaders/grid_position.vert.glsl b/yt_idv/shaders/grid_position.vert.glsl
index 8c18267c..202f24f0 100644
--- a/yt_idv/shaders/grid_position.vert.glsl
+++ b/yt_idv/shaders/grid_position.vert.glsl
@@ -2,22 +2,20 @@ in vec4 model_vertex; // The location of the vertex in model space
 in vec3 in_dx;
 in vec3 in_left_edge;
 in vec3 in_right_edge;
-flat out vec4 vv_model;
-flat out mat4 vinverse_proj;
-flat out mat4 vinverse_mvm;
-flat out mat4 vinverse_pmvm;
-flat out vec3 vdx;
-flat out vec3 vleft_edge;
-flat out vec3 vright_edge;
+flat out vec4 vv_model;  // spherical
+flat out mat4 vinverse_proj;  // cartesian
+flat out mat4 vinverse_mvm; // cartesian
+flat out mat4 vinverse_pmvm; // cartesian
+flat out vec3 vdx; // spherical
+flat out vec3 vleft_edge; // spherical
+flat out vec3 vright_edge; // spherical
 
 vec3 transform_vec3(vec3 v) {
     if (is_spherical) {
-        int theta = 2;
-        int phi = 1;
         return vec3(
-            v[0] * sin(v[phi]) * cos(v[theta]),
-            v[0] * sin(v[phi]) * sin(v[theta]),
-            v[0] * cos(v[phi])
+            v[id_r] * sin(v[id_phi]) * cos(v[id_theta]),
+            v[id_r] * sin(v[id_phi]) * sin(v[id_theta]),
+            v[id_r] * cos(v[id_phi])
         );
     } else {
         return v;
@@ -31,13 +29,15 @@ vec4 transform_vec4(vec4 v) {
 
 void main()
 {
-    vv_model = model_vertex;
+    // camera uniforms:
+    // projection, modelview
+    vv_model = model_vertex; //spherical coords
     vinverse_proj = inverse(projection);
     // inverse model-view-matrix
     vinverse_mvm = inverse(modelview);
     vinverse_pmvm = inverse(projection * modelview);
-    gl_Position = projection * modelview * model_vertex;
-    vdx = vec3(in_dx);
-    vleft_edge = vec3(in_left_edge);
-    vright_edge = vec3(in_right_edge);
+    gl_Position = projection * modelview * transform_vec4(model_vertex);
+    vdx = vec3(in_dx);  // spherical coords
+    vleft_edge = vec3(in_left_edge); // spherical coords
+    vright_edge = vec3(in_right_edge); // spheerical coords
 }
diff --git a/yt_idv/shaders/known_uniforms.inc.glsl b/yt_idv/shaders/known_uniforms.inc.glsl
index 694d7072..5a600603 100644
--- a/yt_idv/shaders/known_uniforms.inc.glsl
+++ b/yt_idv/shaders/known_uniforms.inc.glsl
@@ -63,3 +63,6 @@ uniform float iso_max;
 
 // spherical coordinates
 uniform bool is_spherical;
+uniform int id_theta;
+uniform int id_r;
+uniform int id_phi;
diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index e0220986..709b7f5b 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -1,10 +1,10 @@
-in vec4 v_model;
-flat in vec3 dx;
-flat in vec3 left_edge;
-flat in vec3 right_edge;
-flat in mat4 inverse_proj;
-flat in mat4 inverse_mvm;
-flat in mat4 inverse_pmvm;
+in vec4 v_model; // spherical
+flat in vec3 dx;  // spherical
+flat in vec3 left_edge;  // spherical
+flat in vec3 right_edge;  // spherical
+flat in mat4 inverse_proj;  // cartesian
+flat in mat4 inverse_mvm; // cartesian
+flat in mat4 inverse_pmvm; // cartesian
 flat in ivec3 texture_offset;
 out vec4 output_color;
 
@@ -17,12 +17,10 @@ bool within_bb(vec3 pos)
 
 vec3 transform_vec3(vec3 v) {
     if (is_spherical) {
-        int theta = 2;
-        int phi = 1;
         return vec3(
-            v[0] * sin(v[phi]) * cos(v[theta]),
-            v[0] * sin(v[phi]) * sin(v[theta]),
-            v[0] * cos(v[phi])
+            v[id_r] * sin(v[id_phi]) * cos(v[id_theta]),
+            v[id_r] * sin(v[id_phi]) * sin(v[id_theta]),
+            v[id_r] * cos(v[id_phi])
         );
     } else {
         return v;
@@ -46,20 +44,22 @@ vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
 
 void sph_main()
 {
-    vec3 ray_position = v_model.xyz;
+    vec3 ray_position = v_model.xyz; // now spherical
+
     vec3 p0 = camera_pos.xyz;
     vec4 curr_color = vec4(0.0);
 
-    vec3 temp_left_edge = transform_vec3(left_edge);
-    vec3 temp_right_edge = transform_vec3(right_edge);
-    vec3 c_left_edge = min(temp_left_edge, temp_right_edge);
-    vec3 c_right_edge = max(temp_left_edge, temp_right_edge);
-
-    vec3 range = (c_right_edge + dx/2.0) - (c_left_edge - dx/2.0);
+    vec3 range = (right_edge + dx/2.0) - (left_edge - dx/2.0);
     vec3 nzones = range / dx;
     vec3 ndx = 1.0/nzones;
 
-    vec3 tex_curr_pos = (ray_position - c_left_edge) / range;
+    vec3 tex_curr_pos = (ray_position - left_edge) / range;
+
+    // need a bounds check too
+
+    // should get texture position in spherical space
+    // cartesian_to_spherical(ray position)
+
     tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
     bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
     if (sampled) {

From c9a5bffc6c4cb2035c27532d4ea338736874f417 Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Thu, 1 Dec 2022 15:08:48 -0600
Subject: [PATCH 10/35] working? now only converting to cartesian for GLpos,
 etc.

---
 examples/amr_spherical_volume_rendering.py | 4 +++-
 yt_idv/opengl_support.py                   | 4 +---
 yt_idv/scene_components/base_component.py  | 9 +++++++++
 yt_idv/scene_data/block_collection.py      | 3 +++
 4 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/examples/amr_spherical_volume_rendering.py b/examples/amr_spherical_volume_rendering.py
index c84ac5d1..342be9cd 100644
--- a/examples/amr_spherical_volume_rendering.py
+++ b/examples/amr_spherical_volume_rendering.py
@@ -15,5 +15,7 @@
 rc = yt_idv.render_context(height=800, width=800, gui=True)
 dd = ds.all_data()
 dd.max_level = 1
-sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
+# sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
+sg = rc.add_scene(ds, ("index", "theta"), no_ghost=True)
+# sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
 rc.run()
diff --git a/yt_idv/opengl_support.py b/yt_idv/opengl_support.py
index 1d13eaf2..316579d6 100644
--- a/yt_idv/opengl_support.py
+++ b/yt_idv/opengl_support.py
@@ -353,9 +353,7 @@ def bind(self, program=None):
                 _ = GL.glEnableVertexAttribArray(loc)
         _ = GL.glBindBuffer(GL.GL_ARRAY_BUFFER, self.id)
         if loc >= 0:
-            each_n = self.each
-            print((loc, each_n))
-            GL.glVertexAttribPointer(loc, each_n, self.opengl_type, False, 0, None)
+            GL.glVertexAttribPointer(loc, self.each, self.opengl_type, False, 0, None)
         yield
         if loc >= 0:
             GL.glDisableVertexAttribArray(loc)
diff --git a/yt_idv/scene_components/base_component.py b/yt_idv/scene_components/base_component.py
index 4b9817a0..dcd51bfa 100644
--- a/yt_idv/scene_components/base_component.py
+++ b/yt_idv/scene_components/base_component.py
@@ -240,6 +240,15 @@ def run_program(self, scene):
                     p._set_uniform("id_theta", axis_id["theta"])
                     p._set_uniform("id_r", axis_id["r"])
                     p._set_uniform("id_phi", axis_id["phi"])
+                    # print("setting axis ordering uniforms")
+                    # print({"id_theta": axis_id["theta"],
+                    #        "id_r": axis_id["r"],
+                    #        "id_phi": axis_id["phi"],
+                    #        })
+                else:
+                    p._set_uniform("id_theta", 0)
+                    p._set_uniform("id_r", 0)
+                    p._set_uniform("id_phi", 0)
 
                 scene.camera._set_uniforms(scene, p)
                 self._set_uniforms(scene, p)
diff --git a/yt_idv/scene_data/block_collection.py b/yt_idv/scene_data/block_collection.py
index 49f0c681..79f77a0e 100644
--- a/yt_idv/scene_data/block_collection.py
+++ b/yt_idv/scene_data/block_collection.py
@@ -62,6 +62,9 @@ def add_data(self, field, no_ghost=False):
             vert.append([1.0, 1.0, 1.0, 1.0])
             dds = (block.RightEdge - block.LeftEdge) / block.source_mask.shape
             dx.append(dds.tolist())
+            le.append(block.LeftEdge.tolist())
+            re.append(block.RightEdge.tolist())
+
         for (g, node, (sl, _dims, _gi)) in self.data_source.tiles.slice_traverse():
             block = node.data
             self.blocks_by_grid[g.id - g._id_offset].append((id(block), i))

From 065eb6da0f9b002bd732b421bb20992e31627103 Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Thu, 1 Dec 2022 16:44:35 -0600
Subject: [PATCH 11/35] add a shell fragment to the example, skipped linting

---
 examples/amr_spherical_volume_rendering.py | 50 +++++++++++++++++++---
 yt_idv/shaders/ray_tracing.frag.glsl       | 10 ++---
 2 files changed, 48 insertions(+), 12 deletions(-)

diff --git a/examples/amr_spherical_volume_rendering.py b/examples/amr_spherical_volume_rendering.py
index 342be9cd..5af2aabf 100644
--- a/examples/amr_spherical_volume_rendering.py
+++ b/examples/amr_spherical_volume_rendering.py
@@ -1,21 +1,59 @@
 import numpy as np
 import yt
-
+import unyt
 import yt_idv
 
-fake_data = {"density": np.random.random((256, 256, 256))}
+# yt reminder: phi is the polar angle (0 to 2pi)
+# theta is the angle from north (0 to pi)
+
+bbox = np.array([[0.0, 1.0],  # r
+                 [0.0, np.pi],  # theta
+                 [0.0, 2*np.pi]])  # phi
+sz = (256, 256, 256)
+fake_data = {"density": np.random.random(sz)}
 ds = yt.load_uniform_grid(
     fake_data,
-    [256, 256, 256],
-    bbox=np.array([[0.0, 1.0], [0.0, np.pi], [0.0, 2 * np.pi]]),
+    sz,
+    bbox=bbox,
     nprocs=4096,
-    geometry="spherical",
+    geometry=("spherical", ("r", "theta", "phi")),
+    length_unit="m",
+)
+
+def _shell_fragment(field, data):
+    center_phi = unyt.unyt_quantity(np.pi/2, "")
+    center_theta = unyt.unyt_quantity(np.pi/2, "")
+    center_r = unyt.unyt_quantity(0.5, "m")
+
+    dr = np.abs(data["index", "r"] - center_r)
+    shell = np.exp(-(dr / 0.2))
+    minval = 0.01
+    shell[shell < 0.1] = minval
+
+    # neglecting the periodicity in phi here!
+    phi = data["index", "phi"] # polar 0:2pi angle
+    theta = data["index", "theta"]  # azimuthal  0:pi angle
+    dphi = np.abs(phi - center_phi)
+    dtheta = np.abs(theta - center_theta)
+    dist = 45 * np.pi / 180.
+    shell[dphi > dist / 2.] = 0.05
+    shell[dtheta > dist / 2.] = 0.05
+
+    return shell
+
+
+yt.add_field(
+    name=("stream", "shell_fragment"),
+    function=_shell_fragment,
+    sampling_type="local",
+    units="",
 )
 
 rc = yt_idv.render_context(height=800, width=800, gui=True)
 dd = ds.all_data()
 dd.max_level = 1
 # sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
-sg = rc.add_scene(ds, ("index", "theta"), no_ghost=True)
+# sg = rc.add_scene(ds, ("index", "theta"), no_ghost=True)
 # sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
+sg = rc.add_scene(ds, ("stream", "shell_fragment"), no_ghost=True)
 rc.run()
diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 709b7f5b..8aed1055 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -45,21 +45,19 @@ vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
 void sph_main()
 {
     vec3 ray_position = v_model.xyz; // now spherical
-
+    vec3 ray_position_xyz = transform_vec3(ray_position);
     vec3 p0 = camera_pos.xyz;
+    vec3 dir = -normalize(camera_pos.xyz - ray_position_xyz);
     vec4 curr_color = vec4(0.0);
 
+
+
     vec3 range = (right_edge + dx/2.0) - (left_edge - dx/2.0);
     vec3 nzones = range / dx;
     vec3 ndx = 1.0/nzones;
 
     vec3 tex_curr_pos = (ray_position - left_edge) / range;
 
-    // need a bounds check too
-
-    // should get texture position in spherical space
-    // cartesian_to_spherical(ray position)
-
     tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
     bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
     if (sampled) {

From ddb69fd4694360346d6c1fc311f5c58381e0af64 Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Fri, 2 Dec 2022 10:41:03 -0600
Subject: [PATCH 12/35] swap the angles in the transformation

---
 examples/amr_spherical_volume_rendering.py | 12 ++++++------
 yt_idv/shaders/grid_expand.geom.glsl       |  8 +++++---
 yt_idv/shaders/grid_position.vert.glsl     |  8 +++++---
 yt_idv/shaders/ray_tracing.frag.glsl       |  8 +++++---
 4 files changed, 21 insertions(+), 15 deletions(-)

diff --git a/examples/amr_spherical_volume_rendering.py b/examples/amr_spherical_volume_rendering.py
index 5af2aabf..26672dbe 100644
--- a/examples/amr_spherical_volume_rendering.py
+++ b/examples/amr_spherical_volume_rendering.py
@@ -6,9 +6,9 @@
 # yt reminder: phi is the polar angle (0 to 2pi)
 # theta is the angle from north (0 to pi)
 
-bbox = np.array([[0.0, 1.0],  # r
-                 [0.0, np.pi],  # theta
-                 [0.0, 2*np.pi]])  # phi
+bbox = np.array([[0.3, 1.0],  # r
+                 [0.0, 3*np.pi/2],  # phi
+                 [0.0, np.pi]])  # theta
 sz = (256, 256, 256)
 fake_data = {"density": np.random.random(sz)}
 ds = yt.load_uniform_grid(
@@ -16,7 +16,7 @@
     sz,
     bbox=bbox,
     nprocs=4096,
-    geometry=("spherical", ("r", "theta", "phi")),
+    geometry=("spherical", ("r", "phi", "theta")),
     length_unit="m",
 )
 
@@ -54,6 +54,6 @@ def _shell_fragment(field, data):
 dd.max_level = 1
 # sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
 # sg = rc.add_scene(ds, ("index", "theta"), no_ghost=True)
-# sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
-sg = rc.add_scene(ds, ("stream", "shell_fragment"), no_ghost=True)
+sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
+# sg = rc.add_scene(ds, ("stream", "shell_fragment"), no_ghost=True)
 rc.run()
diff --git a/yt_idv/shaders/grid_expand.geom.glsl b/yt_idv/shaders/grid_expand.geom.glsl
index d7fa2bab..ec5c7821 100644
--- a/yt_idv/shaders/grid_expand.geom.glsl
+++ b/yt_idv/shaders/grid_expand.geom.glsl
@@ -38,10 +38,12 @@ uniform int aindex[14] = int[](6, 7, 4, 5, 1, 7, 3, 6, 2, 4, 0, 1, 2, 3);
 
 vec3 transform_vec3(vec3 v) {
     if (is_spherical) {
+        // yt reminder: phi is the polar angle (0 to 2pi)
+        // theta is the angle from north (0 to pi)
         return vec3(
-            v[id_r] * sin(v[id_phi]) * cos(v[id_theta]),
-            v[id_r] * sin(v[id_phi]) * sin(v[id_theta]),
-            v[id_r] * cos(v[id_phi])
+            v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
+            v[id_r] * sin(v[id_theta]) * sin(v[id_phi]),
+            v[id_r] * cos(v[id_theta])
         );
     } else {
         return v;
diff --git a/yt_idv/shaders/grid_position.vert.glsl b/yt_idv/shaders/grid_position.vert.glsl
index 202f24f0..a72c2485 100644
--- a/yt_idv/shaders/grid_position.vert.glsl
+++ b/yt_idv/shaders/grid_position.vert.glsl
@@ -11,11 +11,13 @@ flat out vec3 vleft_edge; // spherical
 flat out vec3 vright_edge; // spherical
 
 vec3 transform_vec3(vec3 v) {
+    // yt reminder: phi is the polar angle (0 to 2pi)
+    // theta is the angle from north (0 to pi)
     if (is_spherical) {
         return vec3(
-            v[id_r] * sin(v[id_phi]) * cos(v[id_theta]),
-            v[id_r] * sin(v[id_phi]) * sin(v[id_theta]),
-            v[id_r] * cos(v[id_phi])
+            v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
+            v[id_r] * sin(v[id_theta]) * sin(v[id_phi]),
+            v[id_r] * cos(v[id_theta])
         );
     } else {
         return v;
diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 8aed1055..b0154c9e 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -18,9 +18,11 @@ bool within_bb(vec3 pos)
 vec3 transform_vec3(vec3 v) {
     if (is_spherical) {
         return vec3(
-            v[id_r] * sin(v[id_phi]) * cos(v[id_theta]),
-            v[id_r] * sin(v[id_phi]) * sin(v[id_theta]),
-            v[id_r] * cos(v[id_phi])
+        // yt reminder: phi is the polar angle (0 to 2pi)
+        // theta is the angle from north (0 to pi)
+            v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
+            v[id_r] * sin(v[id_theta]) * sin(v[id_phi]),
+            v[id_r] * cos(v[id_theta])
         );
     } else {
         return v;

From 0501c6fb473fdc1350c40cac454db508073acbd1 Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Fri, 2 Dec 2022 10:50:41 -0600
Subject: [PATCH 13/35] simplify the example

---
 examples/amr_spherical_volume_rendering.py | 40 +++-------------------
 1 file changed, 5 insertions(+), 35 deletions(-)

diff --git a/examples/amr_spherical_volume_rendering.py b/examples/amr_spherical_volume_rendering.py
index 26672dbe..02e5e21a 100644
--- a/examples/amr_spherical_volume_rendering.py
+++ b/examples/amr_spherical_volume_rendering.py
@@ -6,9 +6,9 @@
 # yt reminder: phi is the polar angle (0 to 2pi)
 # theta is the angle from north (0 to pi)
 
-bbox = np.array([[0.3, 1.0],  # r
-                 [0.0, 3*np.pi/2],  # phi
-                 [0.0, np.pi]])  # theta
+bbox = np.array([[0.5, 1.0],  # r
+                 [0.0, np.pi/4],  # phi
+                 [np.pi/4, np.pi/2]])  # theta
 sz = (256, 256, 256)
 fake_data = {"density": np.random.random(sz)}
 ds = yt.load_uniform_grid(
@@ -20,40 +20,10 @@
     length_unit="m",
 )
 
-def _shell_fragment(field, data):
-    center_phi = unyt.unyt_quantity(np.pi/2, "")
-    center_theta = unyt.unyt_quantity(np.pi/2, "")
-    center_r = unyt.unyt_quantity(0.5, "m")
-
-    dr = np.abs(data["index", "r"] - center_r)
-    shell = np.exp(-(dr / 0.2))
-    minval = 0.01
-    shell[shell < 0.1] = minval
-
-    # neglecting the periodicity in phi here!
-    phi = data["index", "phi"] # polar 0:2pi angle
-    theta = data["index", "theta"]  # azimuthal  0:pi angle
-    dphi = np.abs(phi - center_phi)
-    dtheta = np.abs(theta - center_theta)
-    dist = 45 * np.pi / 180.
-    shell[dphi > dist / 2.] = 0.05
-    shell[dtheta > dist / 2.] = 0.05
-
-    return shell
-
-
-yt.add_field(
-    name=("stream", "shell_fragment"),
-    function=_shell_fragment,
-    sampling_type="local",
-    units="",
-)
-
 rc = yt_idv.render_context(height=800, width=800, gui=True)
 dd = ds.all_data()
 dd.max_level = 1
-# sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
+sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
 # sg = rc.add_scene(ds, ("index", "theta"), no_ghost=True)
-sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
-# sg = rc.add_scene(ds, ("stream", "shell_fragment"), no_ghost=True)
+# sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
 rc.run()

From 2cb14a854edee933bf08db53cf2c9e3e1d1b8587 Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Fri, 2 Dec 2022 11:30:20 -0600
Subject: [PATCH 14/35] cleanup

---
 yt_idv/scene_components/base_component.py |  9 -----
 yt_idv/shaders/grid_expand.geom.glsl      | 49 +++++++++++------------
 yt_idv/shaders/grid_position.vert.glsl    | 32 ++++++++-------
 yt_idv/shaders/known_uniforms.inc.glsl    |  6 +--
 yt_idv/shaders/ray_tracing.frag.glsl      | 30 +++++++-------
 5 files changed, 60 insertions(+), 66 deletions(-)

diff --git a/yt_idv/scene_components/base_component.py b/yt_idv/scene_components/base_component.py
index dcd51bfa..4b9817a0 100644
--- a/yt_idv/scene_components/base_component.py
+++ b/yt_idv/scene_components/base_component.py
@@ -240,15 +240,6 @@ def run_program(self, scene):
                     p._set_uniform("id_theta", axis_id["theta"])
                     p._set_uniform("id_r", axis_id["r"])
                     p._set_uniform("id_phi", axis_id["phi"])
-                    # print("setting axis ordering uniforms")
-                    # print({"id_theta": axis_id["theta"],
-                    #        "id_r": axis_id["r"],
-                    #        "id_phi": axis_id["phi"],
-                    #        })
-                else:
-                    p._set_uniform("id_theta", 0)
-                    p._set_uniform("id_r", 0)
-                    p._set_uniform("id_phi", 0)
 
                 scene.camera._set_uniforms(scene, p)
                 self._set_uniforms(scene, p)
diff --git a/yt_idv/shaders/grid_expand.geom.glsl b/yt_idv/shaders/grid_expand.geom.glsl
index ec5c7821..f52471cd 100644
--- a/yt_idv/shaders/grid_expand.geom.glsl
+++ b/yt_idv/shaders/grid_expand.geom.glsl
@@ -1,22 +1,24 @@
 layout ( points ) in;
 layout ( triangle_strip, max_vertices = 14 ) out;
 
-flat in vec3 vdx[];  // spherical
-flat in vec3 vleft_edge[]; // spherical
-flat in vec3 vright_edge[]; // spherical
+// note: all in/out variables below are always in native coordinates (e.g.,
+// spherical or cartesian) except when noted.
+flat in vec3 vdx[];
+flat in vec3 vleft_edge[];
+flat in vec3 vright_edge[];
 
 flat out vec3 dx;
 flat out vec3 left_edge;
 flat out vec3 right_edge;
-flat out mat4 inverse_proj;
-flat out mat4 inverse_mvm;
-flat out mat4 inverse_pmvm;
+flat out mat4 inverse_proj; // always cartesian
+flat out mat4 inverse_mvm; // always cartesian
+flat out mat4 inverse_pmvm; // always cartesian
 out vec4 v_model;
 
-flat in mat4 vinverse_proj[];  // cartesian
-flat in mat4 vinverse_mvm[]; // cartesian
-flat in mat4 vinverse_pmvm[]; // cartesian
-flat in vec4 vv_model[];  // spherical
+flat in mat4 vinverse_proj[];  // always cartesian
+flat in mat4 vinverse_mvm[]; // always cartesian
+flat in mat4 vinverse_pmvm[]; // always cartesian
+flat in vec4 vv_model[];
 
 flat out ivec3 texture_offset;
 
@@ -38,8 +40,9 @@ uniform int aindex[14] = int[](6, 7, 4, 5, 1, 7, 3, 6, 2, 4, 0, 1, 2, 3);
 
 vec3 transform_vec3(vec3 v) {
     if (is_spherical) {
-        // yt reminder: phi is the polar angle (0 to 2pi)
-        // theta is the angle from north (0 to pi)
+        // in yt, phi is the polar angle from (0, 2pi), theta is the azimuthal
+        // angle (0, pi). the id_ values below are uniforms that depend on the
+        // yt dataset coordinate ordering
         return vec3(
             v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
             v[id_r] * sin(v[id_theta]) * sin(v[id_phi]),
@@ -52,31 +55,27 @@ vec3 transform_vec3(vec3 v) {
 
 void main() {
 
-    vec4 center = gl_in[0].gl_Position;  // cartesian
+    vec4 center = gl_in[0].gl_Position;  // always cartesian
 
-    vec3 width = vright_edge[0] - vleft_edge[0];  // spherical
+    vec3 width = vright_edge[0] - vleft_edge[0];
 
     vec4 newPos;
     vec3 current_pos;
 
     for (int i = 0; i < 14; i++) {
-        // walks through each vertex of the triangle strip, emits it
+        // walks through each vertex of the triangle strip, emit it. need to
+        // emit gl_Position in cartesian, but pass native coords out in v_model
         current_pos = vec3(vleft_edge[0] + width * arrangement[aindex[i]]);
-        newPos = vec4(transform_vec3(current_pos), 1.0);
-        // need to emit gl_positionin cartesian, but pass native coords out in v_model
+        newPos = vec4(transform_vec3(current_pos), 1.0); // cartesian
         gl_Position = projection * modelview * newPos;  // cartesian
-        left_edge = vleft_edge[0];  // spherical
-        right_edge = vright_edge[0]; // spherical
+        left_edge = vleft_edge[0];
+        right_edge = vright_edge[0];
         inverse_pmvm = vinverse_pmvm[0];
         inverse_proj = vinverse_proj[0];
         inverse_mvm = vinverse_mvm[0];
-        dx = vdx[0];  // spherical
-        v_model = vec4(current_pos, 1.0);  // now spherical!
+        dx = vdx[0];
+        v_model = vec4(current_pos, 1.0);
         texture_offset = ivec3(0);
         EmitVertex();
-        // EmitVertex emits the current values of output variables to the
-        // current output primitive on the first (and possibly only) primitive
-        // stream
-        // https://registry.khronos.org/OpenGL-Refpages/gl4/html/EmitVertex.xhtml
     }
 }
diff --git a/yt_idv/shaders/grid_position.vert.glsl b/yt_idv/shaders/grid_position.vert.glsl
index a72c2485..abf2ee21 100644
--- a/yt_idv/shaders/grid_position.vert.glsl
+++ b/yt_idv/shaders/grid_position.vert.glsl
@@ -1,19 +1,23 @@
-in vec4 model_vertex; // The location of the vertex in model space
+// note: all in/out variables below are always in native coordinates (e.g.,
+// spherical or cartesian) except when noted.
+
+in vec4 model_vertex;
 in vec3 in_dx;
 in vec3 in_left_edge;
 in vec3 in_right_edge;
-flat out vec4 vv_model;  // spherical
-flat out mat4 vinverse_proj;  // cartesian
-flat out mat4 vinverse_mvm; // cartesian
-flat out mat4 vinverse_pmvm; // cartesian
-flat out vec3 vdx; // spherical
-flat out vec3 vleft_edge; // spherical
-flat out vec3 vright_edge; // spherical
+flat out vec4 vv_model;
+flat out mat4 vinverse_proj;  // always cartesian
+flat out mat4 vinverse_mvm;  // always cartesian
+flat out mat4 vinverse_pmvm; // always cartesian
+flat out vec3 vdx;
+flat out vec3 vleft_edge;
+flat out vec3 vright_edge;
 
 vec3 transform_vec3(vec3 v) {
-    // yt reminder: phi is the polar angle (0 to 2pi)
-    // theta is the angle from north (0 to pi)
     if (is_spherical) {
+        // in yt, phi is the polar angle from (0, 2pi), theta is the azimuthal
+        // angle (0, pi). the id_ values below are uniforms that depend on the
+        // yt dataset coordinate ordering
         return vec3(
             v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
             v[id_r] * sin(v[id_theta]) * sin(v[id_phi]),
@@ -33,13 +37,13 @@ void main()
 {
     // camera uniforms:
     // projection, modelview
-    vv_model = model_vertex; //spherical coords
+    vv_model = model_vertex;
     vinverse_proj = inverse(projection);
     // inverse model-view-matrix
     vinverse_mvm = inverse(modelview);
     vinverse_pmvm = inverse(projection * modelview);
     gl_Position = projection * modelview * transform_vec4(model_vertex);
-    vdx = vec3(in_dx);  // spherical coords
-    vleft_edge = vec3(in_left_edge); // spherical coords
-    vright_edge = vec3(in_right_edge); // spheerical coords
+    vdx = vec3(in_dx);
+    vleft_edge = vec3(in_left_edge);
+    vright_edge = vec3(in_right_edge);
 }
diff --git a/yt_idv/shaders/known_uniforms.inc.glsl b/yt_idv/shaders/known_uniforms.inc.glsl
index 5a600603..e5a15b6e 100644
--- a/yt_idv/shaders/known_uniforms.inc.glsl
+++ b/yt_idv/shaders/known_uniforms.inc.glsl
@@ -63,6 +63,6 @@ uniform float iso_max;
 
 // spherical coordinates
 uniform bool is_spherical;
-uniform int id_theta;
-uniform int id_r;
-uniform int id_phi;
+uniform int id_theta;  // azimuthal angle (0 to pi) index in the yt dataset
+uniform int id_r;  // radial index in the yt dataset
+uniform int id_phi;  // polar angle (0 to 2pi) indexi n the yt dataset
diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index b0154c9e..4fea666f 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -1,10 +1,10 @@
-in vec4 v_model; // spherical
-flat in vec3 dx;  // spherical
-flat in vec3 left_edge;  // spherical
-flat in vec3 right_edge;  // spherical
-flat in mat4 inverse_proj;  // cartesian
-flat in mat4 inverse_mvm; // cartesian
-flat in mat4 inverse_pmvm; // cartesian
+in vec4 v_model;
+flat in vec3 dx;
+flat in vec3 left_edge;
+flat in vec3 right_edge;
+flat in mat4 inverse_proj;  // always cartesian
+flat in mat4 inverse_mvm; // always cartesian
+flat in mat4 inverse_pmvm; // always cartesian
 flat in ivec3 texture_offset;
 out vec4 output_color;
 
@@ -18,8 +18,9 @@ bool within_bb(vec3 pos)
 vec3 transform_vec3(vec3 v) {
     if (is_spherical) {
         return vec3(
-        // yt reminder: phi is the polar angle (0 to 2pi)
-        // theta is the angle from north (0 to pi)
+            // in yt, phi is the polar angle from (0, 2pi), theta is the azimuthal
+            // angle (0, pi). the id_ values below are uniforms that depend on the
+            // yt dataset coordinate ordering
             v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
             v[id_r] * sin(v[id_theta]) * sin(v[id_phi]),
             v[id_r] * cos(v[id_theta])
@@ -44,16 +45,15 @@ vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
 //   void (vec3 tex_curr_pos, inout vec4 curr_color, float tdelta, float t,
 //         vec3 direction);
 
-void sph_main()
+void spherical_coord_shortcircuit()
 {
     vec3 ray_position = v_model.xyz; // now spherical
-    vec3 ray_position_xyz = transform_vec3(ray_position);
-    vec3 p0 = camera_pos.xyz;
+    vec3 ray_position_xyz = transform_vec3(ray_position); // cartesian
+    vec3 p0 = camera_pos.xyz; // cartesian
     vec3 dir = -normalize(camera_pos.xyz - ray_position_xyz);
     vec4 curr_color = vec4(0.0);
 
-
-
+    // do the texture evaluation in the native coordinate space
     vec3 range = (right_edge + dx/2.0) - (left_edge - dx/2.0);
     vec3 nzones = range / dx;
     vec3 ndx = 1.0/nzones;
@@ -78,7 +78,7 @@ void main()
     // Obtain screen coordinates
     // https://www.opengl.org/wiki/Compute_eye_space_from_window_space#From_gl_FragCoord
     if (is_spherical) {
-        sph_main();
+        spherical_coord_shortcircuit();
         return;
     }
 

From fb80cea1529336402847c2c08e5a24a1cabbbe3f Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Fri, 2 Dec 2022 11:39:34 -0600
Subject: [PATCH 15/35] bad merge: add back the constant color shader

---
 yt_idv/shaders/shaderlist.yaml | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/yt_idv/shaders/shaderlist.yaml b/yt_idv/shaders/shaderlist.yaml
index f583b1c9..4b0a9cba 100644
--- a/yt_idv/shaders/shaderlist.yaml
+++ b/yt_idv/shaders/shaderlist.yaml
@@ -229,6 +229,13 @@ component_shaders:
       first_fragment: slice_sample
       second_vertex: passthrough
       second_fragment: apply_colormap
+    constant:
+      description: Constant Color
+      first_vertex: grid_position
+      first_geometry: grid_expand
+      first_fragment: constant
+      second_vertex: passthrough
+      second_fragment: apply_colormap
   octree_block_rendering:
     default_value: max_intensity
     max_intensity:

From 65a6abfa1baf5c9d16bca973c3f53bbf1636a51d Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Fri, 2 Dec 2022 11:46:10 -0600
Subject: [PATCH 16/35] more cleaning

---
 examples/amr_spherical_volume_rendering.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/examples/amr_spherical_volume_rendering.py b/examples/amr_spherical_volume_rendering.py
index 02e5e21a..56a3994e 100644
--- a/examples/amr_spherical_volume_rendering.py
+++ b/examples/amr_spherical_volume_rendering.py
@@ -1,14 +1,13 @@
 import numpy as np
 import yt
-import unyt
+
 import yt_idv
 
 # yt reminder: phi is the polar angle (0 to 2pi)
 # theta is the angle from north (0 to pi)
 
-bbox = np.array([[0.5, 1.0],  # r
-                 [0.0, np.pi/4],  # phi
-                 [np.pi/4, np.pi/2]])  # theta
+# coord ordering here will be r, phi, theta
+bbox = np.array([[0.5, 1.0], [0.0, np.pi / 4], [np.pi / 4, np.pi / 2]])
 sz = (256, 256, 256)
 fake_data = {"density": np.random.random(sz)}
 ds = yt.load_uniform_grid(

From f43d5772d12a0e7687e67f20b2a7c71c10b4e86f Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Tue, 6 Dec 2022 16:36:44 -0600
Subject: [PATCH 17/35] move the spherical uniforms to child

---
 yt_idv/scene_components/base_component.py | 9 ---------
 yt_idv/scene_components/blocks.py         | 9 +++++++++
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/yt_idv/scene_components/base_component.py b/yt_idv/scene_components/base_component.py
index d19ad219..563d6b5c 100644
--- a/yt_idv/scene_components/base_component.py
+++ b/yt_idv/scene_components/base_component.py
@@ -245,15 +245,6 @@ def run_program(self, scene):
             return
         with self.fb.bind(True):
             with self.program1.enable() as p:
-                p._set_uniform(
-                    "is_spherical", self.data.data_source.ds.geometry == "spherical"
-                )
-                if self.data.data_source.ds.geometry == "spherical":
-                    axis_id = self.data.data_source.ds.coordinates.axis_id
-                    p._set_uniform("id_theta", axis_id["theta"])
-                    p._set_uniform("id_r", axis_id["r"])
-                    p._set_uniform("id_phi", axis_id["phi"])
-
                 scene.camera._set_uniforms(scene, p)
                 self._set_uniforms(scene, p)
                 p._set_uniform("iso_tolerance", float(self.iso_tolerance))
diff --git a/yt_idv/scene_components/blocks.py b/yt_idv/scene_components/blocks.py
index 72ae4321..6513105e 100644
--- a/yt_idv/scene_components/blocks.py
+++ b/yt_idv/scene_components/blocks.py
@@ -139,6 +139,15 @@ def draw(self, scene, program):
                         GL.glDrawArrays(GL.GL_POINTS, tex_ind * each, each)
 
     def _set_uniforms(self, scene, shader_program):
+        shader_program._set_uniform(
+            "is_spherical", self.data.data_source.ds.geometry == "spherical"
+        )
+        if self.data.data_source.ds.geometry == "spherical":
+            axis_id = self.data.data_source.ds.coordinates.axis_id
+            shader_program._set_uniform("id_theta", axis_id["theta"])
+            shader_program._set_uniform("id_r", axis_id["r"])
+            shader_program._set_uniform("id_phi", axis_id["phi"])
+
         shader_program._set_uniform("box_width", self.box_width)
         shader_program._set_uniform("sample_factor", self.sample_factor)
         shader_program._set_uniform("ds_tex", 0)

From d9a9771b52895b70c611bc197fc19167c62371bc Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Fri, 16 Dec 2022 10:30:25 -0600
Subject: [PATCH 18/35] add a bbox command line arg for example

---
 .gitignore                                 |  2 +
 examples/amr_spherical_volume_rendering.py | 56 +++++++++++++++-------
 2 files changed, 41 insertions(+), 17 deletions(-)

diff --git a/.gitignore b/.gitignore
index c347c913..000fe673 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,5 @@
+# temp files
+*.swp
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]
diff --git a/examples/amr_spherical_volume_rendering.py b/examples/amr_spherical_volume_rendering.py
index 56a3994e..be781183 100644
--- a/examples/amr_spherical_volume_rendering.py
+++ b/examples/amr_spherical_volume_rendering.py
@@ -1,3 +1,5 @@
+import sys
+
 import numpy as np
 import yt
 
@@ -6,23 +8,43 @@
 # yt reminder: phi is the polar angle (0 to 2pi)
 # theta is the angle from north (0 to pi)
 
+
 # coord ordering here will be r, phi, theta
-bbox = np.array([[0.5, 1.0], [0.0, np.pi / 4], [np.pi / 4, np.pi / 2]])
+
+bbox_options = {
+    "partial": np.array([[0.5, 1.0], [0.0, np.pi / 4], [np.pi / 4, np.pi / 2]]),
+    "whole": np.array([[0.1, 1.0], [0.0, 2 * np.pi], [0, np.pi]]),
+    "north_hemi": np.array([[0.1, 1.0], [0.0, 2 * np.pi], [0, 0.5 * np.pi]]),
+    "south_hemi": np.array([[0.1, 1.0], [0.0, 2 * np.pi], [0.5 * np.pi, np.pi]]),
+    "ew_hemi": np.array([[0.1, 1.0], [0.0, np.pi], [0.0, np.pi]]),
+}
+
+
 sz = (256, 256, 256)
 fake_data = {"density": np.random.random(sz)}
-ds = yt.load_uniform_grid(
-    fake_data,
-    sz,
-    bbox=bbox,
-    nprocs=4096,
-    geometry=("spherical", ("r", "phi", "theta")),
-    length_unit="m",
-)
-
-rc = yt_idv.render_context(height=800, width=800, gui=True)
-dd = ds.all_data()
-dd.max_level = 1
-sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
-# sg = rc.add_scene(ds, ("index", "theta"), no_ghost=True)
-# sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
-rc.run()
+
+if __name__ == "__main__":
+    if len(sys.argv) > 1:
+        bbox_type = sys.argv[1]
+    else:
+        bbox_type = "partial"
+
+    bbox = bbox_options[bbox_type]
+
+    ds = yt.load_uniform_grid(
+        fake_data,
+        sz,
+        bbox=bbox,
+        nprocs=4096,
+        geometry=("spherical", ("r", "phi", "theta")),
+        length_unit="m",
+    )
+
+    rc = yt_idv.render_context(height=800, width=800, gui=True)
+    dd = ds.all_data()
+    dd.max_level = 1
+    sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
+    # sg = rc.add_scene(ds, ("index", "theta"), no_ghost=True)
+    # sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
+    sg.camera.focus = [0.0, 0.0, 0.0]
+    rc.run()

From 0474454f482b13661f3f13ad5d1aff31ff645d86 Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Fri, 6 Jan 2023 14:13:45 -0600
Subject: [PATCH 19/35] in progress: spherical ray intersections

---
 yt_idv/scene_data/block_collection.py  |  48 ++++++++++
 yt_idv/shaders/grid_expand.geom.glsl   |   7 ++
 yt_idv/shaders/grid_position.vert.glsl |  10 ++
 yt_idv/shaders/ray_tracing.frag.glsl   | 121 +++++++++++++++++++++++++
 4 files changed, 186 insertions(+)

diff --git a/yt_idv/scene_data/block_collection.py b/yt_idv/scene_data/block_collection.py
index 79f77a0e..1498f114 100644
--- a/yt_idv/scene_data/block_collection.py
+++ b/yt_idv/scene_data/block_collection.py
@@ -41,6 +41,39 @@ def add_data(self, field, no_ghost=False):
         # Every time we change our data source, we wipe all existing ones.
         # We now set up our vertices into our current data source.
         vert, dx, le, re = [], [], [], []
+
+        # spherical-only: SHOULD ONLY DO THIS IF DATASET IS SPHERICAL
+        # normal vectors and plane constants for min/max phi face
+        phi_plane_le = []
+        phi_plane_re = []
+
+        axis_id = self.data_source.ds.coordinates.axis_id
+
+        def calculate_phi_normal_plane(le_or_re):
+            # calculates the parameters for a plane parallel to the z-axis and
+            # passing through the x-y values defined by the polar angle phi
+            # the result plane defines the boundary on +/-phi faces.
+            #
+            # eq of plane:
+            # X dot N = d where N is a normal vector, d is a constant,
+            #             X is a position vector
+            # this function returns N and d.
+            phi = le_or_re[axis_id["phi"]]
+            theta = le_or_re[axis_id["theta"]]
+            r = le_or_re[axis_id["r"]]
+
+            # prob can/should use a yt func here to ensure consistency....
+            z = r * np.cos(theta)
+            r_xy = r * np.sin(theta)
+            x = r_xy * np.cos(phi)
+            y = r_xy * np.sin(phi)
+            pt = np.array([x, y, z])
+
+            z_hat = np.array([0, 0, 1])
+            normal_vec = np.cross(pt, z_hat)
+            d = np.dot(normal_vec, pt)
+            return normal_vec, d
+
         self.min_val = +np.inf
         self.max_val = -np.inf
         if self.scale:
@@ -65,6 +98,12 @@ def add_data(self, field, no_ghost=False):
             le.append(block.LeftEdge.tolist())
             re.append(block.RightEdge.tolist())
 
+            normal, const = calculate_phi_normal_plane(le[-1])
+            phi_plane_le.append(np.array([*normal, const]))
+
+            normal, const = calculate_phi_normal_plane(re[-1])
+            phi_plane_re.append(np.array([*normal, const]))
+
         for (g, node, (sl, _dims, _gi)) in self.data_source.tiles.slice_traverse():
             block = node.data
             self.blocks_by_grid[g.id - g._id_offset].append((id(block), i))
@@ -95,6 +134,15 @@ def add_data(self, field, no_ghost=False):
             VertexAttribute(name="in_right_edge", data=re)
         )
 
+        phi_plane_re = np.array(phi_plane_re, dtype="f4")
+        phi_plane_le = np.array(phi_plane_le, dtype="f4")
+        self.vertex_array.attributes.append(
+            VertexAttribute(name="phi_plane_le", data=phi_plane_le)
+        )
+        self.vertex_array.attributes.append(
+            VertexAttribute(name="phi_plane_re", data=phi_plane_re)
+        )
+
         # Now we set up our textures
         self._load_textures()
 
diff --git a/yt_idv/shaders/grid_expand.geom.glsl b/yt_idv/shaders/grid_expand.geom.glsl
index f52471cd..0a2a24f6 100644
--- a/yt_idv/shaders/grid_expand.geom.glsl
+++ b/yt_idv/shaders/grid_expand.geom.glsl
@@ -6,6 +6,8 @@ layout ( triangle_strip, max_vertices = 14 ) out;
 flat in vec3 vdx[];
 flat in vec3 vleft_edge[];
 flat in vec3 vright_edge[];
+flat in vec4 vphi_plane_le[];
+flat in vec4 vphi_plane_re[];
 
 flat out vec3 dx;
 flat out vec3 left_edge;
@@ -22,6 +24,9 @@ flat in vec4 vv_model[];
 
 flat out ivec3 texture_offset;
 
+flat out vec4 phi_plane_le;
+flat out vec4 phi_plane_re;
+
 // https://stackoverflow.com/questions/28375338/cube-using-single-gl-triangle-strip
 // suggests that the triangle strip we want for the cube is
 
@@ -73,6 +78,8 @@ void main() {
         inverse_pmvm = vinverse_pmvm[0];
         inverse_proj = vinverse_proj[0];
         inverse_mvm = vinverse_mvm[0];
+        phi_plane_le = vphi_plane_le[0];
+        phi_plane_re = vphi_plane_re[0];
         dx = vdx[0];
         v_model = vec4(current_pos, 1.0);
         texture_offset = ivec3(0);
diff --git a/yt_idv/shaders/grid_position.vert.glsl b/yt_idv/shaders/grid_position.vert.glsl
index abf2ee21..d61ed56c 100644
--- a/yt_idv/shaders/grid_position.vert.glsl
+++ b/yt_idv/shaders/grid_position.vert.glsl
@@ -5,6 +5,10 @@ in vec4 model_vertex;
 in vec3 in_dx;
 in vec3 in_left_edge;
 in vec3 in_right_edge;
+
+in vec4 phi_plane_le;  // first 3 elements are the normal vector, final is constant
+in vec4 phi_plane_re;
+
 flat out vec4 vv_model;
 flat out mat4 vinverse_proj;  // always cartesian
 flat out mat4 vinverse_mvm;  // always cartesian
@@ -13,6 +17,9 @@ flat out vec3 vdx;
 flat out vec3 vleft_edge;
 flat out vec3 vright_edge;
 
+flat out vec4 vphi_plane_le;
+flat out vec4 vphi_plane_re;
+
 vec3 transform_vec3(vec3 v) {
     if (is_spherical) {
         // in yt, phi is the polar angle from (0, 2pi), theta is the azimuthal
@@ -46,4 +53,7 @@ void main()
     vdx = vec3(in_dx);
     vleft_edge = vec3(in_left_edge);
     vright_edge = vec3(in_right_edge);
+
+    vphi_plane_le = vec4(phi_plane_le);
+    vphi_plane_re = vec4(phi_plane_re);
 }
diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 4fea666f..1f9311a6 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -8,6 +8,11 @@ flat in mat4 inverse_pmvm; // always cartesian
 flat in ivec3 texture_offset;
 out vec4 output_color;
 
+flat in vec4 phi_plane_le;
+flat in vec4 phi_plane_re;
+
+const float INFINITY = 1. / 0.;
+const float PI = 3.1415926535897932384626433832795;
 bool within_bb(vec3 pos)
 {
     bvec3 left =  greaterThanEqual(pos, left_edge);
@@ -45,6 +50,106 @@ vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
 //   void (vec3 tex_curr_pos, inout vec4 curr_color, float tdelta, float t,
 //         vec3 direction);
 
+float get_ray_plane_intersection(vec3 p_normal, float p_constant, vec3 ray_origin, vec3 ray_dir)
+{
+    float n_dot_u = dot(p_normal, ray_dir);
+    float n_dot_ro = dot(p_normal, ray_origin);
+    // add check for n_dot_u == 0 (ray is parallel to plane)
+    if (n_dot_u == float(0.0))
+    {
+        // the ray is parallel to the plane. there are either no intersections
+        // or infinite intersections. In the second case, we are guaranteed
+        // to intersect one of the other faces, so we can drop this plane.
+        return INFINITY;
+    }
+
+    return (p_constant - n_dot_ro) / n_dot_u;
+}
+
+vec2 get_ray_sphere_intersection(float r, vec3 ray_origin, vec3 ray_dir)
+{
+    float dir_dot_dir = dot(ray_dir, ray_dir);
+    float ro_dot_ro = dot(ray_origin, ray_origin);
+    float dir_dot_ro = dot(ray_dir, ray_origin);
+    float rsq = r * r; // could be calculated in vertex shader
+
+    float a_2 = 2.0 * dir_dot_dir;
+    float b = 2.0 * dir_dot_ro;
+    float c =  ro_dot_ro - rsq;
+    float determinate = b*b - 2.0 * a_2 * c;
+    float cutoff_val = 0.0;
+    if (determinate < cutoff_val)
+    {
+        return vec2(INFINITY, INFINITY);
+    }
+    else if (determinate == cutoff_val)
+    {
+        return vec2(-b / a_2, INFINITY);
+    }
+    else
+    {
+        return vec2((-b - sqrt(determinate))/ a_2, (-b + sqrt(determinate))/ a_2);
+    }
+
+}
+
+vec2 get_ray_cone_intersection(float theta, vec3 ray_origin, vec3 ray_dir)
+{
+    // note: cos(theta) and vhat could be calculated in vertex shader
+    float costheta;
+    vec3 vhat;
+    if (theta > PI/2.0)
+    {
+        // if theta is past PI/2, the cone will point in negative z and the
+        // half angle should be measured from the -z axis, not +z.
+        vhat = vec3(0.0, 0.0, -1.0);
+        costheta = cos(PI - theta);
+    }
+    else
+    {
+        vhat = vec3(0.0, 0.0, 1.0);
+        costheta = cos(theta);
+    }
+    float cos2t = pow(costheta, 2);
+    // note: theta = PI/2.0 is well defined. determinate = 0 in that case and
+    // the cone becomes a plane in x-y.
+
+    float dir_dot_vhat = dot(ray_dir, vhat);
+    float dir_dot_dir = dot(ray_dir, ray_dir);
+    float ro_dot_vhat = dot(ray_origin, vhat);
+    float ro_dot_dir = dot(ray_origin, ray_dir);
+    float ro_dot_ro = dot(ray_origin, ray_dir);
+
+    float a_2 = 2.0*(pow(dir_dot_vhat, 2) - dir_dot_dir * cos2t);
+    float b = 2.0 * (ro_dot_vhat * dir_dot_vhat - ro_dot_dir*cos2t);
+    float c = pow(ro_dot_vhat, 2) - ro_dot_ro*cos2t;
+    float determinate = b*b - 2.0 * a_2 * c;
+    if (determinate < 0.0)
+    {
+        return vec2(INFINITY, INFINITY);
+    }
+    else if (determinate == 0.0)
+    {
+        return vec2(-b / a_2, INFINITY);
+    }
+    else
+    {
+        vec2 t = vec2((-b - sqrt(determinate))/ a_2, (-b + sqrt(determinate))/ a_2);
+
+        // check that t values are not for the shadow cone
+        float cutoff_t = (-1.0 * ro_dot_vhat) / dir_dot_vhat;
+        if (t[0] < cutoff_t)
+        {
+            t[0] = INFINITY;
+        }
+        else if (t[1] < cutoff_t)
+        {
+            t[1] = INFINITY;
+        }
+        return t;
+    }
+}
+
 void spherical_coord_shortcircuit()
 {
     vec3 ray_position = v_model.xyz; // now spherical
@@ -53,6 +158,22 @@ void spherical_coord_shortcircuit()
     vec3 dir = -normalize(camera_pos.xyz - ray_position_xyz);
     vec4 curr_color = vec4(0.0);
 
+    // intersections
+    vec2 t_sphere_outer = get_ray_sphere_intersection(right_edge[id_r], ray_position_xyz, dir);
+    if (t_sphere_outer[0] == INFINITY && t_sphere_outer[1] == INFINITY)
+    {
+        // if there are no intersections with the outer sphere, then there
+        // will be no intersections with the remaining faces and we can stop
+        // looking.
+        discard;
+    }
+
+    vec2 t_sphere_inner = get_ray_sphere_intersection(left_edge[id_r], ray_position_xyz, dir);
+    float t_p_1 = get_ray_plane_intersection(vec3(phi_plane_le), phi_plane_le[3], ray_position_xyz, dir);
+    float t_p_2 = get_ray_plane_intersection(vec3(phi_plane_re), phi_plane_re[3], ray_position_xyz, dir);
+    vec2 t_cone_outer = get_ray_cone_intersection(right_edge[id_theta], ray_position_xyz, dir);
+    vec2 t_cone_inner= get_ray_cone_intersection(left_edge[id_theta], ray_position_xyz, dir);
+
     // do the texture evaluation in the native coordinate space
     vec3 range = (right_edge + dx/2.0) - (left_edge - dx/2.0);
     vec3 nzones = range / dx;

From 0031ee881dffb8474ab9216f5bf87cab748c0192 Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Mon, 9 Jan 2023 12:13:55 -0600
Subject: [PATCH 20/35] vectorize normal plane calc, only do it if spherical

---
 yt_idv/scene_data/_geometry_utils.py  | 37 ++++++++++++
 yt_idv/scene_data/block_collection.py | 81 ++++++++++-----------------
 2 files changed, 67 insertions(+), 51 deletions(-)
 create mode 100644 yt_idv/scene_data/_geometry_utils.py

diff --git a/yt_idv/scene_data/_geometry_utils.py b/yt_idv/scene_data/_geometry_utils.py
new file mode 100644
index 00000000..a646f6fc
--- /dev/null
+++ b/yt_idv/scene_data/_geometry_utils.py
@@ -0,0 +1,37 @@
+import numpy as np
+
+
+def phi_normal_planes(edge_coordinates, axis_id, cast_type: str = None):
+    # for spherical geometries, calculates the cartesian normals and constants
+    # defining the planes normal to a fixed-phi value. The phi normal plane for
+    # a given spherical coordinate (r, theta, phi) will contain the given
+    # coordinate and the z-axis.
+    #
+    # edge_coordinates: 3D array of shape (N, 3) containing the spherical
+    #                   coordinates for which we want the phi-normal.
+    # axis_id: dictionary that maps the spherical coordinate axis names to the
+    #          index number.
+
+    phi = edge_coordinates[:, axis_id["phi"]]
+    theta = edge_coordinates[:, axis_id["theta"]]
+    r = edge_coordinates[:, axis_id["r"]]
+
+    # get the cartesian values of the coordinates
+    z = r * np.cos(theta)
+    r_xy = r * np.sin(theta)
+    x = r_xy * np.cos(phi)
+    y = r_xy * np.sin(phi)
+    xyz = np.column_stack([x, y, z])
+
+    # construct the planes
+    z_hat = np.array([0, 0, 1])
+    # cross product is vectorized, result is shape (N, 3):
+    normal_vec = np.cross(xyz, z_hat)
+    # dot product is not vectorized, do it manually via an elemntwise multiplication
+    # then summation. result will have shape (N,)
+    d = (normal_vec * xyz).sum(axis=1)  # manual dot product
+
+    normals_d = np.column_stack([normal_vec, d])
+    if cast_type is not None:
+        normals_d = normals_d.astype(cast_type)
+    return normals_d
diff --git a/yt_idv/scene_data/block_collection.py b/yt_idv/scene_data/block_collection.py
index 1498f114..19e3987b 100644
--- a/yt_idv/scene_data/block_collection.py
+++ b/yt_idv/scene_data/block_collection.py
@@ -5,6 +5,7 @@
 from yt.data_objects.data_containers import YTDataContainer
 
 from yt_idv.opengl_support import Texture3D, VertexArray, VertexAttribute
+from yt_idv.scene_data._geometry_utils import phi_normal_planes
 from yt_idv.scene_data.base_data import SceneData
 
 
@@ -42,38 +43,6 @@ def add_data(self, field, no_ghost=False):
         # We now set up our vertices into our current data source.
         vert, dx, le, re = [], [], [], []
 
-        # spherical-only: SHOULD ONLY DO THIS IF DATASET IS SPHERICAL
-        # normal vectors and plane constants for min/max phi face
-        phi_plane_le = []
-        phi_plane_re = []
-
-        axis_id = self.data_source.ds.coordinates.axis_id
-
-        def calculate_phi_normal_plane(le_or_re):
-            # calculates the parameters for a plane parallel to the z-axis and
-            # passing through the x-y values defined by the polar angle phi
-            # the result plane defines the boundary on +/-phi faces.
-            #
-            # eq of plane:
-            # X dot N = d where N is a normal vector, d is a constant,
-            #             X is a position vector
-            # this function returns N and d.
-            phi = le_or_re[axis_id["phi"]]
-            theta = le_or_re[axis_id["theta"]]
-            r = le_or_re[axis_id["r"]]
-
-            # prob can/should use a yt func here to ensure consistency....
-            z = r * np.cos(theta)
-            r_xy = r * np.sin(theta)
-            x = r_xy * np.cos(phi)
-            y = r_xy * np.sin(phi)
-            pt = np.array([x, y, z])
-
-            z_hat = np.array([0, 0, 1])
-            normal_vec = np.cross(pt, z_hat)
-            d = np.dot(normal_vec, pt)
-            return normal_vec, d
-
         self.min_val = +np.inf
         self.max_val = -np.inf
         if self.scale:
@@ -98,12 +67,6 @@ def calculate_phi_normal_plane(le_or_re):
             le.append(block.LeftEdge.tolist())
             re.append(block.RightEdge.tolist())
 
-            normal, const = calculate_phi_normal_plane(le[-1])
-            phi_plane_le.append(np.array([*normal, const]))
-
-            normal, const = calculate_phi_normal_plane(re[-1])
-            phi_plane_re.append(np.array([*normal, const]))
-
         for (g, node, (sl, _dims, _gi)) in self.data_source.tiles.slice_traverse():
             block = node.data
             self.blocks_by_grid[g.id - g._id_offset].append((id(block), i))
@@ -117,35 +80,51 @@ def calculate_phi_normal_plane(le_or_re):
         LE = np.array([b.LeftEdge for i, b in self.blocks.values()]).min(axis=0)
         RE = np.array([b.RightEdge for i, b in self.blocks.values()]).max(axis=0)
         self.diagonal = np.sqrt(((RE - LE) ** 2).sum())
+
         # Now we set up our buffer
         vert = np.array(vert, dtype="f4")
         dx = np.array(dx, dtype="f4")
-        le = np.array(le, dtype="f4")
-        re = np.array(re, dtype="f4")
+        le = np.array(le)  # delay type conversion until after _set_geometry_attributes
+        re = np.array(re)
+        self._set_geometry_attributes(le, re)
 
         self.vertex_array.attributes.append(
             VertexAttribute(name="model_vertex", data=vert)
         )
         self.vertex_array.attributes.append(VertexAttribute(name="in_dx", data=dx))
         self.vertex_array.attributes.append(
-            VertexAttribute(name="in_left_edge", data=le)
+            VertexAttribute(name="in_left_edge", data=le.astype("f4"))
         )
         self.vertex_array.attributes.append(
-            VertexAttribute(name="in_right_edge", data=re)
-        )
-
-        phi_plane_re = np.array(phi_plane_re, dtype="f4")
-        phi_plane_le = np.array(phi_plane_le, dtype="f4")
-        self.vertex_array.attributes.append(
-            VertexAttribute(name="phi_plane_le", data=phi_plane_le)
-        )
-        self.vertex_array.attributes.append(
-            VertexAttribute(name="phi_plane_re", data=phi_plane_re)
+            VertexAttribute(name="in_right_edge", data=re.astype("f4"))
         )
 
         # Now we set up our textures
         self._load_textures()
 
+    def _set_geometry_attributes(self, le, re):
+        # set any vertex_array attributes that depend on the yt geometry type
+
+        # note: these comparisons should change to literal comparisons to enum
+        # members with future yt, see https://github.com/yt-project/yt/pull/4244
+        yt_geom = self.data_source.ds.geometry
+        if yt_geom == "cartesian":
+            return
+        elif yt_geom == "spherical":
+            axis_id = self.data_source.ds.coordinates.axis_id
+            phi_plane_le = phi_normal_planes(le, axis_id, cast_type="f4")
+            phi_plane_re = phi_normal_planes(re, axis_id, cast_type="f4")
+            self.vertex_array.attributes.append(
+                VertexAttribute(name="phi_plane_le", data=phi_plane_le)
+            )
+            self.vertex_array.attributes.append(
+                VertexAttribute(name="phi_plane_re", data=phi_plane_re)
+            )
+        else:
+            raise NotImplementedError(
+                f"{self.name} does not implement {yt_geom} geometries."
+            )
+
     def viewpoint_iter(self, camera):
         for block in self.data_source.tiles.traverse(viewpoint=camera.position):
             vbo_i, _ = self.blocks[id(block)]

From 4e4fdecc3b628fbd40e251c694266188a05921e5 Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Wed, 12 Apr 2023 11:27:28 -0500
Subject: [PATCH 21/35] lingering changes, they look OK

---
 yt_idv/shaders/ray_tracing.frag.glsl | 20 ++++++--------------
 1 file changed, 6 insertions(+), 14 deletions(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 1f9311a6..d536df76 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -134,19 +134,11 @@ vec2 get_ray_cone_intersection(float theta, vec3 ray_origin, vec3 ray_dir)
     }
     else
     {
-        vec2 t = vec2((-b - sqrt(determinate))/ a_2, (-b + sqrt(determinate))/ a_2);
-
-        // check that t values are not for the shadow cone
-        float cutoff_t = (-1.0 * ro_dot_vhat) / dir_dot_vhat;
-        if (t[0] < cutoff_t)
-        {
-            t[0] = INFINITY;
-        }
-        else if (t[1] < cutoff_t)
-        {
-            t[1] = INFINITY;
-        }
-        return t;
+        // note: it is possible to have real solutions that intersect the shadow cone
+        // and not the actual cone. those values should be discarded. But they will
+        // fail subsequent bounds checks for interesecting the volume, so we can
+        // just handle it there instead of adding another check here.
+        return vec2((-b - sqrt(determinate))/ a_2, (-b + sqrt(determinate))/ a_2);
     }
 }
 
@@ -160,7 +152,7 @@ void spherical_coord_shortcircuit()
 
     // intersections
     vec2 t_sphere_outer = get_ray_sphere_intersection(right_edge[id_r], ray_position_xyz, dir);
-    if (t_sphere_outer[0] == INFINITY && t_sphere_outer[1] == INFINITY)
+    if (isinf(t_sphere_outer[0]) && isinf(t_sphere_outer[1]))
     {
         // if there are no intersections with the outer sphere, then there
         // will be no intersections with the remaining faces and we can stop

From afcef53ac0f2b2d1501f1aee811e8a005826681d Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Tue, 6 Jun 2023 15:16:59 -0500
Subject: [PATCH 22/35] add a simple test

---
 tests/test_spherical_vol_rendering.py | 44 +++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)
 create mode 100644 tests/test_spherical_vol_rendering.py

diff --git a/tests/test_spherical_vol_rendering.py b/tests/test_spherical_vol_rendering.py
new file mode 100644
index 00000000..69495fad
--- /dev/null
+++ b/tests/test_spherical_vol_rendering.py
@@ -0,0 +1,44 @@
+import os
+
+import numpy as np
+import pytest
+import yt
+
+import yt_idv
+
+
+@pytest.fixture()
+def osmesa_fake_spherical():
+    """Return an OSMesa context that has a "fake" AMR dataset added, with "radius"
+    as the field.
+    """
+
+    sz = (10, 10, 10)
+    fake_data = {"density": np.random.random(sz)}
+
+    bbox = np.array([[0.1, 1.0], [0.0, 2 * np.pi], [0, np.pi]])
+
+    ds = yt.load_uniform_grid(
+        fake_data,
+        sz,
+        bbox=bbox,
+        nprocs=1,
+        geometry=("spherical", ("r", "phi", "theta")),
+        length_unit="m",
+    )
+    dd = ds.all_data()
+    rc = yt_idv.render_context("osmesa", width=1024, height=1024)
+    rc.add_scene(dd, "density", no_ghost=True)
+    yield rc
+    rc.osmesa.OSMesaDestroyContext(rc.context)
+
+
+def test_spherical(osmesa_fake_spherical, tmp_path):
+    # just checking that it runs here
+    image = osmesa_fake_spherical.run()
+
+    outdir = tmp_path / "snapshots"
+    outdir.mkdir()
+    fn = str(outdir / "testout.png")
+    _ = yt.write_bitmap(image, fn)
+    assert os.path.exists(fn)

From 4f73cdd61dfab86347f69188a36e9f1f1b08de53 Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Thu, 8 Jun 2023 06:34:50 -0500
Subject: [PATCH 23/35] only rescale for cartesian

---
 yt_idv/scene_data/block_collection.py | 30 +++++++++++++++++----------
 1 file changed, 19 insertions(+), 11 deletions(-)

diff --git a/yt_idv/scene_data/block_collection.py b/yt_idv/scene_data/block_collection.py
index 1d5b6a87..206b0918 100644
--- a/yt_idv/scene_data/block_collection.py
+++ b/yt_idv/scene_data/block_collection.py
@@ -82,11 +82,15 @@ def add_data(self, field, no_ghost=False):
 
         # Now we set up our buffer
         vert = np.array(vert, dtype="f4")
-        units = self.data_source.ds.units
-        ratio = (units.code_length / units.unitary).base_value
-        dx = np.array(dx, dtype="f4") * ratio
-        le = np.array(le, dtype="f4") * ratio
-        re = np.array(re, dtype="f4") * ratio
+        dx = np.array(dx, dtype="f4")
+        le = np.array(le, dtype="f4")
+        re = np.array(re, dtype="f4")
+        if self._yt_geom_str == "cartesian":
+            units = self.data_source.ds.units
+            ratio = (units.code_length / units.unitary).base_value
+            dx = dx * ratio
+            le = le * ratio
+            re = re * ratio
         self._set_geometry_attributes(le, re)
         self.vertex_array.attributes.append(
             VertexAttribute(name="model_vertex", data=vert)
@@ -102,15 +106,19 @@ def add_data(self, field, no_ghost=False):
         # Now we set up our textures
         self._load_textures()
 
+    @property
+    def _yt_geom_str(self):
+        # note: casting to string for compatibility with new and old geometry
+        # attributes (now an enum member in latest yt),
+        # see https://github.com/yt-project/yt/pull/4244
+        return str(self.data_source.ds.geometry)
+
     def _set_geometry_attributes(self, le, re):
         # set any vertex_array attributes that depend on the yt geometry type
 
-        # note: these comparisons should change to literal comparisons to enum
-        # members with future yt, see https://github.com/yt-project/yt/pull/4244
-        yt_geom = self.data_source.ds.geometry
-        if yt_geom == "cartesian":
+        if self._yt_geom_str == "cartesian":
             return
-        elif yt_geom == "spherical":
+        elif self._yt_geom_str == "spherical":
             axis_id = self.data_source.ds.coordinates.axis_id
             phi_plane_le = phi_normal_planes(le, axis_id, cast_type="f4")
             phi_plane_re = phi_normal_planes(re, axis_id, cast_type="f4")
@@ -122,7 +130,7 @@ def _set_geometry_attributes(self, le, re):
             )
         else:
             raise NotImplementedError(
-                f"{self.name} does not implement {yt_geom} geometries."
+                f"{self.name} does not implement {self._yt_geom_str} geometries."
             )
 
     def viewpoint_iter(self, camera):

From ac9cd0086507f8430dbdde17f2c1cbd16625711d Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Thu, 8 Jun 2023 06:40:18 -0500
Subject: [PATCH 24/35] move shader constants to include

---
 yt_idv/shaders/header.inc.glsl       | 3 +++
 yt_idv/shaders/ray_tracing.frag.glsl | 2 --
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/yt_idv/shaders/header.inc.glsl b/yt_idv/shaders/header.inc.glsl
index 5ae7f438..5a3e146a 100644
--- a/yt_idv/shaders/header.inc.glsl
+++ b/yt_idv/shaders/header.inc.glsl
@@ -1 +1,4 @@
 #version 330 core
+
+const float INFINITY = 1. / 0.;
+const float PI = 3.1415926535897932384626433832795;
diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index d536df76..1cafec10 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -11,8 +11,6 @@ out vec4 output_color;
 flat in vec4 phi_plane_le;
 flat in vec4 phi_plane_re;
 
-const float INFINITY = 1. / 0.;
-const float PI = 3.1415926535897932384626433832795;
 bool within_bb(vec3 pos)
 {
     bvec3 left =  greaterThanEqual(pos, left_edge);

From 30a8cf466ac0a24ec94949babdaee91c09ab4e04 Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Wed, 14 Jun 2023 22:11:56 -0500
Subject: [PATCH 25/35] temp in progress, broken

---
 yt_idv/shaders/ray_tracing.frag.glsl | 23 ++++++++++++-----------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 1cafec10..6187701a 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -48,7 +48,7 @@ vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
 //   void (vec3 tex_curr_pos, inout vec4 curr_color, float tdelta, float t,
 //         vec3 direction);
 
-float get_ray_plane_intersection(vec3 p_normal, float p_constant, vec3 ray_origin, vec3 ray_dir)
+float get_ray_plane_intersection(inout vec2 t_points, vec3 p_normal, float p_constant, vec3 ray_origin, vec3 ray_dir)
 {
     float n_dot_u = dot(p_normal, ray_dir);
     float n_dot_ro = dot(p_normal, ray_origin);
@@ -64,8 +64,9 @@ float get_ray_plane_intersection(vec3 p_normal, float p_constant, vec3 ray_origi
     return (p_constant - n_dot_ro) / n_dot_u;
 }
 
-vec2 get_ray_sphere_intersection(float r, vec3 ray_origin, vec3 ray_dir)
+void get_ray_sphere_intersection(inout vec2 t_points, float r, vec3 ray_origin, vec3 ray_dir)
 {
+    // gets called first
     float dir_dot_dir = dot(ray_dir, ray_dir);
     float ro_dot_ro = dot(ray_origin, ray_origin);
     float dir_dot_ro = dot(ray_dir, ray_origin);
@@ -75,18 +76,16 @@ vec2 get_ray_sphere_intersection(float r, vec3 ray_origin, vec3 ray_dir)
     float b = 2.0 * dir_dot_ro;
     float c =  ro_dot_ro - rsq;
     float determinate = b*b - 2.0 * a_2 * c;
-    float cutoff_val = 0.0;
-    if (determinate < cutoff_val)
+    float t_point;
+    if (determinate == 0.0)
     {
-        return vec2(INFINITY, INFINITY);
+        t_point = -b / a_2;
+        if within_bbox(ge)
     }
-    else if (determinate == cutoff_val)
+    else if (determinate > 0.0)
     {
-        return vec2(-b / a_2, INFINITY);
-    }
-    else
-    {
-        return vec2((-b - sqrt(determinate))/ a_2, (-b + sqrt(determinate))/ a_2);
+        t_points[0] = (-b - sqrt(determinate))/ a_2;
+        t_points[1] = (-b + sqrt(determinate))/ a_2;
     }
 
 }
@@ -149,6 +148,8 @@ void spherical_coord_shortcircuit()
     vec4 curr_color = vec4(0.0);
 
     // intersections
+    vec2 t_points = vec2(INFINITY, INFINITY);
+
     vec2 t_sphere_outer = get_ray_sphere_intersection(right_edge[id_r], ray_position_xyz, dir);
     if (isinf(t_sphere_outer[0]) && isinf(t_sphere_outer[1]))
     {

From 8006d1babb7eb590caa7e37d7c4aebe1a6c2400e Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Wed, 14 Jun 2023 22:40:28 -0500
Subject: [PATCH 26/35] this... this works?

---
 yt_idv/shaders/ray_tracing.frag.glsl | 132 +++++++++++++++++++++------
 1 file changed, 102 insertions(+), 30 deletions(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 6187701a..6ee1d42d 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -13,8 +13,8 @@ flat in vec4 phi_plane_re;
 
 bool within_bb(vec3 pos)
 {
-    bvec3 left =  greaterThanEqual(pos, left_edge);
-    bvec3 right = lessThanEqual(pos, right_edge);
+    bvec3 left =  greaterThanEqual(pos, left_edge-0.001);
+    bvec3 right = lessThanEqual(pos, right_edge+0.001);
     return all(left) && all(right);
 }
 
@@ -33,6 +33,26 @@ vec3 transform_vec3(vec3 v) {
     }
 }
 
+vec3 reverse_transform_vec3(vec3 v) {
+    if (is_spherical) {
+        vec3 vout = vec3(0.);
+        vout[id_r] = length(v);
+        vout[id_phi] = atan(v[1], v[0]);
+        vout[id_theta] = atan(sqrt(v[1]*v[1] + v[0]*v[0]), v[2]);
+        return vout;
+    } else {
+        return v;
+    }
+}
+
+vec3 ray_position_to_native_coords(float t, vec3 ray_origin, vec3 ray_dir) {
+    vec3 xyz = ray_origin + t * ray_dir;
+    if (is_spherical){
+        return reverse_transform_vec3(xyz);
+    }
+    return xyz;
+}
+
 vec3 get_offset_texture_position(sampler3D tex, vec3 tex_curr_pos)
 {
     ivec3 texsize = textureSize(tex, 0); // lod (mipmap level) always 0?
@@ -48,23 +68,27 @@ vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
 //   void (vec3 tex_curr_pos, inout vec4 curr_color, float tdelta, float t,
 //         vec3 direction);
 
-float get_ray_plane_intersection(inout vec2 t_points, vec3 p_normal, float p_constant, vec3 ray_origin, vec3 ray_dir)
+void get_ray_plane_intersection(inout vec2 t_points, inout int n_points, vec3 p_normal, float p_constant, vec3 ray_origin, vec3 ray_dir)
 {
     float n_dot_u = dot(p_normal, ray_dir);
     float n_dot_ro = dot(p_normal, ray_origin);
     // add check for n_dot_u == 0 (ray is parallel to plane)
-    if (n_dot_u == float(0.0))
+    if (n_dot_u != float(0.0))
     {
-        // the ray is parallel to the plane. there are either no intersections
-        // or infinite intersections. In the second case, we are guaranteed
-        // to intersect one of the other faces, so we can drop this plane.
-        return INFINITY;
+        float t_point = (p_constant - n_dot_ro) / n_dot_u;
+        if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+        {
+            t_points[n_points] = t_point;
+            n_points = n_points + 1;
+        }
     }
+    // otherwise, the ray is parallel to the plane. there are either no intersections
+    // or infinite intersections. In the second case, we are guaranteed
+    // to intersect the other faces, so we can drop this plane.
 
-    return (p_constant - n_dot_ro) / n_dot_u;
 }
 
-void get_ray_sphere_intersection(inout vec2 t_points, float r, vec3 ray_origin, vec3 ray_dir)
+void get_ray_sphere_intersection(inout vec2 t_points, inout int n_points, float r, vec3 ray_origin, vec3 ray_dir)
 {
     // gets called first
     float dir_dot_dir = dot(ray_dir, ray_dir);
@@ -80,17 +104,33 @@ void get_ray_sphere_intersection(inout vec2 t_points, float r, vec3 ray_origin,
     if (determinate == 0.0)
     {
         t_point = -b / a_2;
-        if within_bbox(ge)
+        if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+        {
+            t_points[n_points] = t_point;
+            n_points = n_points + 1;
+        }
     }
     else if (determinate > 0.0)
     {
-        t_points[0] = (-b - sqrt(determinate))/ a_2;
-        t_points[1] = (-b + sqrt(determinate))/ a_2;
+
+        t_point = (-b - sqrt(determinate))/ a_2;
+        if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+        {
+            t_points[n_points] = t_point;
+            n_points = n_points + 1;
+        }
+
+        t_point = (-b + sqrt(determinate))/ a_2;
+        if (n_points < 2 && within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+        {
+            t_points[n_points] = t_point;
+            n_points = n_points + 1;
+        }
     }
 
 }
 
-vec2 get_ray_cone_intersection(float theta, vec3 ray_origin, vec3 ray_dir)
+void get_ray_cone_intersection(inout vec2 t_points, inout int n_points, float theta, vec3 ray_origin, vec3 ray_dir)
 {
     // note: cos(theta) and vhat could be calculated in vertex shader
     float costheta;
@@ -121,13 +161,19 @@ vec2 get_ray_cone_intersection(float theta, vec3 ray_origin, vec3 ray_dir)
     float b = 2.0 * (ro_dot_vhat * dir_dot_vhat - ro_dot_dir*cos2t);
     float c = pow(ro_dot_vhat, 2) - ro_dot_ro*cos2t;
     float determinate = b*b - 2.0 * a_2 * c;
+    float t_point;
     if (determinate < 0.0)
     {
-        return vec2(INFINITY, INFINITY);
+        return;
     }
     else if (determinate == 0.0)
     {
-        return vec2(-b / a_2, INFINITY);
+        t_point = -b / a_2;
+        if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+        {
+            t_points[n_points] = t_point;
+            n_points = n_points + 1;
+        }
     }
     else
     {
@@ -135,7 +181,19 @@ vec2 get_ray_cone_intersection(float theta, vec3 ray_origin, vec3 ray_dir)
         // and not the actual cone. those values should be discarded. But they will
         // fail subsequent bounds checks for interesecting the volume, so we can
         // just handle it there instead of adding another check here.
-        return vec2((-b - sqrt(determinate))/ a_2, (-b + sqrt(determinate))/ a_2);
+        t_point = (-b - sqrt(determinate))/ a_2;
+        if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+        {
+            t_points[n_points] = t_point;
+            n_points = n_points + 1;
+        }
+
+        t_point = (-b + sqrt(determinate))/ a_2;
+        if (n_points < 2 && within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+        {
+            t_points[n_points] = t_point;
+            n_points = n_points + 1;
+        }
     }
 }
 
@@ -149,27 +207,41 @@ void spherical_coord_shortcircuit()
 
     // intersections
     vec2 t_points = vec2(INFINITY, INFINITY);
+    int n_points = 0; // number of intersections found
 
-    vec2 t_sphere_outer = get_ray_sphere_intersection(right_edge[id_r], ray_position_xyz, dir);
-    if (isinf(t_sphere_outer[0]) && isinf(t_sphere_outer[1]))
-    {
-        // if there are no intersections with the outer sphere, then there
-        // will be no intersections with the remaining faces and we can stop
-        // looking.
-        discard;
+    get_ray_sphere_intersection(t_points, n_points, right_edge[id_r], ray_position_xyz, dir);
+    if (n_points < 2){
+        get_ray_sphere_intersection(t_points, n_points, left_edge[id_r], ray_position_xyz, dir);
+    }
+    if (n_points < 2){
+        get_ray_plane_intersection(t_points, n_points, vec3(phi_plane_le), phi_plane_le[3], ray_position_xyz, dir);
+    }
+    if (n_points < 2){
+        get_ray_plane_intersection(t_points, n_points, vec3(phi_plane_re), phi_plane_re[3], ray_position_xyz, dir);
+    }
+    if (n_points < 2){
+        get_ray_cone_intersection(t_points, n_points, right_edge[id_theta], ray_position_xyz, dir);
+    }
+    if (n_points < 2){
+        get_ray_cone_intersection(t_points, n_points, left_edge[id_theta], ray_position_xyz, dir);
     }
 
-    vec2 t_sphere_inner = get_ray_sphere_intersection(left_edge[id_r], ray_position_xyz, dir);
-    float t_p_1 = get_ray_plane_intersection(vec3(phi_plane_le), phi_plane_le[3], ray_position_xyz, dir);
-    float t_p_2 = get_ray_plane_intersection(vec3(phi_plane_re), phi_plane_re[3], ray_position_xyz, dir);
-    vec2 t_cone_outer = get_ray_cone_intersection(right_edge[id_theta], ray_position_xyz, dir);
-    vec2 t_cone_inner= get_ray_cone_intersection(left_edge[id_theta], ray_position_xyz, dir);
-
+    if (n_points < 1) {
+        discard;
+    }
     // do the texture evaluation in the native coordinate space
     vec3 range = (right_edge + dx/2.0) - (left_edge - dx/2.0);
     vec3 nzones = range / dx;
     vec3 ndx = 1.0/nzones;
 
+    // take those t values, get the spherical position for sampling
+    float t_mid;
+    if (n_points == 1){
+        t_mid = t_points[0];
+    } else {
+        t_mid = (t_points[1] + t_points[0])/2.;
+    }
+    ray_position = ray_position_to_native_coords(t_mid, ray_position, dir);
     vec3 tex_curr_pos = (ray_position - left_edge) / range;
 
     tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;

From c35dd52a57e78b60c940cb4ddc502f81dfa23d0b Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Thu, 15 Jun 2023 09:28:09 -0500
Subject: [PATCH 27/35] re-order the intersection checks

---
 yt_idv/shaders/ray_tracing.frag.glsl | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 6ee1d42d..c3389ff2 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -209,19 +209,29 @@ void spherical_coord_shortcircuit()
     vec2 t_points = vec2(INFINITY, INFINITY);
     int n_points = 0; // number of intersections found
 
+    // outer sphere
     get_ray_sphere_intersection(t_points, n_points, right_edge[id_r], ray_position_xyz, dir);
-    if (n_points < 2){
-        get_ray_sphere_intersection(t_points, n_points, left_edge[id_r], ray_position_xyz, dir);
-    }
+    // note: the order of geometry intersections is important. we are not
+    // explicitly handling the case of a 3 or 4-point intersection where a ray
+    // intersects the two phi-planes and the inner sphere. But by checking the
+    // phi-planes first, we will at least sample some of those points.
+    // left edge fixed phi-plane
     if (n_points < 2){
         get_ray_plane_intersection(t_points, n_points, vec3(phi_plane_le), phi_plane_le[3], ray_position_xyz, dir);
     }
+    // right edge fixed phi-plane
     if (n_points < 2){
         get_ray_plane_intersection(t_points, n_points, vec3(phi_plane_re), phi_plane_re[3], ray_position_xyz, dir);
     }
+    // inner sphere
+    if (n_points < 2){
+        get_ray_sphere_intersection(t_points, n_points, left_edge[id_r], ray_position_xyz, dir);
+    }
+    // outer cone surface
     if (n_points < 2){
         get_ray_cone_intersection(t_points, n_points, right_edge[id_theta], ray_position_xyz, dir);
     }
+    // inner cone surface
     if (n_points < 2){
         get_ray_cone_intersection(t_points, n_points, left_edge[id_theta], ray_position_xyz, dir);
     }

From 90f7f650eacea73665b0c2eb9bf2272825762781 Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Thu, 15 Jun 2023 09:51:01 -0500
Subject: [PATCH 28/35] ray tracing... kind of works? not super stable

---
 yt_idv/shaders/ray_tracing.frag.glsl | 75 ++++++++++++++++++++++++----
 1 file changed, 65 insertions(+), 10 deletions(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index c3389ff2..38d201bb 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -246,21 +246,76 @@ void spherical_coord_shortcircuit()
 
     // take those t values, get the spherical position for sampling
     float t_mid;
+
     if (n_points == 1){
+
+        // single intersection: on cusp. just return a single sample.
         t_mid = t_points[0];
-    } else {
-        t_mid = (t_points[1] + t_points[0])/2.;
+
+        ray_position = ray_position_to_native_coords(t_mid, ray_position, dir);
+        vec3 tex_curr_pos = (ray_position - left_edge) / range;
+
+        tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
+        bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
+        if (sampled) {
+            output_color = curr_color;
+        } else {
+            output_color = vec4(0.0);
+        }
+        return;
     }
-    ray_position = ray_position_to_native_coords(t_mid, ray_position, dir);
-    vec3 tex_curr_pos = (ray_position - left_edge) / range;
 
-    tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
-    bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
-    if (sampled) {
-        output_color = curr_color;
-    } else {
-        output_color = vec4(0.0);
+    // sample once at the midpoint (works?)
+//    t_mid = (t_points[1] + t_points[0])/2.;
+//    ray_position = ray_position_to_native_coords(t_mid, ray_position, dir);
+//    vec3 tex_curr_pos = (ray_position - left_edge) / range;
+//
+//    tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
+//    bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
+//    if (sampled) {
+//        output_color = curr_color;
+//    } else {
+//        output_color = vec4(0.0);
+//    }
+//    return;
+
+    // ray tracing here
+    vec3 ray_origin = ray_position; // should this be p0, camera instead?
+    float t_start = min(t_points[0], t_points[1]);
+    float t_end = max(t_points[0], t_points[1]);
+    float n_samples = 4.0;
+    float dt = (t_end - t_start)/n_samples;
+    float current_t = t_start;
+    bool ever_sampled = false;
+    bool sampled;
+    vec4 v_clip_coord;
+    float f_ndc_depth;
+    float depth = 1.0;
+    while (current_t <= t_end ){
+        ray_position = ray_position_to_native_coords(current_t, ray_origin, dir);
+        vec3 tex_curr_pos = (ray_position - left_edge) / range;
+
+
+        current_t = current_t + dt;
+        tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
+        sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
+//        if (sampled) {
+//            ever_sampled = true;
+//            v_clip_coord = projection * modelview * vec4(transform_vec3(ray_position), 1.0);
+//            f_ndc_depth = v_clip_coord.z / v_clip_coord.w;
+//            depth = min(depth, (1.0 - 0.0) * 0.5 * f_ndc_depth + (1.0 + 0.0) * 0.5);
+//
+//        }
+
     }
+
+    output_color = cleanup_phase(curr_color, dir, t_start, t_end);
+
+//    if (ever_sampled) {
+//        gl_FragDepth = depth;
+//    }
+
+
 }
 
 void main()

From 95c73ebbc0d5a933492073d768d73518dd29cac4 Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Thu, 15 Jun 2023 10:18:00 -0500
Subject: [PATCH 29/35] switch to p0/camera as ray origin

---
 yt_idv/shaders/ray_tracing.frag.glsl | 103 +++++++++++++--------------
 1 file changed, 50 insertions(+), 53 deletions(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 38d201bb..8a8b2182 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -206,34 +206,34 @@ void spherical_coord_shortcircuit()
     vec4 curr_color = vec4(0.0);
 
     // intersections
-    vec2 t_points = vec2(INFINITY, INFINITY);
+    vec2 t_points;
     int n_points = 0; // number of intersections found
 
     // outer sphere
-    get_ray_sphere_intersection(t_points, n_points, right_edge[id_r], ray_position_xyz, dir);
+    get_ray_sphere_intersection(t_points, n_points, right_edge[id_r], p0, dir);
     // note: the order of geometry intersections is important. we are not
     // explicitly handling the case of a 3 or 4-point intersection where a ray
     // intersects the two phi-planes and the inner sphere. But by checking the
     // phi-planes first, we will at least sample some of those points.
     // left edge fixed phi-plane
     if (n_points < 2){
-        get_ray_plane_intersection(t_points, n_points, vec3(phi_plane_le), phi_plane_le[3], ray_position_xyz, dir);
+        get_ray_plane_intersection(t_points, n_points, vec3(phi_plane_le), phi_plane_le[3], p0, dir);
     }
     // right edge fixed phi-plane
     if (n_points < 2){
-        get_ray_plane_intersection(t_points, n_points, vec3(phi_plane_re), phi_plane_re[3], ray_position_xyz, dir);
+        get_ray_plane_intersection(t_points, n_points, vec3(phi_plane_re), phi_plane_re[3], p0, dir);
     }
     // inner sphere
     if (n_points < 2){
-        get_ray_sphere_intersection(t_points, n_points, left_edge[id_r], ray_position_xyz, dir);
+        get_ray_sphere_intersection(t_points, n_points, left_edge[id_r], p0, dir);
     }
     // outer cone surface
     if (n_points < 2){
-        get_ray_cone_intersection(t_points, n_points, right_edge[id_theta], ray_position_xyz, dir);
+        get_ray_cone_intersection(t_points, n_points, right_edge[id_theta], p0, dir);
     }
     // inner cone surface
     if (n_points < 2){
-        get_ray_cone_intersection(t_points, n_points, left_edge[id_theta], ray_position_xyz, dir);
+        get_ray_cone_intersection(t_points, n_points, left_edge[id_theta], p0, dir);
     }
 
     if (n_points < 1) {
@@ -245,14 +245,10 @@ void spherical_coord_shortcircuit()
     vec3 ndx = 1.0/nzones;
 
     // take those t values, get the spherical position for sampling
-    float t_mid;
-
     if (n_points == 1){
 
         // single intersection: on cusp. just return a single sample.
-        t_mid = t_points[0];
-
-        ray_position = ray_position_to_native_coords(t_mid, ray_position, dir);
+        ray_position = ray_position_to_native_coords(t_points[0], p0, dir);
         vec3 tex_curr_pos = (ray_position - left_edge) / range;
 
         tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
@@ -266,50 +262,51 @@ void spherical_coord_shortcircuit()
     }
 
     // sample once at the midpoint (works?)
-//    t_mid = (t_points[1] + t_points[0])/2.;
-//    ray_position = ray_position_to_native_coords(t_mid, ray_position, dir);
-//    vec3 tex_curr_pos = (ray_position - left_edge) / range;
+    float t_mid = (t_points[1] + t_points[0])/2.;
+    ray_position = ray_position_to_native_coords(t_mid, p0, dir);
+    vec3 tex_curr_pos = (ray_position - left_edge) / range;
+
+    tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
+    bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
+    if (sampled) {
+        output_color = curr_color;
+    } else {
+        output_color = vec4(0.0);
+    }
+    return;
+
+    // ray tracing here... hmm... not super stable
+//    vec3 ray_origin = ray_position; // should this be p0, camera instead?
+//    float t_start = min(t_points[0], t_points[1]);
+//    float t_end = max(t_points[0], t_points[1]);
+//    float n_samples = 4.0;
+//    float dt = (t_end - t_start)/n_samples;
+//    float current_t = t_start;
+////    bool ever_sampled = false;
+//    bool sampled;
+//    vec3 tex_curr_pos;
+////    vec4 v_clip_coord;
+////    float f_ndc_depth;
+////    float depth = 1.0;
+//    while (current_t <= t_end ){
+//        ray_position = ray_position_to_native_coords(current_t, p0, dir);
+//        tex_curr_pos = (ray_position - left_edge) / range;
+////        if (tex_curr_pos[0] >0 && tex_curr_pos[1] > 0 && tex_curr_pos[2] >0){
+//        tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
+//        sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
 //
-//    tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
-//    bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
-//    if (sampled) {
-//        output_color = curr_color;
-//    } else {
-//        output_color = vec4(0.0);
+//            //        if (sampled) {
+//            //            ever_sampled = true;
+//            //            v_clip_coord = projection * modelview * vec4(transform_vec3(ray_position), 1.0);
+//            //            f_ndc_depth = v_clip_coord.z / v_clip_coord.w;
+//            //            depth = min(depth, (1.0 - 0.0) * 0.5 * f_ndc_depth + (1.0 + 0.0) * 0.5);
+//            //
+//            //        }
+////        }
+//        current_t = current_t + dt;
 //    }
-//    return;
-
-    // ray tracing here
-    vec3 ray_origin = ray_position; // should this be p0, camera instead?
-    float t_start = min(t_points[0], t_points[1]);
-    float t_end = max(t_points[0], t_points[1]);
-    float n_samples = 4.0;
-    float dt = (t_end - t_start)/n_samples;
-    float current_t = t_start;
-    bool ever_sampled = false;
-    bool sampled;
-    vec4 v_clip_coord;
-    float f_ndc_depth;
-    float depth = 1.0;
-    while (current_t <= t_end ){
-        ray_position = ray_position_to_native_coords(current_t, ray_origin, dir);
-        vec3 tex_curr_pos = (ray_position - left_edge) / range;
-
-
-        current_t = current_t + dt;
-        tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
-        sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
-//        if (sampled) {
-//            ever_sampled = true;
-//            v_clip_coord = projection * modelview * vec4(transform_vec3(ray_position), 1.0);
-//            f_ndc_depth = v_clip_coord.z / v_clip_coord.w;
-//            depth = min(depth, (1.0 - 0.0) * 0.5 * f_ndc_depth + (1.0 + 0.0) * 0.5);
 //
-//        }
-
-    }
-
-    output_color = cleanup_phase(curr_color, dir, t_start, t_end);
+//    output_color = cleanup_phase(curr_color, dir, t_start, t_end);
 
 //    if (ever_sampled) {
 //        gl_FragDepth = depth;

From 1f40671e76956ab65c85de90089a374e231ed40b Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Thu, 15 Jun 2023 10:28:57 -0500
Subject: [PATCH 30/35] set depth buffer

---
 yt_idv/shaders/ray_tracing.frag.glsl | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 8a8b2182..7015055b 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -202,6 +202,7 @@ void spherical_coord_shortcircuit()
     vec3 ray_position = v_model.xyz; // now spherical
     vec3 ray_position_xyz = transform_vec3(ray_position); // cartesian
     vec3 p0 = camera_pos.xyz; // cartesian
+//    vec3 p0 = ray_position_xyz;
     vec3 dir = -normalize(camera_pos.xyz - ray_position_xyz);
     vec4 curr_color = vec4(0.0);
 
@@ -244,6 +245,10 @@ void spherical_coord_shortcircuit()
     vec3 nzones = range / dx;
     vec3 ndx = 1.0/nzones;
 
+    vec4 v_clip_coord;
+    float f_ndc_depth;
+    float depth = 1.0;
+
     // take those t values, get the spherical position for sampling
     if (n_points == 1){
 
@@ -255,6 +260,9 @@ void spherical_coord_shortcircuit()
         bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
         if (sampled) {
             output_color = curr_color;
+            v_clip_coord = projection * modelview * vec4(transform_vec3(ray_position), 1.0);
+            f_ndc_depth = v_clip_coord.z / v_clip_coord.w;
+            gl_FragDepth = min(depth, (1.0 - 0.0) * 0.5 * f_ndc_depth + (1.0 + 0.0) * 0.5);
         } else {
             output_color = vec4(0.0);
         }
@@ -270,6 +278,9 @@ void spherical_coord_shortcircuit()
     bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
     if (sampled) {
         output_color = curr_color;
+        v_clip_coord = projection * modelview * vec4(transform_vec3(ray_position), 1.0);
+        f_ndc_depth = v_clip_coord.z / v_clip_coord.w;
+        gl_FragDepth = min(depth, (1.0 - 0.0) * 0.5 * f_ndc_depth + (1.0 + 0.0) * 0.5);
     } else {
         output_color = vec4(0.0);
     }

From d229ba80f21782fe1dd279aa61f3286492530168 Mon Sep 17 00:00:00 2001
From: Chris Havlin <chris.havlin@gmail.com>
Date: Thu, 15 Jun 2023 15:33:34 -0500
Subject: [PATCH 31/35] more experimenting more chaos

---
 yt_idv/shaders/ray_tracing.frag.glsl | 109 +++++++++++++++++++++------
 1 file changed, 86 insertions(+), 23 deletions(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 7015055b..8774cbc4 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -13,8 +13,8 @@ flat in vec4 phi_plane_re;
 
 bool within_bb(vec3 pos)
 {
-    bvec3 left =  greaterThanEqual(pos, left_edge-0.001);
-    bvec3 right = lessThanEqual(pos, right_edge+0.001);
+    bvec3 left =  greaterThanEqual(pos, left_edge);
+    bvec3 right = lessThanEqual(pos, right_edge);
     return all(left) && all(right);
 }
 
@@ -24,8 +24,8 @@ vec3 transform_vec3(vec3 v) {
             // in yt, phi is the polar angle from (0, 2pi), theta is the azimuthal
             // angle (0, pi). the id_ values below are uniforms that depend on the
             // yt dataset coordinate ordering
-            v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
             v[id_r] * sin(v[id_theta]) * sin(v[id_phi]),
+            v[id_r] * sin(v[id_theta]) * cos(v[id_phi]),
             v[id_r] * cos(v[id_theta])
         );
     } else {
@@ -35,10 +35,40 @@ vec3 transform_vec3(vec3 v) {
 
 vec3 reverse_transform_vec3(vec3 v) {
     if (is_spherical) {
+        // glsl atan docs.
+//        atan returns either the angle whose trigonometric arctangent is
+//        y/x or y_over_x, depending on which overload is invoked.
+//        In the first overload, the signs of y and x are used to determine
+//        the quadrant that the angle lies in. The value returned by atan in
+//        this case is in the range [−π,π]. The result is undefined if x=0
+//
+//.
+//
+//       For the second overload, atan returns the angle whose tangent
+//        is y_over_x. The value returned in this case is in the range
+//       [−π/2,π/2]
+//.
         vec3 vout = vec3(0.);
         vout[id_r] = length(v);
-        vout[id_phi] = atan(v[1], v[0]);
-        vout[id_theta] = atan(sqrt(v[1]*v[1] + v[0]*v[0]), v[2]);
+        // phi: the 0 to 2pi angle
+        float xy = sqrt(v[0]*v[0] + v[1] * v[1]);
+//        vout[id_phi] = acos(v[1]/xy);// + PI;  // want 0 to 2pi
+//        vout[id_phi] = atan(v[1], v[0]) + PI;//
+        // theta: the polar 0 to pi angle
+        vout[id_theta] = acos(v[2]/vout[id_r]);//atan(sqrt(v[1]*v[1] + v[0]*v[0]) , v[2]);// + PI/2;
+
+        float phi = acos(abs(v[0]/xy));
+        if (v[0] < 0. && v[1] >0){
+            phi += PI/2;
+        }
+        if (v[0] < 0. && v[1] < 0){
+            phi += PI;
+        }
+        if (v[0] > 0 && v[1] <0){
+            phi = 2 * PI - phi;
+        }
+        vout[id_phi] = phi;
+
         return vout;
     } else {
         return v;
@@ -73,15 +103,18 @@ void get_ray_plane_intersection(inout vec2 t_points, inout int n_points, vec3 p_
     float n_dot_u = dot(p_normal, ray_dir);
     float n_dot_ro = dot(p_normal, ray_origin);
     // add check for n_dot_u == 0 (ray is parallel to plane)
-    if (n_dot_u != float(0.0))
+    if (n_dot_u == float(0.0))
     {
-        float t_point = (p_constant - n_dot_ro) / n_dot_u;
-        if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
-        {
-            t_points[n_points] = t_point;
-            n_points = n_points + 1;
-        }
+        return;
+    }
+
+    float t_point = (p_constant - n_dot_ro) / n_dot_u;
+    if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+    {
+        t_points[n_points] = t_point;
+        n_points = n_points + 1;
     }
+
     // otherwise, the ray is parallel to the plane. there are either no intersections
     // or infinite intersections. In the second case, we are guaranteed
     // to intersect the other faces, so we can drop this plane.
@@ -101,7 +134,10 @@ void get_ray_sphere_intersection(inout vec2 t_points, inout int n_points, float
     float c =  ro_dot_ro - rsq;
     float determinate = b*b - 2.0 * a_2 * c;
     float t_point;
-    if (determinate == 0.0)
+    if (determinate < 0.0) {
+        return;
+    }
+    else if (determinate == 0.0)
     {
         t_point = -b / a_2;
         if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
@@ -110,11 +146,11 @@ void get_ray_sphere_intersection(inout vec2 t_points, inout int n_points, float
             n_points = n_points + 1;
         }
     }
-    else if (determinate > 0.0)
+    else
     {
 
         t_point = (-b - sqrt(determinate))/ a_2;
-        if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+        if (n_points<2 && within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
         {
             t_points[n_points] = t_point;
             n_points = n_points + 1;
@@ -178,11 +214,10 @@ void get_ray_cone_intersection(inout vec2 t_points, inout int n_points, float th
     else
     {
         // note: it is possible to have real solutions that intersect the shadow cone
-        // and not the actual cone. those values should be discarded. But they will
-        // fail subsequent bounds checks for interesecting the volume, so we can
-        // just handle it there instead of adding another check here.
+        // and not the actual cone. those values will fail the bounds checks for
+        // interesecting the volume so no special treatment is needed.
         t_point = (-b - sqrt(determinate))/ a_2;
-        if (within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
+        if (n_points < 2 && within_bb(ray_position_to_native_coords(t_point, ray_origin, ray_dir)))
         {
             t_points[n_points] = t_point;
             n_points = n_points + 1;
@@ -201,8 +236,8 @@ void spherical_coord_shortcircuit()
 {
     vec3 ray_position = v_model.xyz; // now spherical
     vec3 ray_position_xyz = transform_vec3(ray_position); // cartesian
-    vec3 p0 = camera_pos.xyz; // cartesian
-//    vec3 p0 = ray_position_xyz;
+//    vec3 p0 = camera_pos.xyz; // cartesian
+    vec3 p0 = ray_position_xyz;
     vec3 dir = -normalize(camera_pos.xyz - ray_position_xyz);
     vec4 curr_color = vec4(0.0);
 
@@ -237,9 +272,37 @@ void spherical_coord_shortcircuit()
         get_ray_cone_intersection(t_points, n_points, left_edge[id_theta], p0, dir);
     }
 
-    if (n_points < 1) {
-        discard;
+//    if (n_points == 0) {
+//        discard;
+//    }
+    // hmm, losing theta > pi/2 values... NO. loosing phi values
+    if (n_points == 1){
+        output_color = vec4(1.,0., 0., 1.);
+    } else if (n_points == 2){
+        output_color = vec4(0.,1., 0., 1.);
+    } else if (n_points == 0){
+        output_color = vec4(0., 0., 1., 1.);
     }
+//    vec3 roundtrip = reverse_transform_vec3(ray_position_xyz);
+//    output_color = vec4(0., 0., 1., 1.);
+    // correct: id_r, id_theta
+    // incorrect: id_phi
+//    if (abs(roundtrip[id_phi] - ray_position[id_phi]) < 0.01){
+//        output_color = vec4(1., 0, 0., 1.);
+//    }
+//    if (roundtrip[id_phi] > 2 * PI){
+//        output_color = vec4(1., 0., 0., 1);
+//    }
+//    // all y < 0 are being missed.
+//    if (ray_position_xyz[1] > 0) {
+//        output_color = vec4(1., 0., 0., 1);
+//    } else {
+//        output_color = vec4(0., 1., 0., 1);
+//    }
+//    output_color = vec4(ray_position_xyz[1], 0., 0., 1);
+
+    return;
+
     // do the texture evaluation in the native coordinate space
     vec3 range = (right_edge + dx/2.0) - (left_edge - dx/2.0);
     vec3 nzones = range / dx;

From ac5c38246ec540a403e2b68624fc54b21709ba54 Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Fri, 16 Jun 2023 10:11:20 -0500
Subject: [PATCH 32/35] incremental cleanup

---
 yt_idv/shaders/ray_tracing.frag.glsl | 46 ++++++++++++++--------------
 1 file changed, 23 insertions(+), 23 deletions(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 8774cbc4..2ed46896 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -35,28 +35,26 @@ vec3 transform_vec3(vec3 v) {
 
 vec3 reverse_transform_vec3(vec3 v) {
     if (is_spherical) {
-        // glsl atan docs.
-//        atan returns either the angle whose trigonometric arctangent is
-//        y/x or y_over_x, depending on which overload is invoked.
-//        In the first overload, the signs of y and x are used to determine
-//        the quadrant that the angle lies in. The value returned by atan in
-//        this case is in the range [−π,π]. The result is undefined if x=0
-//
-//.
-//
-//       For the second overload, atan returns the angle whose tangent
-//        is y_over_x. The value returned in this case is in the range
-//       [−π/2,π/2]
-//.
+        // glsl atan docs:
+        //        atan returns either the angle whose trigonometric arctangent is
+        //        y/x or y_over_x, depending on which overload is invoked.
+        //        In the first overload, the signs of y and x are used to determine
+        //        the quadrant that the angle lies in. The value returned by atan in
+        //        this case is in the range [−π,π]. The result is undefined if x=0
+        //
+        //       For the second overload, atan returns the angle whose tangent
+        //        is y_over_x. The value returned in this case is in the range
+        //       [−π/2,π/2]
+        //
         vec3 vout = vec3(0.);
         vout[id_r] = length(v);
-        // phi: the 0 to 2pi angle
-        float xy = sqrt(v[0]*v[0] + v[1] * v[1]);
-//        vout[id_phi] = acos(v[1]/xy);// + PI;  // want 0 to 2pi
-//        vout[id_phi] = atan(v[1], v[0]) + PI;//
         // theta: the polar 0 to pi angle
-        vout[id_theta] = acos(v[2]/vout[id_r]);//atan(sqrt(v[1]*v[1] + v[0]*v[0]) , v[2]);// + PI/2;
-
+        vout[id_theta] = acos(v[2]/vout[id_r]);
+        // phi: the 0 to 2pi angle : leaving
+//        vout[id_phi] = atan(v[1], v[0]) + PI;
+        // note: stopped trusting the glsl atan, decided to write my own method
+        // here but the glsl atan is probably fine. need to figure out other things now...
+        float xy = sqrt(v[0]*v[0] + v[1] * v[1]);
         float phi = acos(abs(v[0]/xy));
         if (v[0] < 0. && v[1] >0){
             phi += PI/2;
@@ -272,9 +270,10 @@ void spherical_coord_shortcircuit()
         get_ray_cone_intersection(t_points, n_points, left_edge[id_theta], p0, dir);
     }
 
-//    if (n_points == 0) {
-//        discard;
-//    }
+    if (n_points == 0) {
+        // when discarding, loosing phi values corresponding to y<0...
+        discard;
+    }
     // hmm, losing theta > pi/2 values... NO. loosing phi values
     if (n_points == 1){
         output_color = vec4(1.,0., 0., 1.);
@@ -349,7 +348,8 @@ void spherical_coord_shortcircuit()
     }
     return;
 
-    // ray tracing here... hmm... not super stable
+    // ray tracing here... hmm... not super stable. some lag, occasional crashes.
+    // might be the shader is doing too much or might be something else...
 //    vec3 ray_origin = ray_position; // should this be p0, camera instead?
 //    float t_start = min(t_points[0], t_points[1]);
 //    float t_end = max(t_points[0], t_points[1]);

From 232c21ba11599b8759a86dc56491e3353d647eae Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Fri, 16 Jun 2023 10:12:37 -0500
Subject: [PATCH 33/35] and another increment of cleanup

---
 yt_idv/shaders/ray_tracing.frag.glsl | 27 ---------------------------
 1 file changed, 27 deletions(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 2ed46896..e9868f21 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -274,33 +274,6 @@ void spherical_coord_shortcircuit()
         // when discarding, loosing phi values corresponding to y<0...
         discard;
     }
-    // hmm, losing theta > pi/2 values... NO. loosing phi values
-    if (n_points == 1){
-        output_color = vec4(1.,0., 0., 1.);
-    } else if (n_points == 2){
-        output_color = vec4(0.,1., 0., 1.);
-    } else if (n_points == 0){
-        output_color = vec4(0., 0., 1., 1.);
-    }
-//    vec3 roundtrip = reverse_transform_vec3(ray_position_xyz);
-//    output_color = vec4(0., 0., 1., 1.);
-    // correct: id_r, id_theta
-    // incorrect: id_phi
-//    if (abs(roundtrip[id_phi] - ray_position[id_phi]) < 0.01){
-//        output_color = vec4(1., 0, 0., 1.);
-//    }
-//    if (roundtrip[id_phi] > 2 * PI){
-//        output_color = vec4(1., 0., 0., 1);
-//    }
-//    // all y < 0 are being missed.
-//    if (ray_position_xyz[1] > 0) {
-//        output_color = vec4(1., 0., 0., 1);
-//    } else {
-//        output_color = vec4(0., 1., 0., 1);
-//    }
-//    output_color = vec4(ray_position_xyz[1], 0., 0., 1);
-
-    return;
 
     // do the texture evaluation in the native coordinate space
     vec3 range = (right_edge + dx/2.0) - (left_edge - dx/2.0);

From e1dbdfa6b21920737aca667a063b6a86f8040342 Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Fri, 16 Jun 2023 10:27:55 -0500
Subject: [PATCH 34/35] more cleanup, currently at best config...

---
 examples/amr_spherical_volume_rendering.py |   4 +-
 yt_idv/shaders/ray_tracing.frag.glsl       | 111 ++++++++++-----------
 2 files changed, 57 insertions(+), 58 deletions(-)

diff --git a/examples/amr_spherical_volume_rendering.py b/examples/amr_spherical_volume_rendering.py
index be781183..37b46958 100644
--- a/examples/amr_spherical_volume_rendering.py
+++ b/examples/amr_spherical_volume_rendering.py
@@ -43,8 +43,8 @@
     rc = yt_idv.render_context(height=800, width=800, gui=True)
     dd = ds.all_data()
     dd.max_level = 1
-    sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
+    # sg = rc.add_scene(ds, ("index", "r"), no_ghost=True)
     # sg = rc.add_scene(ds, ("index", "theta"), no_ghost=True)
-    # sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
+    sg = rc.add_scene(ds, ("index", "phi"), no_ghost=True)
     sg.camera.focus = [0.0, 0.0, 0.0]
     rc.run()
diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index e9868f21..4e32b709 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -13,8 +13,8 @@ flat in vec4 phi_plane_re;
 
 bool within_bb(vec3 pos)
 {
-    bvec3 left =  greaterThanEqual(pos, left_edge);
-    bvec3 right = lessThanEqual(pos, right_edge);
+    bvec3 left =  greaterThanEqual(pos, left_edge-.01);
+    bvec3 right = lessThanEqual(pos, right_edge+.01);
     return all(left) && all(right);
 }
 
@@ -51,21 +51,21 @@ vec3 reverse_transform_vec3(vec3 v) {
         // theta: the polar 0 to pi angle
         vout[id_theta] = acos(v[2]/vout[id_r]);
         // phi: the 0 to 2pi angle : leaving
-//        vout[id_phi] = atan(v[1], v[0]) + PI;
+        vout[id_phi] = atan(v[0], v[1]);// + PI;
         // note: stopped trusting the glsl atan, decided to write my own method
         // here but the glsl atan is probably fine. need to figure out other things now...
-        float xy = sqrt(v[0]*v[0] + v[1] * v[1]);
-        float phi = acos(abs(v[0]/xy));
-        if (v[0] < 0. && v[1] >0){
-            phi += PI/2;
-        }
-        if (v[0] < 0. && v[1] < 0){
-            phi += PI;
-        }
-        if (v[0] > 0 && v[1] <0){
-            phi = 2 * PI - phi;
-        }
-        vout[id_phi] = phi;
+//        float xy = sqrt(v[0]*v[0] + v[1] * v[1]);
+//        float phi = acos(abs(v[0]/xy));
+//        if (v[0] < 0. && v[1] >0){
+//            phi += PI/2;
+//        }
+//        if (v[0] < 0. && v[1] < 0){
+//            phi += PI;
+//        }
+//        if (v[0] > 0 && v[1] <0){
+//            phi = 2 * PI - phi;
+//        }
+//        vout[id_phi] = phi;
 
         return vout;
     } else {
@@ -234,8 +234,8 @@ void spherical_coord_shortcircuit()
 {
     vec3 ray_position = v_model.xyz; // now spherical
     vec3 ray_position_xyz = transform_vec3(ray_position); // cartesian
-//    vec3 p0 = camera_pos.xyz; // cartesian
-    vec3 p0 = ray_position_xyz;
+    vec3 p0 = camera_pos.xyz; // cartesian
+//    vec3 p0 = ray_position_xyz;
     vec3 dir = -normalize(camera_pos.xyz - ray_position_xyz);
     vec4 curr_color = vec4(0.0);
 
@@ -283,6 +283,8 @@ void spherical_coord_shortcircuit()
     vec4 v_clip_coord;
     float f_ndc_depth;
     float depth = 1.0;
+    vec3 tex_curr_pos;
+    bool sampled;
 
     // take those t values, get the spherical position for sampling
     if (n_points == 1){
@@ -292,7 +294,7 @@ void spherical_coord_shortcircuit()
         vec3 tex_curr_pos = (ray_position - left_edge) / range;
 
         tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
-        bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
+        sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
         if (sampled) {
             output_color = curr_color;
             v_clip_coord = projection * modelview * vec4(transform_vec3(ray_position), 1.0);
@@ -304,13 +306,13 @@ void spherical_coord_shortcircuit()
         return;
     }
 
-    // sample once at the midpoint (works?)
+    // sample once at the midpoint just to try... works? doesnt crash at least.
+    // but still missing that phi hemisphere.
     float t_mid = (t_points[1] + t_points[0])/2.;
     ray_position = ray_position_to_native_coords(t_mid, p0, dir);
-    vec3 tex_curr_pos = (ray_position - left_edge) / range;
-
+    tex_curr_pos = (ray_position - left_edge) / range;
     tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
-    bool sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
+    sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
     if (sampled) {
         output_color = curr_color;
         v_clip_coord = projection * modelview * vec4(transform_vec3(ray_position), 1.0);
@@ -323,41 +325,38 @@ void spherical_coord_shortcircuit()
 
     // ray tracing here... hmm... not super stable. some lag, occasional crashes.
     // might be the shader is doing too much or might be something else...
-//    vec3 ray_origin = ray_position; // should this be p0, camera instead?
-//    float t_start = min(t_points[0], t_points[1]);
-//    float t_end = max(t_points[0], t_points[1]);
-//    float n_samples = 4.0;
-//    float dt = (t_end - t_start)/n_samples;
-//    float current_t = t_start;
-////    bool ever_sampled = false;
-//    bool sampled;
+    float t_start = min(t_points[0], t_points[1]);
+    float t_end = max(t_points[0], t_points[1]);
+    float n_samples = 4.0;
+    float dt = (t_end - t_start)/n_samples;
+    float current_t = t_start;
+    bool ever_sampled = false;
+//    bool sampled;  // defined above, but would define here if we could delete above.
 //    vec3 tex_curr_pos;
-////    vec4 v_clip_coord;
-////    float f_ndc_depth;
-////    float depth = 1.0;
-//    while (current_t <= t_end ){
-//        ray_position = ray_position_to_native_coords(current_t, p0, dir);
-//        tex_curr_pos = (ray_position - left_edge) / range;
-////        if (tex_curr_pos[0] >0 && tex_curr_pos[1] > 0 && tex_curr_pos[2] >0){
-//        tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
-//        sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
-//
-//            //        if (sampled) {
-//            //            ever_sampled = true;
-//            //            v_clip_coord = projection * modelview * vec4(transform_vec3(ray_position), 1.0);
-//            //            f_ndc_depth = v_clip_coord.z / v_clip_coord.w;
-//            //            depth = min(depth, (1.0 - 0.0) * 0.5 * f_ndc_depth + (1.0 + 0.0) * 0.5);
-//            //
-//            //        }
-////        }
-//        current_t = current_t + dt;
-//    }
-//
-//    output_color = cleanup_phase(curr_color, dir, t_start, t_end);
-
-//    if (ever_sampled) {
-//        gl_FragDepth = depth;
-//    }
+//    vec4 v_clip_coord;
+//    float f_ndc_depth;
+//    float depth = 1.0;
+    while (current_t <= t_end ){
+        ray_position = ray_position_to_native_coords(current_t, p0, dir);
+        tex_curr_pos = (ray_position - left_edge) / range;
+        tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
+        sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));
+
+        if (sampled) {
+            ever_sampled = true;
+            v_clip_coord = projection * modelview * vec4(transform_vec3(ray_position), 1.0);
+            f_ndc_depth = v_clip_coord.z / v_clip_coord.w;
+            depth = min(depth, (1.0 - 0.0) * 0.5 * f_ndc_depth + (1.0 + 0.0) * 0.5);
+
+        }
+        current_t = current_t + dt;
+    }
+
+    output_color = cleanup_phase(curr_color, dir, t_start, t_end);
+
+    if (ever_sampled) {
+        gl_FragDepth = depth;
+    }
 
 
 }

From d4311d6a8cbd002c2d0cd4d142e1b2212b7f43fb Mon Sep 17 00:00:00 2001
From: chrishavlin <chris.havlin@gmail.com>
Date: Fri, 16 Jun 2023 10:33:39 -0500
Subject: [PATCH 35/35] remove redundant declaration

---
 yt_idv/shaders/ray_tracing.frag.glsl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/yt_idv/shaders/ray_tracing.frag.glsl b/yt_idv/shaders/ray_tracing.frag.glsl
index 4e32b709..6f3c2273 100644
--- a/yt_idv/shaders/ray_tracing.frag.glsl
+++ b/yt_idv/shaders/ray_tracing.frag.glsl
@@ -291,7 +291,7 @@ void spherical_coord_shortcircuit()
 
         // single intersection: on cusp. just return a single sample.
         ray_position = ray_position_to_native_coords(t_points[0], p0, dir);
-        vec3 tex_curr_pos = (ray_position - left_edge) / range;
+        tex_curr_pos = (ray_position - left_edge) / range;
 
         tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
         sampled = sample_texture(tex_curr_pos, curr_color, 0.0, 0.0, vec3(0.0));