#!/usr/bin/env python
# coding: utf-8

# In[414]:


from IPython.core.display import HTML
import numpy as np
import math


# [x3d](https://pypi.org/project/x3d) is a [PyPI](https://pypi.org) distributed Python package maintained by the 
# [Web3D Consortium](https://web3d.org). See also a local help file [x3d.x3d](x3d.x3d.html) explaining the import statement 
# 
# `from x3d import x3d` .

# In[415]:


from x3d import x3d


# [common_geometry](common_geometry.html) is a Python package, local to this Jupyter notebook folder, which defines convenience classes for generating X3D elements with the geometric definitition pertinent to this notebook. See [common_geometry.py](common_geometry.py) for Python source code; [common_geometry](common_geometry.html) for documentation autogenerated by the [pydoc](https://docs.python.org/3/library/pydoc.html) documentation system.
# 

# In[416]:


import common_geometry


# In[417]:


alpha_degrees = 30.0 # the half angle of the cone
beta_degrees  = 30.0 # the angle the intersecting plane makes with horizontal
a = 1.0              # distance of intersecting plane to cone vertex
f = 2.5              # total height of cone

if (alpha_degrees + beta_degrees >= 90.0):
    raise AssertionError("choice of alpha , beta angles geometrically invalid")




# In[418]:


class InscribedSphere:
    """
    data structure with named fields describing sphere inscribed (internally tangent)
    to a cone with intersection with plane
    """
    def __init__(self):
        self.center = np.array((0.,0.,0.))
        self.radius = 0.0
        # focus is the point of tangency with the plane
        self.focus = np.array((0.,0.,0.))
        # contact_circle is the circular section where the sphere is in tangent contact with the cone
        self.contact_circle_center = ((0.,0.,0.))
        self.contact_circle_radius = 0.
        
    def contact_circle_point(self, theta):
        """
        returns a (3,) array of a point
        on the circle, theta a cclockwise rotation (relative to y axis) from the x axis
        """
        return self.contact_circle_center + self.contact_circle_radius * np.array((math.cos(theta),0.0, math.sin(theta)))


# In[419]:


# convert angles to radians
alpha = math.radians(alpha_degrees)
beta  = math.radians(beta_degrees)


# In[420]:


# coordinates of points and direction vectors
pointC = np.array((0.,0.,0.))
axis   = np.array((0.,-1.,0.))  # axis of cone, pointing to open end
pointB   = pointC - a * axis

z      = np.array((0.,0.,1.0))  # global z axis
normal = np.array((math.sin(beta), math.cos(beta),0.0))

ua     = np.cross(z, normal) # a unit vector in inclined plane, largely in -x direction for beta > 0
                             # will be along semimajor axis of ellipse
ub     = np.cross(normal, ua) # will be along semiminor axis


# ### Evaluation of geometry of the upper inscribed sphere
#     
# ∡CBA = alpha
# 
# ∡ACB = pi/2 - beta
# 
# ∣AD is bisector of ∡BAC
# 
# DE ⊥ AC
# 
# 
# 
# 
# 
# 

# ![Drawing of upper sphere](top_drawing.png)

# In[421]:


angCBA=alpha                       # alpha is half angle of cone

angBAC= math.pi/2 + beta - alpha   # angles in triangle sum to pi
angDAC= angBAC/2                   # AD segment is angle bisector of angle BAC
angACD= math.pi/2 - beta           # beta is angle that segment AC makes with horizontal
angCDA= math.pi - angDAC - angACD  # angles in triangle sum to pi

# now apply law of sines
disBC = a                          # definition of cone dimension a
disAC = math.sin(angCBA) * disBC / math.sin(angBAC) # law of sines
disCD = math.sin(angDAC) * disAC / math.sin(angCDA) # law of sines
disBD = disBC - disCD
disDH = disBD * math.sin(angCBA)
disDF = disDH * math.sin(angCBA)
disFH = disDH * math.cos(angCBA)

disDE = disCD * math.sin(angACD)
disCE = disCD * math.cos(angACD)


# In[ ]:





# In[422]:


upperSphere = InscribedSphere()
upperSphere.radius = disDE
upperSphere.center = pointC - disCD * axis
upperSphere.focus = pointC + disCE * ua
upperSphere.contact_circle_center = upperSphere.center - disDF * axis
upperSphere.contact_circle_radius = disFH


# In[ ]:





# #### Evaluation of geometry of the lower inscribed sphere
# 
# ∡CBA = alpha
# 
# ∡ACB = pi/2 - beta
# 
# ∣AD is bisector of ∡BAC
# 
# DE ⊥ AC
# 

