from __future__ import division

import matplotlib as mpl

from mpl_toolkits.mplot3d import Axes3D

import numpy as np

import matplotlib.pyplot as plt

from matplotlib import cm

from matplotlib.patches import FancyArrowPatch

from mpl_toolkits.mplot3d import proj3d

from matplotlib import pylab

from matplotlib.text import Annotation

from sympy.vector import *

from sympy import *

class Arrow3D(FancyArrowPatch):

FancyArrowPatch.draw(self, renderer)

class Annotation3D(Annotation):

Annotation.draw(self, renderer)

def annotate3D(ax, s, *args, **kwargs):

ax.add_artist(tag)

def PlotBasisCartesian(ax, font_size=20):

textcoords='offset points', ha='right',va='bottom')

def plot_arrow(ax,o, v1, e, name,col):

textcoords='offset points', ha='right',va='bottom')

def plot_basis(ax,o,v1,v2,v3,e,col):

plot_arrow(ax,o,v1,e,r'$ \hat{e_1}$',col)

plot_arrow(ax,o,v2,e,r'$ \hat{e_2}$',col)

plot_arrow(ax,o,v3,e,r'$ \hat{e_3}$',col)

def unit_tangent_natural(r,s):

return diff(r,s)

def unit_tangent_time(r,t):

drdt = unit_tangent_natural(r,t)

return drdt / drdt.magnitude()

def tangent_line_natural(r,s,c):

tangent = unit_tangent_natural(r,s)

return r + c * tangent

def tangent_line_natural_at_point(r,s,s0,c,c0):

return tangent_line_natural(r,s,c).subs({s:s0,c:c0})

def tangent_line_time(r,t,c):

tangent = unit_tangent_time(r,t)

return r + c * tangent

def tangent_line_time_at_point(r,t,t0,c, c0):

return tangent_line_time(r,t,c).subs({t:t0,c:c0})

def space_curve(ax, e, r, t, tt, col = 'k'):

ax.plot(fx(tt),fy(tt),fz(tt),color = col)

def space_curve_from_vector_matrix(ax, r, t, tt,col = 'k'):

ax.plot(fx(tt),fy(tt),fz(tt),color = col)

def curve_2d_vector_matrix(ax, r,t,tt,col = 'k'):

ax.plot(fx(tt),fy(tt),color = col)

def curv_vec_natural(r, s):

tangent = unit_tangent_natural(r,s)

return diff(tangent, s)

def curv_vec_from_tangent(tangent, s):

return diff(tangent,s)

def curvature_natural(r,s):

return curv_vec_natural(r,s).magnitude()

def curv_mag(k):

return k.magnitude()

def rad_of_curv_natural(r,s):

return 1/curvature_natural(r,s)

def rad_of_curv_from_curv_vec(k):

return 1/curv_mag(k)

def curv_vec_time(r,t):

return dTdt / drdt_mag

def normal_vec_time(r,t):

k = curve_vec_time(r,t)

return k/k.magnitude()

def principle_normal_natural(r,s):

k = curv_vec_natural(r,s)

return k/k.magnitude()

def principle_normal_from_curv_vec(k):

return k / k.magnitude()

def principle_normal_from_tangent(tangent,s):

k = diff(tangent,s)

return k / k.magnitude()

def curvature_time(r,t):

rprime1 = diff(r,t)

rprime2 = diff(rprime1,t)

return (rprime1.cross(rprime2)).magnitude() / (rprime1.magnitude()**3)

def normal_plane_natural(r, s, s0):

tangent = unit_tangent_natural(r,s)

return (r-r.subs({s:s0})).dot(tangent.subs({s:s0}))

def normal_plane_time(r,t,t0):

tangent = unit_tangent_time(r,t)

return (r-r.subs({t:t0})).dot(tangent.subs({t:t0}))

def osculating_plane_natural(y,r,s):

return (y-r).dot(bi)

def plot_plane(ax,e,r1,tn1,nm1,X,Y,col='b', plane_alpha=0.1):

ax.plot_surface( float(r1.dot(e.i)) + X * float(tn1.dot(e.i)) + Y * float(nm1.dot(e.i)),

float(r1.dot(e.j)) + X * float(tn1.dot(e.j)) + Y * float(nm1.dot(e.j)),

float(r1.dot(e.k))+ X * float(tn1.dot(e.k)) + Y * float(nm1.dot(e.k)), color = col,alpha = plane_alpha)

def plot_plane_from_point_normal(ax, e, P, N):

float(P.dot(e.k))+ X * float(tn1.dot(e.k)) + Y * float(nm1.dot(e.k)), alpha = 0.1)

def binormal_natural(r,s):

return tangent.cross(normal)

def binormal_natural_from_tangent(tangent,s):

norm = diff(tangent,s)

norm = norm/norm.magnitude()

return tangent.cross(normal)

def binormal_from_tangent_normal(tanget,normal):

return tangent.cross(normal)

def binormal_line(r,s,c):

bi = binormal_natural(r,s)

return r + c*bi

def binormal_line_at_point(r,s,s0,c,c0):

return binormal_line(r,s,c).subs({s:s0,c:c0})

def torsion_natural(r,s):

return -dbds.dot(normal)

def torsion_time(r,t):

return num / (den**2)

def torsion_from_tangent_normal(tangent,normal,s):

binormal = tangent.cross(normal)

dbds = diff(binormal, s)

return -dbds.dot(normal)

def torsion_from_normal_binormal(normal,binormal,s):

dbds = diff(binormal, s)

return -dbds.dot(normal)

def dbds_natural(r,s):

return diff(binormal, s)

