out ff;
%off nat;
for all x let cos(x)**2+sin(x)**2=1;
for all x,y let sqrt(x)*sqrt(y)=sqrt(x*y);

%-------------- 3D parabolic field ------------
v(1):=-sin(s); 
v(2):=cos(s);
v(3):=1;

%v(1):=s; 
%v(2):=s**3;
%v(3):=1;

dv:=sqrt(v(1)**2 + v(2)**2 + v(3)**2);
%-------------- first derivative of field ------------
vv(1):=df(v(1),s);
vv(2):=df(v(2),s);
vv(3):=df(v(3),s);
dvv:=sqrt(vv(1)**2 + vv(2)**2 + vv(3)**2);

%-------------- second derivative of field ------------
vvv(1):=df(vv(1),s);
vvv(2):=df(vv(2),s);
vvv(3):=df(vv(3),s);
%-------------- tangent unit field ------------
tv(1):=v(1)/dv;
tv(2):=v(2)/dv;
tv(3):=v(3)/dv;
%-------------- binormal field ------------
VxVV(1):=v(2)*vv(3) - vv(2)*v(3);
VxVV(2):=v(3)*vv(1) - vv(3)*v(1);
VxVV(3):=v(1)*vv(2) - vv(1)*v(2);
dVxVV:=sqrt(VxVV(1)**2 + VxVV(2)**2 + VxVV(3)**2);
bv(1):=VxVV(1)/dVxVV;
bv(2):=VxVV(2)/dVxVV;
bv(3):=VxVV(3)/dVxVV;
%-------------- normal field ------------
pnv(1):=bv(2)*tv(3) - tv(2)*bv(3);
pnv(2):=bv(3)*tv(1) - tv(3)*bv(1);
pnv(3):=bv(1)*tv(2) - tv(1)*bv(2);
dpnv:=sqrt(pnv(1)**2 + pnv(2)**2 + pnv(3)**2);
nv(1):=pnv(1)/dpnv;
nv(2):=pnv(2)/dpnv;
nv(3):=pnv(3)/dpnv;

%-------------- curvature ------------
cur:=dVxVV/(dV**3);

tormat:=mat
(
(v(1),v(2),v(3)),
(vv(1),vv(2),vv(3)),
(vvv(1),vvv(2),vvv(3))
);

det(tormat);

tor:=det(tormat)/(dVxVV**2);



%-------------- frene equations ------------
(1/dv)*df(tV(1),S) - cur*nv(1);
(1/dv)*df(tV(2),S) - cur*nv(2);
(1/dv)*df(tV(3),S) - cur*nv(3);

(1/dv)*df(nV(1),S) + cur*tv(1) - tor*bv(1);
(1/dv)*df(nV(2),S) + cur*tv(2) - tor*bv(2);
(1/dv)*df(nV(3),S) + cur*tv(3) - tor*bv(3);


(1/dv)*df(bV(1),S) + cur*nv(1);
(1/dv)*df(bV(2),S) + cur*nv(2);
(1/dv)*df(bV(3),S) + cur*nv(3);



shut ff;
bye;
end;