A geodesic from the origin to any point

This program will draw a curve from the orgin to any point.  The height of a curve in Nil space is equal to the area contained by the projection of the curve and the line connecting the ends of the curve. So the projection of the geodesic is an arc of the circle where the area of the circular segment is equal to the final z cordinate. This function finds the circle given the circular area (the z cordinate) and the chord length (c = Sqrt[xcor^2 +ycor^2]). That plots that arc and lifts it into Nil space. However this function does not work for points on the z axis. Which is why the second function is included.

Clear[t,x,y,r,f,g,h,a,b,c, points];
(*Input a point*)
xcor = 1;
ycor = 1;
zcor = 10;

(*User should not change any of the code below this point*)
c = Sqrt[xcor^2 +ycor^2];
area = zcor; 
case = (xcor^2 +ycor^2)*Pi/8;
If[case<zcor, g = r^2*Pi-(1/2*r^2*(2*ArcSin[c/(2*r)]-Sin[2*ArcSin[c/(2*r)]]))-area, g = (1/2*r^2*(2*ArcSin[c/(2*r)]-Sin[2*ArcSin[c/(2*r)]]))-area];
rad =FindRoot[g, {r, c}]
If[case <zcor,  τ = 2*Pi-2*ArcSin[c/(2*r)]/.rad,τ = 2*ArcSin[c/(2*r)]/.rad ];
If[xcor <0, θ =τ/2+ ArcTan[ycor/xcor], θ = Pi+ τ/2 + ArcTan[ycor/xcor] ];  
tend = τ 
Clear[t,x,y,z,dx,dy,dz,f];
oldx=r *-Sin[t]/.rad;
oldy=r*(-Cos[t]+1)/.rad;
x =(Cos[θ]oldx - Sin[θ]oldy);
y = (Sin[θ] oldx + Cos[θ]oldy);
dz=-.5*(x*D[y,t]-D[x,t]*y);
dz=TrigReduce[dz];
z=Integrate[dz,t];
temp =z;
t = 0;
d = temp;
Clear[t]
z = z -d;
Print["The projection on the xy plane"];
ParametricPlot[{x,y},{t,0,tend}]
Print["The curve in Nil space"];
ParametricPlot3D[{x,y,z},{t,0,tend}]
points = Table[{x,y,z}, {t, 0, tend, .01}];

Print["See how the tangent plane changes on the curve"];
Animate[{f,g,h}= points[[i]];
u={1,0,g/2};
v={0,1,-f/2};
e=Cross[u,v];
{a,b,c}=e;
d=e.{f,g,h};
Show[ParametricPlot3D[{x,y,z},{t,0,tend},BoxRatios->{1,1,2}],Plot3D[(-a/c)*xx+(-b/c)*yy+d/c,{xx,-(xcor+2*r),xcor+ 2*r},{yy,-(ycor+2*r),ycor +2*r},PlotRange -> {0, zcor+1},Mesh->{{{f, Thick}}, {{g,Thick}}},BoxRatios->{1,1,2},ClippingStyle->None]], {i, 1, tend/.01,1}]/.rad


dx=D[x,t];
dy=D[y,t];
f=dx^2+dy^2+dz^2;
f=Sqrt[TrigReduce[f]];
Print["The length of the curve is:"];
Print[NIntegrate[f,{t,0.00,tend},WorkingPrecision->15]];
t=tend;
Print["The final position is :" ]
Print[{x,y,z}]

Clear[t,x,y,f,g,h,a,b,c];

Geodesic to a point on the z axis

Clear[t];
zcor = 12

α = Sqrt[zcor/Pi];

x=α*Sin[t];
y=(α)*(-Cos[t]+1);
dz=.5*(x*D[y,t]-D[x,t]*y);
dz=TrigReduce[dz];
z=Integrate[dz,t];
Print["The projected ellipse"];
ParametricPlot[{x,y},{t,0,2*Pi}]
Print["The curve in Nil space"];
ParametricPlot3D[{x,y,z},{t,0,2*Pi}]

points = Table[{x,y,z}, {t, 0, 2*Pi, .01}];

Print["See how the tangent plane changes on the curve"];
Animate[{f,g,h}= points[[i]];
u={1,0,-g/2};
v={0,1,f/2};
e=Cross[u,v];
{a,b,c}=e;
d=e.{f,g,h};
Show[ParametricPlot3D[{x,y,z},{t,0,2*Pi}, BoxRatios->{1,1,1}], Plot3D[(-a/c)*xx+(-b/c)*yy+d/c,{xx,-10,10},{yy,-10,10},PlotRange -> {0, zcor}, Mesh->{{{f, Thick}}, {{g,Thick}}},BoxRatios->{1,1,1},ClippingStyle->None]],{i, 1,2*Pi/.01,1}]

dx=D[x,t];
dy=D[y,t];
f=dx^2+dy^2+dz^2;
f=Sqrt[TrigReduce[f]];
Print["The length of the curve is:"];
Print[NIntegrate[f,{t,0.00,2.*Pi},WorkingPrecision->15]];
t=2*Pi;
Print["The final position is :" ]
Print[z]

|

|