def torsion_from_dbds_normal(dbds,normal):

return -dbds.dot(normal)

def plot_frenet_frame(ax,o,T,N,e,col):

plot_arrow(ax,o,B,e,r'$ \hat{B}$',col)

def plot_frenet_frame_2(ax, r, tangent, normal, binormal,e,time_point, X, Y):

float(r1.dot(e.k))+ X * float(bn1.dot(e.k)) + Y * float(nm1.dot(e.k)), color = 'r',alpha = 0.1)

def monge_patch_z(x,y,fxy, e):

return x*e.i + y*e.j + fxy * e.k

def space_surface(ax, f, x, y, X,Y,rs,cs):

fxy = lambdify( (x,y), f, "numpy" )

# plot the lambda funcs

ax.plot_surface(X,Y, fxy(X,Y),rstride=rs, cstride=cs,color='g',alpha=1,linewidth=0.0,antialiased=False)

def computePeaks( X, Y ):

return 3*((1-X)**2)*exp(-(X**2) - (Y+1)**2) - 10*(X/5 - X**3 - Y**5)*exp(-X**2-Y**2) - (1/3)*exp(-(X+1)**2 - Y**2);

def surface_normal(r,u,v):

ru_cross_rv = diff(r,u).cross(diff(r,v))

return ru_cross_rv

def surface_unit_normal(r,u,v):

ru_cross_rv = diff(r,u).cross(diff(r,v))

return ru_cross_rv/ru_cross_rv.magnitude()

def compute_E(x,u):

xu = diff(x,u)

return xu.dot(xu)

def compute_F(x,u,v):

xu = diff(x,u)

xv = diff(x,v)

return xu.dot(xv)

def compute_G(x,v):

xv = diff(x,v)

return xv.dot(xv)

def first_fundamental_form(x,u,v):

return EFG

def surface_area_patch( x, u, v, u_min, u_max, v_min, v_max ):

return N(integrate( int_1,(v,v_min,v_max)).doit())

def compute_L_1(x,u,v):

return -xu.dot(nu)

def compute_M_1(x,u,v):

return (S(1)/2)*( -xu.dot(nv)-xv.dot(nu))

def compute_N_1(x,u,v):

return -xv.dot(nv)

def compute_L_2(x,u,v):

n = surface_unit_normal(x,u,v)

xuu = diff(x,u,2)

return xuu.dot(n)

def compute_M_2(x,u,v):

n = surface_unit_normal(x,u,v)

xuv = diff(diff(x,v),u)

return xuv.dot(n)

def compute_N_2(x,u,v):

n = surface_unit_normal(x,u,v)

xvv = diff(diff(x,v),v)

return xvv.dot(n)

def second_fundamental_form(x,u,v):

return LMN

def second_fundamental_form_version2(x,u,v):

return LMN

def quadratic_formula(a,b,c):

return sols

def rect_vector_of_sphere_coords(e, Y):

return Y[0] * cos (Y[2]) * sin(Y[1]) * e.i+ Y[0] * sin (Y[2]) * sin(Y[1]) * e.j + Y[0] * cos(Y[1]) * e.k

def rect_vector_of_cylinder_coords(e, Y):

return Y[0] * cos (Y[1]) * e.i+ Y[0] * sin (Y[1]) * e.j + Y[2] * e.k

def rect_vector_of_parab_cylinder_coords(e, Y):

return (S(1)/2) (Y[0]**2 - Y[1]**2 ) * e.i+ Y[0] * Y[1] * e.j + Y[2] * e.k

def rect_vector_of_paraboloidal_coords(e, Y):

return Y[0]*Y[1]*cos(Y[2]) * e.i+ Y[0]*Y[1]*sin(Y[2]) * e.j + (S(1)/2) (Y[0]**2 - Y[1]**2 )* e.k

def rect_vector_of_elliptic_cylindrical_coords(e, Y):

return cosh(Y[0]) *cos(Y[1])* e.i+ sinh(Y[0]) *sin(Y[1])* e.j + Y[2]* e.k

def rect_vector_of_elliptic_cylindrical_coords(e, Y):

return cosh(Y[0]) *cos(Y[1])* e.i+ sinh(Y[0]) *sin(Y[1])* e.j + Y[2]* e.k

def vector_to_matrix_form(e,v):

return Matrix([v.dot(e.i), v.dot(e.j),v.dot(e.k)])

def basis_vectors(r,u1,u2,u3):

return [drdu1, drdu2, drdu3]

def matrix_from_vector2(e,v1,v2,v3):

return Matrix([[v1.dot(e.i),v1.dot(e.j),v1.dot(e.k)],

[v2.dot(e.i),v2.dot(e.j),v2.dot(e.k)],

[v3.dot(e.i),v3.dot(e.j),v3.dot(e.k)]])

def vector_from_matrix(e,v,i):

return v[0]*e.i + v[1]*e.j + v[2] * e.k

def Jacobian_WRT_coords(V,C):

return V.jacobian(C)

def metric_from_jacobian(A):

return A.T * A

def metric_from_coords(V,C):

return metric_from_jacobian(Jacobian_WRT_coords(V,C))

def christoffel_symbol_2_2(G, p, a, b, X, d, G_inv):

return output

def compute_christoffel_symbols_2(metric, d, basis, metric_inv):

return gamma2

def compute_christoffel_symbols_3(metric, d, basis, metric_inv):

return gamma2

from sympy import Array

def covariant_derivative_covariant_vector(T, Y, gamma2, d):

return d_r_T_i

def covariant_derivative_contravariant_vector(T, Y, gamma2, d):