# Implementation of module {rootray_cart}.
# Last edited on 2021-02-18 22:24:28 by jstolfi

import rootray
import hacks
import rn
import sys
from math import sqrt, sin, cos, inf, nan, pi

def halfspace(p, q, k):
  assert p != q, "points must be distinct"
  d = len(p)
  assert 0 <= k and k < d
  # Equation of the
  # Coefficients of equation of the hyperplane {r(t)[k] = A t + B}:
  pk = p[k]
  qk = q[k]
  A =  qk - pk
  B = pk
  if A == 0 or abs(B/A) > 1.0e+100:
    # Line is nearly parallel to the hyperplane:
    if B > 0:
      # Pretend whole line is outside:
      f = rootray.vacuum()
    elif B < 0:
      # Pretend whole line is inside:
      f = rootray.plenum()
    else:
      # Line practically lies on hyperplane. Pretend it enters at {t = 0}:
      f = (+1, [ 0 ])
  else:
    # Compute intersection
    t = -B/A
    # Decide whether {r(-oo)} is inside or outside:
    if A < 0:
      # {r(-oo)} is outside:
      f = (+1, [t])
    else:
      f = (-1, [t])
  return f

def slab(p, q, k):
  assert p != q, "points must be distinct"
  d = len(p)
  assert 0 <= k and k < d
  pk = p[k]
  qk = q[k]
  f1 = halfspace((+pk - 1,), (+qk - 1,), 0)
  f2 = halfspace((-pk - 1,), (-qk - 1,), 0)
  f = rootray.intersection(f1, f2)
  return f

def cube(p, q):
  assert p != q, "points must be distinct"
  d = len(p)
  f = rootray.plenum()
  for k in range(d):
    fk = slab(p, q, k)
    f = rootray.intersection(f, fk)
  return f
  
def ball(p, q):
  assert p != q, "points must be distinct"
  d = len(p)
  # Let {r(t)} be {t*p + (1-t)*q}. 
  # Determine the coefficients of the quadratic formula
  # {A t^2 + B t + C} for the characteristic function
  # {|r(t)|^2 - 1}
  v = rn.sub(q, p)
  A = rn.dot(v,v)
  B = 2*rn.dot(v,p)
  C = rn.dot(p,p) - 1
  assert A > 0
  t0, t1 = hacks.real_quadratic_roots(A, B, C)
  if t0 == None and t1 == None:
    # Line does not cross the circle:
    f = rootray.vacuum()
  else:
    assert t0 != None and t1 != None
    assert t0 < t1
    if t0 == -inf and t1 == +inf:
      # Roots are very far apart:
      f = rootray.plenum()
    elif t0 == -inf:
      # Line entrance is very far in the past:
      f = (-1, [t1])
    elif t1 == +inf:
      # Line exit is very far in the future:
      f = (+1, [t0])
    else: 
      f = (+1, [t0, t1])
  return f
  
def cylinder(p, q, m):
  d = len(p)
  assert 0 <= m and m <= d
  if m == 0:
    f = rootray.plenum()
  elif m == d:
    f = ball(p, q)
  else:
    pm = p[0:m];
    qm = q[0:m];
    if pm == qm:
      # Line {p--q} is parallel to the cylinder:
      Rp = rn.norm_sqr(pm)
      if Rp > 1:
        # Line is outside:
        f = rootray.vacuum()
      elif RP < 1:
        # Line is inside:
        f = rootray.plenum()
      else:
        # Line is on surface: pretend it crosses at {t = 0}:
        f = (+1, [0])
    else:
      f = ball(pm, qm)
  return f 

def cone(p, q, m):
  assert p != q, "points must be distinct"
  d = len(p)
  assert 0 <= m and m <= d
  if m == 0:
    f = rootray.plenum()
  elif m == d:
    f = rootray.vacuum()
  else:
    # Let {r(t)} be {t*p + (1-t)*q}. 
    # Determine the coefficients of the quadratic formula
    # {A t^2 + B t + C} for the characteristic function
    # {|r(t)|^2 - 1}
    v = rn.sub(q, p)
    vm = v[0:m]; vn = v[m:]
    pm = p[0:m]; pn = p[m:]
    A = rn.dot(vm,vm) - rn.dot(vn,vn)
    B = 2*(rn.dot(vm,pm) - rn.dot(vn,pn))
    C = rn.dot(pm,pm) - rn.dot(pn,pn)
    if A == 0 or abs(A) < 1.0e-16:
      # Line is on the surface?
      if B == 0 or abs(B) < 1.0e-16:
        if C > 0:
          f = rootray.vacuum()
        elif C < 0:
          f = rootray.plenum()
        else:
          f = (+1, [0])
      else:
        s = -1 if B > 0 else +1
        f = (s, [ -C/B ])
    else:
      t0, t1 = real_quadratic_roots(A, B, C)
      if t0 == None and t1 == None:
        # Line does not cross the cone:
        f = rootray.vacuum()
      else:
        assert t0 != None and t1 != None
        assert t0 < t1
        s = -1 if A < 0 else +1
        if t0 == -inf and t1 == +inf:
          # Roots are very far apart:
          f = (-s, [])
        elif t0 == -inf:
          # Line entrance is very far in the past:
          f = (-s, [t1])
        elif t1 == +inf:
          # Line exit is very far in the future:
          f = (+s, [t0])
        else: 
          f = (+s, [t0, t1])
    return f