/* define the surface */ // ellipsoid // length parameters #declare a=1.5; #declare b=1.5; #declare c=2; // parameterization in terms of generalized coordinates q1 and q2 /* ellipsoid */ /*#declare xs=function(q1,q2) { a*cos(q1)*cos(q2)} #declare zs=function(q1,q2) { b*sin(q1)*cos(q2)} #declare ys=function(q1,q2) { c*sin(q2)} #declare fns = function(xp,yp,zp) {xp*xp/(a*a)+zp*zp/(b*b)+yp*yp/(c*c) -1} //isosurface */ /* hyperbolic paraboloid */ #declare xs=function(q1,q2) { a*q1} #declare ys=function(q1,q2) { (q1*q1-q2*q2)/c} #declare zs=function(q1,q2) { b*q2} #declare fns = function(xp,yp,zp) {xp*xp/(a*a)-zp*zp/(b*b)-c*yp} //isosurface //derivative functions #declare delter=1E-6; // 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*(gin11(q1,q2)*g11_2(q1,q2)+gin12(q1,q2)*g22_1(q1,q2))} #declare C122= function(q1,q2) { .5*(gin11(q1,q2)*(2*g12_2(q1,q2)-g22_1(q1,q2))+gin12(q1,q2)*g22_2(q1,q2))} #declare C211= function(q1,q2) { .5*(gin21(q1,q2)*g11_1(q1,q2)+gin22(q1,q2)*(2*g21_1(q1,q2)-g11_2(q1,q2)))} #declare C212= function(q1,q2) { .5*(gin21(q1,q2)*g11_2(q1,q2)+gin22(q1,q2)*g22_1(q1,q2))} #declare C221= function(q1,q2) { .5*(gin21(q1,q2)*g11_2(q1,q2)+gin22(q1,q2)*g22_1(q1,q2))} #declare C222= function(q1,q2) { .5*(gin21(q1,q2)*(2*g12_2(q1,q2)-g22_1(q1,q2))+gin22(q1,q2)*g22_2(q1,q2))}