geodesicsParamSurf.sagews

from scipy.integrate import odeint

u, v, t, du, dv = var('u v t du dv')

def fading_line3d(points, rgbcolor1, rgbcolor2, *args, **kwds):
    L = len(points)
    vcolor1 = vector(RDF, rgbcolor1)
    vcolor2 = vector(RDF, rgbcolor2)
    return sum(line3d(points[j:j+2], 
                      rgbcolor = tuple( ((L-j)/L)*vcolor1 + (j/L)*vcolor2 ), 
                      *args, **kwds) 
               for j in srange(L-1))

steps = 100

@interact
def _(x = input_box(3*sin(u)*cos(v), 'x'),
      y = input_box(sin(u)*sin(v), 'y'),
      z = input_box(2*cos(u), 'z'),
      int_u = input_grid(1, 2, default = [[0,pi]], label = 'u -interval'), 
      int_v = input_grid(1, 2, default = [[-pi,pi]], label = 'v -interval'),
      init_point = input_grid(1, 2, default = [[-pi/4,pi/8]], label = 'coordinates of \ninitial point'),
      init_vector = input_grid(1, 2, default = [[1,0]], label = 'coordinates of \ninitial vector'),
      int_s = slider(0, 10, 1/10, 
                           default = pi/2, 
                           label = 'geodesic interval'),
      sliding_color = checkbox(True,'change color along the geodesic')):

    int_u = int_u[0]
    int_v = int_v[0]
    u_0, v_0 = init_point[0]
    V_u, V_v = init_vector[0]
    
    F = vector([x, y, z])

    S_plot = parametric_plot3d( F, 
                                (u, int_u[0], int_u[1]), 
                                (v, int_v[0], int_v[1]))
   
    dFu = F.diff(u)
    dFv = F.diff(v)
    
    Fu = fast_float(dFu, u, v)
    Fv = fast_float(dFv, u, v)
    
    ufunc = function('ufunc', t)
    vfunc = function('vfunc', t)
    
    dFtt = F(u=ufunc, v=vfunc).diff(t, t)
    
    ec1 = dFtt.dot_product(dFu(u=ufunc, v=vfunc))
    ec2 = dFtt.dot_product(dFv(u=ufunc, v=vfunc))
    
    dv, ddv, du, ddu = var('dv, ddv, du, ddu')
    
    diffec1 = ec1.subs_expr(diff(ufunc, t) == du, 
                            diff(ufunc, t, t) == ddu,     
                            diff(vfunc, t) == dv, 
                            diff(vfunc, t, t) == ddv, 
                            ufunc == u, vfunc == v)
    diffec2 = ec2.subs_expr(diff(ufunc, t) == du, 
                            diff(ufunc, t, t) == ddu,     
                            diff(vfunc, t) == dv, 
                            diff(vfunc, t, t) == ddv, 
                            ufunc == u, vfunc == v)
    sols = solve([diffec1 == 0 , diffec2 == 0], ddu, ddv)
    
    ddu_rhs = (sols[0][0]).rhs().full_simplify()
    ddv_rhs = (sols[0][1]).rhs().full_simplify()
        
    ddu_ff = fast_float(ddu_rhs, du, dv, u, v)
    ddv_ff = fast_float(ddv_rhs, du, dv, u, v)
    
    def func(y,t):
        v = list(y)
        return [ddu_ff(*v), ddv_ff(*v), v[0], v[1]]
        
    Point = [u_0, v_0]
    velocity = [V_u, V_v]
    Point = map(float, Point)
    velocity = map(float, velocity)
    
    geo2D_aux = odeint(func,
                       y0 = [velocity[0], velocity[1], Point[0], Point[1]],
                       t = srange(0, int_s, 0.01))
    
    geo3D = [F(u=l,v=r) for [j, k, l, r] in geo2D_aux]
    
    if sliding_color:
        g_plot = fading_line3d(geo3D, rgbcolor1 = (1, 0, 0), rgbcolor2 = (0, 1, 0), thickness=4)
    else:
        g_plot = line3d(geo3D, rgbcolor=(0, 1, 0), thickness=4)
    
    P = F(u=Point[0], v=Point[1])
    P_plot = point3d((P[0], P[1], P[2]), rgbcolor = (0, 0, 0), pointsize = 30)
    V = velocity[0] * Fu(u = Point[0], v = Point[1]) + \
        velocity[1] * Fv(u= Point[0], v = Point[1])
    V_plot = arrow3d(P, P + V, color = 'black')
    
    show(g_plot + S_plot + V_plot + P_plot,aspect_ratio = [1, 1, 1])

Interact: please open in CoCalc