# In[423]:


disBC = a
angCBA=alpha
angACB= math.pi/2 + beta
angBAC= math.pi - angCBA - angACB

angCAD = (math.pi-angBAC)/2
angDCA = math.pi/2 - beta
angADC = math.pi - angCAD - angDCA

disAC = math.sin(angCBA) * disBC/math.sin(angBAC)
disCD = math.sin(angCAD) * disAC/math.sin(angADC)
disAD = math.sin(angDCA) * disAC/math.sin(angADC)
disDE = disAD * math.sin(angCAD)
disCE = disCD * math.cos(angDCA)

disBD = disBC + disCD
disDH = disBD * math.sin(angCBA)
disDF = disDH * math.sin(angCBA)
disFH = disDH * math.cos(angCBA)



# In[424]:


lowerSphere = InscribedSphere()
lowerSphere.radius = disDE
lowerSphere.center = pointC + disCD * axis
lowerSphere.focus = pointC - disCE * ua
lowerSphere.contact_circle_center = lowerSphere.center - disDF * axis
lowerSphere.contact_circle_radius = disFH


# In[425]:


sc = x3d.Scene()
sc.children.append( x3d.Background(skyColor = [(1,1,1)]) )

cone_shape = common_geometry.Cone(
    vertex_height = a,
    total_height  = f,
    half_angle    = alpha,
    appearance = x3d.Appearance(
        material = x3d.Material( diffuseColor=(0.0,0.5,0.5),transparency=0.75)
    )
)


plane_shape = common_geometry.TiltedPlane(
    length = 2.5,
    width  = 2.5, 
    beta   = beta, 
    appearance= x3d.Appearance( material = x3d.Material( emissiveColor=(0.8,0.6,0.6), 
                                                        diffuseColor=(0.8,0.6,0.6),transparency=0.2)) )

upper_sphere = common_geometry.Sphere( 
        center = upperSphere.center,
        radius = upperSphere.radius * 0.995 ,
        appearance = x3d.Appearance( material = x3d.Material( diffuseColor=(1,1,0), transparency=0.0) )
)


upper_focus = common_geometry.Sphere( 
        center = upperSphere.focus,
        radius = 0.02 ,
        appearance = x3d.Appearance( material = x3d.Material( diffuseColor=(0,0,0), transparency=0.0) )
)
contact_circle_appearance = x3d.Appearance( 
        lineProperties = x3d.LineProperties( linewidthScaleFactor=1 ),
        material= x3d.Material( emissiveColor = (0,0,1)) )

upper_contact_circle = common_geometry.Circle(
    center = upperSphere.contact_circle_center,
    radius = upperSphere.contact_circle_radius,
    appearance = contact_circle_appearance
)

lower_sphere = common_geometry.Sphere( 
        center = lowerSphere.center,
        radius = lowerSphere.radius * 0.995,
        appearance = x3d.Appearance( material = x3d.Material( diffuseColor=(1,1,0.0)) )
)

lower_focus = common_geometry.Sphere( 
        center = lowerSphere.focus,
        radius = 0.02 ,
        appearance = x3d.Appearance( material = x3d.Material( diffuseColor=(0,0,0), transparency=0.0) )
)

lower_contact_circle = common_geometry.Circle(
    center = lowerSphere.contact_circle_center,
    radius = lowerSphere.contact_circle_radius,
    appearance = contact_circle_appearance
)


figure = x3d.Transform(
    rotation=(0, -1, 0 ,math.radians(20.0)),
    children = [
        cone_shape,
        plane_shape,
        upper_sphere,
        upper_contact_circle,
        upper_focus,
        lower_sphere,
        lower_contact_circle,
        lower_focus
    ]
)
sc.children.append(figure)


# Explicitly construct the ellipse which is the intersection between the plane and the cone.
# 
# Will start with a circle, then stretch it by different scales in x and z axis to be the appropriate
# ellipse, then rotate it around the z axis

# In[426]:


ellipse_center = 0.5 * (lowerSphere.focus + upperSphere.focus)
# c is 1/2 distance between foci
c = 0.5 * np.sqrt( np.square( lowerSphere.focus - upperSphere.focus).sum())

semimajor = 0.5 * np.sqrt( np.square( lowerSphere.contact_circle_center - upperSphere.contact_circle_center).sum())/ math.cos(alpha)
semiminor = np.sqrt( np.square(semimajor) - np.square(c))

