# Parametric surfaces in 3-dimensional space and tangent planes
# Author: Thomas Krainer (Penn State Altoona)
# Prior Version: 08/15/2011
# Version: 02/11/2016
#
# Minor changes.

html('<h1>Parametric Surfaces and Tangent Planes</h1>')
html('<p>Specify a surface in ${\\mathbb R}^3$ by parametric equations \
$$ S : \\left\\{\\begin{aligned} x &= f_1(u,v) \\\ y &= f_2(u,v) \\\ z &= f_3(u,v) \\end{aligned} \\right. \\textrm{ for } (u,v) \\in [u_{\\min},u_{\\max}] \\times [v_{\\min},v_{\\max}]. $$')
html('Via the slider bars you can specify parameter values $(u_0,v_0)$ which in turn determine a reference point $(x_0,y_0,z_0)$ on $S$. The program computes the tangent plane and a normal vector to $S$ at that point, and they are shown in the plot along with the surface.</p>')

@interact # (layout={'top' : [['xcoord'],['ycoord'],['zcoord'],['uborders'],['vborders'],['uslider'],['vslider'],['Plane_size'],['Normalvector_size']],'bottom' : [['show_point','show_tangentplane','show_normalvector'],['show_computations','ndigits']]})
def param_surface(xcoord=input_box(default='(3+cos(u/2)*sin(v)-sin(u/2)*sin(2*v))*cos(u)',type=str,label='$x =$'),\
                  ycoord=input_box(default='(3+cos(u/2)*sin(v)-sin(u/2)*sin(2*v))*sin(u)',type=str,label='$y =$'),\
                  zcoord=input_box(default='sin(u/2)*sin(v)+cos(u/2)*sin(2*v)',type=str,label='$z =$'),\
                  uborders=input_box('[0,2*pi]',type=str,label='$u \\in [u_{\\min},u_{\\max}] =$',width=20),\
                  vborders=input_box(default='[0,2*pi]',type=str,label='$v \\in [v_{\\min},v_{\\max}] =$',width=20),\
                  uslider=slider(0,100,1,default=25,label='$u_0 \\in [u_{\\min},u_{\\max}] :$',display_value=False),\
                  vslider=slider(0,100,1,default=60,label='$v_0 \\in [v_{\\min},v_{\\max}] :$',display_value=False),\
                  Plane_size=slider(0,100,1,label='Tangent plane size:',default=40,display_value=False),\
                  Normalvector_size=slider(0,100,1,label='Normal vector size:',default=10,display_value=False),\
                  show_point=("Show point:",True),\
                  show_tangentplane=("$\\quad$ Show tangent plane:",True),\
                  show_normalvector=("$\\quad$ Show normal vector:",True),\
                  show_computations=("Show computations:",False)):

    # This variable will contain the computations output
    comp_output=''

    ndigits=5     # Numerical accuracy

    # Process input data
    [umin,umax]=sage_eval(uborders)
    [vmin,vmax]=sage_eval(vborders)
    f_1(u,v)=SR(xcoord)
    f_2(u,v)=SR(ycoord)
    f_3(u,v)=SR(zcoord)

    # Parameter values (u_0,v_0) as determined by slider input
    u_0=umin+(uslider/100)*(umax-umin)
    v_0=vmin+(vslider/100)*(vmax-vmin)

    # Compute data needed for tangent plane and normal vector
    uderivative(u,v) = (diff(f_1,u),diff(f_2,u),diff(f_3,u)) # direction vector corresponding to u-variable
    vderivative(u,v) = (diff(f_1,v),diff(f_2,v),diff(f_3,v)) # direction vector corresponding to v-variable
    normalvector = uderivative.cross_product(vderivative) # normal vector to the tangent plane

    if (round(normalvector(u_0,v_0).norm(),ndigits)==0):

       comp_output += 'The surface is not immersed at the designated values for $(u_0,v_0)$.<br>Tangent plane and normal vector are undefined!<br><br>'
       show_tangentplane = False
       show_normalvector = False

    else:

        normednormalvector = (1/normalvector.norm())*normalvector

        # Function that returns a customized tangent plane plot
        def tangentplane(u,v,size):
            vector_1 = (1/uderivative.norm())*uderivative
            vector_aux = vderivative - (vderivative.dot_product(vector_1))*vector_1
            vector_2 = (1/vector_aux.norm())*vector_aux
            # r and s parametrize a quadratic tangent plane piece of size size-by-size in the plot
            r,s = var('r,s')
            return parametric_plot3d(vector((f_1(u,v),f_2(u,v),f_3(u,v)))+r*vector_1(u,v)+s*vector_2(u,v),
            (r,-size/2,size/2),(s,-size/2,size/2),color='green',opacity=0.7,plot_points=[100,100])


    # Assemble the plot according to parameters

    the_plot=parametric_plot3d((f_1(u,v),f_2(u,v),f_3(u,v)),(u,umin,umax),(v,vmin,vmax),aspect_ratio=[1,1,1],plot_points=[100,100])
    plot_size = (vector(the_plot.bounding_box()[1])-vector(the_plot.bounding_box()[0])).norm()

    if (show_tangentplane == True):
        the_plot += tangentplane(u_0,v_0,plot_size*(Plane_size/100))
    if (show_point == True):
        the_plot += point3d((f_1(u_0,v_0),f_2(u_0,v_0),f_3(u_0,v_0)),size=10,color='red')
    if (show_normalvector == True):
        the_plot += arrow(vector((f_1(u_0,v_0),f_2(u_0,v_0),f_3(u_0,v_0))),(vector((f_1(u_0,v_0),f_2(u_0,v_0),f_3(u_0,v_0))) + plot_size*(Normalvector_size/100)*normednormalvector(u_0,v_0)),color='red')


    # Assemble computations output

    if (show_computations == True):

        comp_output += 'Surface:'
        comp_output += '$$ S : \\left\\{\\begin{aligned} x&= %s \\\ y&= %s \\\ z&= %s \\end{aligned}\\right. \\qquad (u,v) \\in [%s,%s] \\times [%s,%s] $$'%(latex(f_1(u,v)),latex(f_2(u,v)),latex(f_3(u,v)),latex(round(umin,ndigits)),latex(round(umax,ndigits)),latex(round(vmin,ndigits)),latex(round(vmax,ndigits)))

        if (show_point == True):
           comp_output += 'Reference point:'
           comp_output += '$$ \\begin{aligned} (u_0,v_0) &= (%s,%s) \\\ (x_0,y_0,z_0) &= (%s,%s,%s) \\end{aligned} $$'%(latex(round(u_0,ndigits)),latex(round(v_0,ndigits)),latex(round(f_1(u_0,v_0),ndigits)),latex(round(f_2(u_0,v_0),ndigits)),latex(round(f_3(u_0,v_0),ndigits)))

        if (show_tangentplane == True):
            comp_output += 'Tangent plane:'
            plane(x,y,z)= round(normalvector(u_0,v_0)[0],ndigits)*(x-round(f_1(u_0,v_0),ndigits)) + round(normalvector(u_0,v_0)[1],ndigits)*(y-round(f_2(u_0,v_0),ndigits)) + round(normalvector(u_0,v_0)[2],ndigits)*(z-round(f_3(u_0,v_0),ndigits))
            comp_output += '$$'+latex(plane(x,y,z)-plane(0,0,0))+'= '+latex(-round(plane(0,0,0),ndigits))+'$$'

        if (show_normalvector == True):
            comp_output += 'Unit normal vector:'
            comp_output += '$$ \\langle %s,%s,%s \\rangle $$'%(latex(round(normednormalvector(u_0,v_0)[0],ndigits)),latex(round(normednormalvector(u_0,v_0)[1],ndigits)),latex(round(normednormalvector(u_0,v_0)[2],ndigits)))


    # Display plot and computations

    html(comp_output)
    show(the_plot)