it compiles! blank though... probably the cartesian scaling...

yt_idv/coordinate_utilities.pyx

@@ -524,3 +524,39 @@ def cartesian_bboxes_edges(MixedCoordBBox bbox_handler,
     xyz_re = (np.asarray(x_re), np.asarray(y_re), np.asarray(z_re))
 
     return xyz_le, xyz_re
+
+
+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

yt_idv/scene_data/block_collection.py

@@ -124,6 +124,7 @@ def _set_geometry_attributes(self, le, re):
             from yt_idv.coordinate_utilities import (
                 SphericalMixedCoordBBox,
                 cartesian_bboxes_edges,
+                phi_normal_planes,
             )
 
             axis_id = self.data_source.ds.coordinates.axis_id
@@ -171,6 +172,17 @@ def _set_geometry_attributes(self, le, re):
             # does not seem that diagonal is used anywhere, but recalculating to
             # be safe...
             self.diagonal = np.sqrt(((re_cart - le_cart) ** 2).sum())
+
+            # pre-calculate the phi-normal planes
+            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 {self._yt_geom_str} geometries."

yt_idv/shaders/grid_expand.geom.glsl

@@ -8,8 +8,12 @@ flat in vec3 vright_edge[];
 #ifdef SPHERICAL_GEOM
 flat in vec3 vleft_edge_cart[];
 flat in vec3 vright_edge_cart[];
+flat in vec4 vphi_plane_le[];
+flat in vec4 vphi_plane_re[];
 flat out vec3 left_edge_cart;
 flat out vec3 right_edge_cart;
+flat out vec4 phi_plane_le;
+flat out vec4 phi_plane_re;
 #endif
 
 flat out vec3 dx;
@@ -75,6 +79,8 @@ void main() {
         #ifdef SPHERICAL_GEOM
         left_edge_cart = vleft_edge_cart[0];
         right_edge_cart = vright_edge_cart[0];
+        phi_plane_le = vphi_plane_le[0];
+        phi_plane_re = vphi_plane_re[0];
         #endif
 
         dx = vdx[0];

yt_idv/shaders/grid_position.vert.glsl

@@ -20,9 +20,15 @@ flat out vec3 vright_edge;
 // pre-computed cartesian le, re
 in vec3 le_cart;
 in vec3 re_cart;
+// pre-computed phi-normal planes
+// first 3 elements are the normal vector, final is constant
+in vec4 phi_plane_le;
+in vec4 phi_plane_re;
 
 flat out vec3 vleft_edge_cart;
 flat out vec3 vright_edge_cart;
+flat out vec4 vphi_plane_le;
+flat out vec4 vphi_plane_re;
 #endif
 
 
@@ -46,5 +52,7 @@ void main()
     // cartesian bounding boxes
     vleft_edge_cart = vec3(le_cart);
     vright_edge_cart = vec3(re_cart);
+    vphi_plane_le = vec4(phi_plane_le);
+    vphi_plane_re = vec4(phi_plane_re);
     #endif
 }

yt_idv/shaders/ray_tracing.frag.glsl

@@ -10,6 +10,8 @@ flat in ivec3 texture_offset;
 #ifdef SPHERICAL_GEOM
 flat in vec3 left_edge_cart;
 flat in vec3 right_edge_cart;
+flat in vec4 phi_plane_le;
+flat in vec4 phi_plane_re;
 #endif
 
 out vec4 output_color;
@@ -22,20 +24,30 @@ bool within_bb(vec3 pos)
 }
 
 #ifdef SPHERICAL_GEOM
+
+vec3 scale_cart(vec3 v) {
+    vec3 vout;
+    vout[0] = v[0] * cart_bbox_max_width + cart_bbox_le[0];
+    vout[1] = v[1] * cart_bbox_max_width + cart_bbox_le[1];
+    vout[2] = v[2] * cart_bbox_max_width + cart_bbox_le[2];
+    return vout;
+}
+
 vec3 cart_to_sphere_vec3(vec3 v) {
     // transform a single point in cartesian coords to spherical
-    vec3 vout = vec3(0.,0.,0.);
+    vec3 vout;
+    vec3 vsc = scale_cart(v); // re-scaled cartesian
 
     // 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, cart_bbox_* variables are also uniforms
     float x = v[0] * cart_bbox_max_width + cart_bbox_le[0];
     float y = v[1] * cart_bbox_max_width + cart_bbox_le[1];
     float z = v[2] * cart_bbox_max_width + cart_bbox_le[2];
-    vout[id_r] = x * x + y * y + z * z;
+    vout[id_r] = vsc[0] * vsc[0] + vsc[1] * vsc[1] + vsc[2] * vsc[2];
     vout[id_r] = sqrt(vout[id_r]);
-    vout[id_theta] = acos(z / vout[id_r]);
-    float phi = atan(y, x);
+    vout[id_theta] = acos(vsc[2] / vout[id_r]);
+    float phi = atan(vsc[1], vsc[0]);
     // atan2 returns -pi to pi, adjust to (0, 2pi)
     if (phi < 0 ){
         phi = phi + 2.0 * PI;
@@ -84,7 +96,8 @@ vec2 quadratic_eval(float b, float a_2, float c){
     b = -1. * b/a_2;
 
     // do the calculation
-    vec2 return_vec = vec2((b - det)*det_is_nonneg + null_val, (b + det)*det_is_nonneg + null_val);
+    vec2 return_vec;
+    return_vec = vec2((b - det)*det_is_nonneg + null_val, (b + det)*det_is_nonneg + null_val);
     // override the second return if det is 0.
     return_vec[1] = return_vec[1] * (1. - det_is_zero)  - 99. * det_is_zero;
 
@@ -222,37 +235,39 @@ void main()
     t0 = max(t0, 0.0);
     if (t1 <= t0) discard;
 
-    vec4 t_control_points = vec4(-99.0)
+    vec4 t_control_points = vec4(-99.0);
     #ifdef SPHERICAL_GEOM
     // here, we check the intersections with the primitive geometries describing the
     // surfaces of the spherical volume element. The intersections are only saved if
     // they fall within the ray parameter range given by the initial slab test
     // there are 0, 1, 2 or 4 intersections possible. 4 intersections is annoying.
-    vec2 t_temp2
-    int n_extra = 0
-    float t2 = -99.0
-    float t3 = -99.0
+    vec2 t_temp2;
+    int n_extra = 0;
+    float t2 = -99.0;
+    float t3 = -99.0;
+    vec3 ray_position_xyz = scale_cart(ray_position);
 
+    // arg... dir, ray_position_xyz sacling needs to be handled probably?
     t_temp2 = get_ray_sphere_intersection(right_edge[id_r], ray_position_xyz, dir);
-    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1)
-    store_temp_intx(n_extra, t_control_points, t_temp2[1], t0, t1)
+    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1);
+    store_temp_intx(n_extra, t_control_points, t_temp2[1], t0, t1);
 
     t_temp2 = get_ray_sphere_intersection(left_edge[id_r], ray_position_xyz, dir);
-    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1)
-    store_temp_intx(n_extra, t_control_points, t_temp2[1], t0, t1)
+    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1);
+    store_temp_intx(n_extra, t_control_points, t_temp2[1], t0, t1);
 
     t_temp2[0] = get_ray_plane_intersection(vec3(phi_plane_le), phi_plane_le[3], ray_position_xyz, dir);
-    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1)
+    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1);
     t_temp2[0] = get_ray_plane_intersection(vec3(phi_plane_re), phi_plane_re[3], ray_position_xyz, dir);
-    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1)
+    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1);
 
     t_temp2 = get_ray_cone_intersection(right_edge[id_theta], ray_position_xyz, dir);
-    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1)
-    store_temp_intx(n_extra, t_control_points, t_temp2[1], t0, t1)
+    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1);
+    store_temp_intx(n_extra, t_control_points, t_temp2[1], t0, t1);
 
     t_temp2 = get_ray_cone_intersection(left_edge[id_theta], ray_position_xyz, dir);
-    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1)
-    store_temp_intx(n_extra, t_control_points, t_temp2[1], t0, t1)
+    store_temp_intx(n_extra, t_control_points, t_temp2[0], t0, t1);
+    store_temp_intx(n_extra, t_control_points, t_temp2[1], t0, t1);
 
     // at this point, t_control_points will contain 1, 2 or 4 intersections that fall within
     // the bounding cartesian box and are guaranteed to be > 0 if they are set.
@@ -266,11 +281,11 @@ void main()
     float full_min = min_of_vec4(t_control_points);
 
     // store the other 2 points
-    vec2 mid_ts = vec2(-99.0);
+    vec2 mid_ts = vec2(-99.0, -99.0);
     int i_mid = 0;
     for (int i=0;i<4;++i)
     {
-        if (t_control_points[i] > full_min and t_control_points[i] < full_max){
+        if (t_control_points[i] > full_min && t_control_points[i] < full_max){
             mid_ts[i_mid] = t_control_points[i];
             i_mid += 1;
         }
@@ -299,9 +314,9 @@ void main()
 
     // initialize ray tracing loop variables
 
-    vec3 p0  // cartesian position at t = t0
-    vec3 p1  // cartesian position at t = t1
-    float t  // current value of ray parameter t
+    vec3 p0;  // cartesian position at t = t0
+    vec3 p1;  // cartesian position at t = t1
+    float t;  // current value of ray parameter t
 
     bool sampled;
     bool ever_sampled = false;
@@ -317,9 +332,10 @@ void main()
 
     // Some more discussion of this here:
     //  http://prideout.net/blog/?p=64
-    for (int ipart =0; 1; ipart++ipart){
-        t0 = t_control_points[ipart + 2 * ipart]  // index 0 or 2
-        t1 = t_control_points[ipart + 1 + 2 * ipart] // index 1 or 3
+    for (int ipart=0;ipart<2;++ipart)
+    {
+        t0 = t_control_points[ipart + 2 * ipart];  // index 0 or 2
+        t1 = t_control_points[ipart + 1 + 2 * ipart]; // index 1 or 3
 
         t0 = max(t0, 0.0); // if t0 = -99., will have t > t1 and loop wont run below.
         p0 = camera_pos.xyz + dir * t0;