print("semiminor axis: %.5f" % semiminor)
print("semimajor axis: %.5f" % semimajor)

def ellipse_point(theta):
    """
    for theta an angle in radians, returns point on ellipse that makes angle theta measured from center
    relative to major axis
    """
    return ellipse_center + \
            math.cos(theta) * semimajor * ua + \
            math.sin(theta) * semiminor * ub


# In[427]:


ellipse_line_properties = x3d.LineProperties( applied=True, linewidthScaleFactor=3)
ellipse_material = x3d.Material(emissiveColor=(0,0,0))

appearance = x3d.Appearance( lineProperties = ellipse_line_properties, material=ellipse_material )

circular_shape = x3d.Transform(
    rotation = (1.0, 0.0, 0.0, math.pi/2),
    children=[
        x3d.Shape(
            geometry = x3d.Circle2D(radius=1.0),
            appearance = appearance
        )
    ]
)

ellipse_shape = x3d.Transform(
    translation = tuple(ellipse_center),
    children=[x3d.Transform(
        rotation=(0 , 0, -1, beta),
        children=[
            x3d.Transform(
                scale = (semimajor,1,semiminor),
                children = [circular_shape]
            )
        ]
    )             ]
)
figure.children.append(ellipse_shape)


# In[428]:


# create polyline: upper focus to ellipse to upper contact circle
theta = math.pi * 0.75
ep = ellipse_point(theta)
phi = math.atan2(ep[2], ep[0])
print("phi %.3f" % phi)
polyline_vertices = x3d.Coordinate(
[tuple(p) for p in \
[
    upperSphere.focus,
    ellipse_point(theta),
    upperSphere.contact_circle_point(phi),
    lowerSphere.focus,
    ellipse_point(theta),
    lowerSphere.contact_circle_point(phi),    
]
]
)
polyline_vertices.DEF='foci_segments'

polyline_indices = [0,1,2,-1,3,4,5,-1]
polyline_colors  = x3d.Color(color=[(0,0,1),(1,0,0)])

polyline = x3d.Shape(
    geometry = x3d.IndexedLineSet(
        colorPerVertex=False,
        coord = polyline_vertices,
        color = polyline_colors,
        coordIndex = polyline_indices,
        colorIndex = [0,1],
    ),
    appearance = x3d.Appearance( lineProperties = ellipse_line_properties, material=ellipse_material )
)

figure.children.append(polyline)


# In[429]:


NumKeyFrame = 32

keys = np.linspace(0.0, 1.0, NumKeyFrame+1, endpoint=True)
key_values = list()
for frac in keys:
    theta = 2. * math.pi * frac
    ep = ellipse_point(theta)
    phi = math.atan2(ep[2], ep[0])
    value = [tuple(p) for p in \
    [
    upperSphere.focus,
    ep,
    upperSphere.contact_circle_point(phi),
    lowerSphere.focus,
    ep,
    lowerSphere.contact_circle_point(phi),    
    ]
    ]
    

    key_values.extend( value )

print(key_values)

interpolator = x3d.CoordinateInterpolator(
    key = list(keys),
    keyValue = key_values,
    DEF = 'foci_segments_interpolator'
)

sc.children.append(interpolator)

route1 = x3d.ROUTE(
    fromNode = 'foci_segments_interpolator',
    fromField = 'value_changed',
    toNode = 'foci_segments',
    toField = 'point'
)

sc.children.append(route1)

time_sensor=x3d.TimeSensor(
    DEF = 'timer1',
    cycleInterval = 5.0,
    enabled=True,
    loop = True,
)
sc.children.append(time_sensor)
sc.children.append(
    x3d.ROUTE(
        fromNode="timer1",
        fromField='fraction_changed',
        toNode = 'foci_segments_interpolator',
        toField = 'set_fraction'
    )
)


# In[430]:


get_ipython().run_cell_magic('html', '', '<script type=\'text/javascript\' charset="UTF-8" src=\'https://x3dom.org/release/x3dom-full.debug.js\'> </script> \n<link rel=\'stylesheet\' type=\'text/css\' href=\'https://x3dom.org/release/x3dom.css\'></link> \n<style>\n\n x3d {\n     width:720px;\n     height:720px;\n     border:2px solid black;\n}\n</style>\n')


# In[431]:


x3dnode = """\
<x3d>
%s
</x3d>
<script>
window.x3dom.reload()
</script>
""" % sc.HTML5()

#print(x3dnode)
HTML(x3dnode)


# In[432]:


if (True):
    print(sc.XML())


# In[ ]:





# In[ ]: