SharedProject.sagewsOpen in CoCalc
#Problem 6.5.2
#sources: http://doc.sagemath.org/pdf/en/reference/plot3d/plot3d.pdf, http://doc.sagemath.org/html/en/reference/plotting/sage/plot/point.html
x=var('x')
y=var('y')
#plot(x)

###############################################
#critical points are when (f_x,f_y) = (0,0)
#critical points happen at y=0 and x= +/- sqrt(2)
pt1=point((0,0,0),color='blue', size=10)
pt2=point((sqrt(2),0,4),color='blue',size=10)
pt3=point((-sqrt(2),0,4),color='blue',size=10)
#plot(x+pt1+pt2+pt3)

###############################################

#tangent planes: http://doc.sagemath.org/html/en/reference/plot3d/sage/plot/plot3d/plot3d.html; http://sagemath.wikispaces.com/plot3d; http://omega.albany.edu:8008/mat214dir/Sage/PointsLinesPlanesSage.pdf

def f(x,y):
return -4*x+2*y #equation of the plane

####tangent plane at the point (+/-sqrt(2),0,4)####
plane=plot3d(f, (-2.2,2.2), (-2,2), opacity=0.5)
poin=point((-1,-1,2),color='red',size=10)
#plot(x+plane+p)

###############################################

####fernet frame at the chosen point####
p=vector([-1,-1,2])
q1=vector([4,-2,1])
t=var('t')
N=parametric_plot3d(p+t*q1, (t,0,0.5)) #normal vector to plane

q2=vector([1.5,1,-4])                  #arbitrary vector in the plane
B1=parametric_plot3d(p+t*q2,(t,0,0.5), color='red')

q3=vector([7,(17.5),7])
B2=parametric_plot3d(p+t*q3,(t,0,0.08),color='black') #cross product of normal and chosen vector
#plot(x+poin+B1+B2+N+plane)
###############################################

pr=vector([sqrt(2),0,4])
vec=vector([1,0,4*(1)^2-(1)^4-(0)^2])
v=parametric_plot3d(pr+t*vec,(t,0,1), color='blue') #vector starting from critical point (sqrt(2), 0, 4) going in the principal direction [1,0]

pr1=vector([sqrt(2),0,4])
vec1=vector([2,0,4*(2)^2-(2)^4-(0)^2])
v1=parametric_plot3d(pr1+t*vec1,(t,0,1), color='red') #same as above with prinicpal direction scaled

pr2=vector([sqrt(2),0,4])
vec2=vector([3,0,4*(3)^2-(3)^4-(0)^2])
v2=parametric_plot3d(pr2+t*vec2,(t,0,0.1))            #same as above wih prinicpal direction scaled again

#k=plot3d((cos(t))^2*(16)+(sin(t))^2*2, (r,-1,1),(t,-4,4), color='purple', adaptive =True ) #Euler's curvature formula
#plot(x+pt2+v+v1+v2)
###############################################
#trying to parametrize surface in terms of one variable
#X=parametric_plot3d((t,t,4*t^2-t^4-t^2),(t,-2.2,2.2),color='blue',adaptive=True, opacity=0.5) #results in just a curve

###############################################
def principal_curve(x,y,v,h):
def f(x,y):
return 4*x^2-x^4-y^2
# V=f(sqrt(2),0)
F=f(x,y)
p=vector([x,y,F])
v1=v[0]
v2=v[1]
fv=f(v1,v2)
Pv=vector([v1,v2,fv])
Nv=p+h*Pv
pt=point((Nv[0],Nv[1],f(Nv[0],Nv[1])), color='red')
#V=parametric_plot3d((p+t*Pv),(t,0,0.5), color='red')

return pt
Pv=vector([1,0])
PNT=principal_curve(sqrt(2),0,vec,0.25)

PNT2=principal_curve(0,0,vec,0.25)
PNT3=principal_curve(0.5,0,vec,0.25)
PNT4=principal_curve(0.25,0,vec,0.25)
PNT5=principal_curve(0,0,vec,0.25)
PNT6=principal_curve(-0.5,0,vec,0.25)
PNT7=principal_curve(-0.25,0,vec,0.25)
PNT8=principal_curve(-1,0,vec,0.25)
PNT9=principal_curve(-1.5,0,vec,0.25)
PNT10=principal_curve(-1.25,0,vec,0.25)
plot(x+pt2+PNT+PNT2+PNT3+PNT4+PNT+PNT6+PNT7+PNT8+PNT9+PNT10)

3D rendering not yet implemented