import numpy as np
import numpy.linalg
import math
from make import rotate_extrude
from geometry import *
from optics import *
from solids import build_8inch_pmt
from transform import rotate, make_rotation_matrix
import sno_pmt_locations

def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20):
    '''Compute the 2D coordinates of the profile of one side of a 
    SNO-style acrylic vessel.  The center of the sphere is at (0,0), with
    the neck extending along the positive y direction.
    
       sphere_radius: Radius of spherical part of profile
       neck_radius: Radius of neck part
       neck_top: y coordinate of top of neck
       angle_step: angular step size (radius) between points on the sphere

       Returns: Tuple of x and y coordinate numpy arrays.
    '''
    if neck_radius >= sphere_radius:
        raise ValueError('neck_radius must be less than sphere_radius')

    intersect_height = (sphere_radius**2 - neck_radius**2)**0.5
    max_angle = math.atan2(intersect_height, neck_radius)

    if neck_top < intersect_height:
        raise ValueError('neck_top must be greater than the y-value where the sphere and cylinder intersect')

    # Start with point at bottom
    angles = np.arange(-math.pi/2, max_angle, angle_step)
    x = list(np.cos(angles) * sphere_radius)
    y = list(np.sin(angles) * sphere_radius)
    x[0] = 0.0 # Round-off error might make cos(-pi/2) not exactly zero
    
    # Neck intersection point
    x.append(neck_radius)
    y.append(intersect_height)

    # Top of neck
    x.append(neck_radius)
    y.append(neck_top)
    
    # Top of neck on axis
    x.append(0.0)
    y.append(neck_top)
    
    return x, y

##### SNO Parts

angle_step = math.pi/20
av_outside_profile = sno_vessel(sphere_radius=6.0604, neck_radius=0.79375, 
                                neck_top=10.50, angle_step=angle_step)
# For simplicity, cap the top of the AV with acrylic
av_inside_profile  = sno_vessel(sphere_radius=6.0053, neck_radius=0.72898,
                                neck_top=10.00, angle_step=angle_step)

av_outside_mesh = rotate_extrude(av_outside_profile[0], av_outside_profile[1],
                                 angle_step)
av_outside_mesh.vertices = rotate(av_outside_mesh.vertices, np.pi/2, (-1,0,0))

av_inside_mesh = rotate_extrude(av_inside_profile[0], av_inside_profile[1],
                                 angle_step)
av_inside_mesh.vertices = rotate(av_inside_mesh.vertices, np.pi/2, (-1,0,0))

def build_sno():
    pmt = build_8inch_pmt()
    geo = Geometry()
    geo.add_solid(Solid(av_outside_mesh, glass, water, color=0x99FFFFFF))
    geo.add_solid(Solid(av_inside_mesh, water, glass, color=0x990000FF))
    geo.pmtids = []
    for x,y,z in zip(sno_pmt_locations.x, sno_pmt_locations.y,
                     sno_pmt_locations.z):
        direction = numpy.array((-x,-y,-z))
        direction /= numpy.linalg.norm(direction)
        y_axis = np.array((0.0,1.0,0.0))
        axis = np.cross(direction, y_axis)
        angle = np.arccos(np.dot(y_axis, direction))

        rotation = make_rotation_matrix(angle, axis)
        displacement = (x/1000.0, y/1000.0, z/1000.0)
        geo.pmtids.append(geo.add_solid(pmt, rotation, displacement))
        
    return geo