/* define the surface */ // ellipsoid // length parameters #define a=1; #define b=1; #define c=1; // parameterization in terms of generalized coordinates q1 and q2 #declare xs=function(q1,q2) { a*cos(q1)*cos(q2)} #declare ys=function(q1,q2) { b*sin(q1)*cos(q2)} #declare zs=function(q1,q2) { c*sin(q2)} //derivative functions #declare delter=1E-8; // a really small "interval" for approximating derivatives #declare x_1=function(q1,q2) { (xs(q1+delter,q2)-xs(q1,q2))/delter} #declare x_2=function(q1,q2) { (xs(q1,q2+delter)-xs(q1,q2))/delter} #declare y_1=function(q1,q2) { (ys(q1+delter,q2)-ys(q1,q2))/delter} #declare y_2=function(q1,q2) { (ys(q1,q2+delter)-ys(q1,q2))/delter} #declare z_1=function(q1,q2) { (zs(q1+delter,q2)-zs(q1,q2))/delter} #declare z_2=function(q1,q2) { (zs(q1,q2+delter)-zs(q1,q2))/delter} //metric functions #declare g11=function(q1,q2) {x_1(q1,q2)*x_1(q1,q2)+y_1(q1,q2)*y_1(q1,q2)+z_1(q1,q2)*z_1(q1,q2) } #declare g12=function(q1,q2) {x_1(q1,q2)*x_2(q1,q2)+y_1(q1,q2)*y_2(q1,q2)+z_1(q1,q2)*z_2(q1,q2) } #declare g21=function(q1,q2) {x_1(q1,q2)*x_2(q1,q2)+y_1(q1,q2)*y_2(q1,q2)+z_1(q1,q2)*z_2(q1,q2) } #declare g22=function(q1,q2) {x_2(q1,q2)*x_2(q1,q2)+y_2(q1,q2)*y_2(q1,q2)+z_2(q1,q2)*z_2(q1,q2) } //metric derivative functions #declare g11_1=function(q1,q2) { (g11(q1+delter,q2)-g11(q1,q2))/delter} #declare g12_1=function(q1,q2) { (g12(q1+delter,q2)-g12(q1,q2))/delter} #declare g21_1=function(q1,q2) { (g21(q1+delter,q2)-g21(q1,q2))/delter} #declare g22_1=function(q1,q2) { (g22(q1+delter,q2)-g22(q1,q2))/delter} #declare g11_2=function(q1,q2) { (g11(q1,q2+delter)-g11(q1,q2))/delter} #declare g12_2=function(q1,q2) { (g12(q1,q2+delter)-g12(q1,q2))/delter} #declare g21_2=function(q1,q2) { (g21(q1,q2+delter)-g21(q1,q2))/delter} #declare g22_2=function(q1,q2) { (g22(q1,q2+delter)-g22(q1,q2))/delter} //metric inverse #declare gfact=function(q1,q2) { g11(q1,q2)*g22(q1,q2)-g12(q1,q2)*g12(q1,q2)} #declare gin11= function(q1,q2) { g22(q1,q2)/gfact(q1,q2) } #declare gin12= function(q1,q2) { -g12(q1,q2)/gfact(q1,q2) } #declare gin21= function(q1,q2) { -g12(q1,q2)/gfact(q1,q2) } #declare gin22= function(q1,q2) { g11(q1,q2)/gfact(q1,q2) } //Christoffel #declare C111= function(q1,q2) { .5*(gin11(q1,q2)*g11_1(q1,q2)+gin12(q1,q2)*(2*g21_1(q1,q2)-g11_2(q1,q2))) } #declare C112= function(q1,q2) { .5*(gin11(q1,q2)*g11_2(q1,q2)+gin12(q1,q2)*g22_1(q1,q2)) } #declare C121= function(q1,q2) { .5*}