File:Illustration of principal curvatures on an oblate spheroid (ellipsoid of revolution).webm

// This work is licensed PD.
#version 3.7;
global_settings { assumed_gamma 1.2 } 

#include "colors.inc"
#include "textures.inc"
#include "shapes.inc"

// camera
camera {
   location <2.5,1.5,-3>*1.25
   angle 45
   right x*image_width/image_height
   look_at 0
}

// light
#declare Dist=50.0;
light_source {
  <25, 25, -25> color White
  fade_distance Dist fade_power 2
  // Uncomment the area lights only if you've got lots of time.
  area_light <-40, 0, -40>, <40, 0, 40>, 4, 4
  adaptive 1
  jitter
}

// rotated ellipse parameters
#declare a  = 1.0;
#declare _e = 0.4; // e�=a�-b�, _e=e/a
#declare b  = a*sqrt( 1 - _e*_e );

#macro ellipsoid_of_rev(col)
sphere {
  0, 1
  no_shadow pigment { color col }

  scale <a,b,a>
}
#end

#macro lon_curvature_radius(p0,_radius,_center)
  //#local p  = vrotate(p0, <0,VAngleD(p0*(x+z), p0.x),0 // untested
  #local p  = p0;

  #local aq = a*a*a*a;
  #local bq = b*b*b*b;
  #local s  = aq*p.y*p.y + bq*p.x*p.x;
  #declare _radius = 1/(aq*bq) * sqrt(s*s*s);

  #local esqr = a*a - b*b; //_e*_e*a*a;
  #declare _center = <esqr * p.x*p.x*p.x / aq, -esqr * p.y*p.y*p.y / bq, p.z>;
#end

#macro lat_curvature_radius(p0,_radius,_center)
  #local phi = VAngle(n0,n0*(x+z)); // geodetic lat at p0
  
  #local lonr = 0;
  #local lonc = 0;
  lon_curvature_radius(p0, lonr, lonc);

  #declare _radius = a / sqrt(1-(_e*a)*(_e*a)*sin(phi)*sin(phi));
  #declare _center = p0 + _radius*vnormalize(lonc-p0);
#end

#macro tplane(col,scl,p0,n0)
box {
  <-1,-1e-4,-1>, <1,1e-4,1>
  no_shadow pigment { color col }

  scale <scl,1,scl>
  Reorient_Trans(y,n0)
  translate p0
}
#end

#macro tex_xy_peel(col)
texture {
  pigment {
    gradient <1,1,0>
    color_map {
      [ 0.0 color rgbf<1,1,1,1> ]
      [ 0.5 color rgbf<1,1,1,1> ]
      [ 0.5 color col ]
      [ 1.0 color col ]
    }
    scale 5e-2
  }
}
#end

// axes, foci
#declare o_axes = union {
  cylinder { -x*a,x*a,4e-3 no_shadow pigment { color White } }
  cylinder { -z*a,z*a,4e-3 no_shadow pigment { color White } }
  cylinder { -y*b,y*b,4e-3 no_shadow pigment { color White } }
  sphere { -_e*a*x,1e-2 pigment { color White } }
  sphere { +_e*a*x,1e-2 pigment { color White } }
  sphere { -_e*a*z,1e-2 pigment { color White } }
  sphere { +_e*a*z,1e-2 pigment { color White } }
};

    

//ellipsoid
#declare o_ellps = object {
  ellipsoid_of_rev(color Grey transmit .5)
  scale 1+5.05e-3 // torus minor radii that vizualize curvature circles
};

//tangential plane in point p0
#declare slat = clock;   // latitude in degrees (this is _not_ the geodetic lat at p0)
#declare slon = 0;       // longitude in degrees, untested for != 0
#declare sn0  = vrotate(vrotate(x, slat*z), slon*-y);
#declare slat = radians(slat);
#declare slon = radians(slon);

#declare p0   = 0;
#declare n0   = 0;
IsectN(0, sn0, ellipsoid_of_rev(Clear), p0, n0)

#declare p0_lat       = VAngle( n0,  n0*(x+z));
#declare p0_slat      = VAngle(sn0, sn0*(x+z));
#declare n0_xz_length = ( abs(sin(slat))<=1e-16 ? b*b/a :
                          ( abs(cos(slat))<=1e-16 ? b :
                            ( vlength(p0) * sin(p0_slat) / sin(pi-p0_lat) )));
#declare p0_lat  = sgn(sin(slat)) * p0_lat;  // geodetic lat at p0
#declare p0_slat = sgn(sin(slat)) * p0_slat; // geocentric lat

#declare o_tan_plane = union {
  object {
    tplane(color DimGray transmit .8, .33, p0 + 5.1e-3*n0, n0)
  }

  cylinder { // inward normal of tplane, cut at xz plane (equatorial plane)
    p0, p0 - n0*n0_xz_length, 5e-3
    pigment { color Gray }
    no_shadow
  }
};

//vertical plane in point p0, plane with minimum curvature
#declare nu = 0;
#declare nu_center = 0;
lat_curvature_radius(p0,nu,nu_center);

//#declare n0v = vcross(n0,-z);
#declare n0v = vcross(nu_center-p0, -z);

#declare o_lat_plane = union {
  intersection {
    object { ellipsoid_of_rev(Clear) }
    object { tplane(color Red transmit .6, 10, p0, n0v) }
  }
  union {
    sphere { // center
      0, 1.66e-2
      pigment { color .5*White+.5*Red }
    }
    cylinder { // radius
      //-nu*n0, nu*n0, 5e-3
      -nu*vnormalize(nu_center-p0),
      +nu*vnormalize(nu_center-p0), 5.4e-3
      tex_xy_peel(color Red transmit -max(-1, min(0, 1.66*sin(slat))))
    }
    torus { // circle
      nu, 5e-3
      pigment { color Red }

      Reorient_Trans(y, n0v)
    }
    difference {
      intersection {
        object { sphere { 0, nu } }
        object { tplane(color .5*White+.5*Red transmit .8, 10, 0, y) }

        Reorient_Trans(y, n0v)
        //translate p0-nu*n0
        translate nu_center
      }
      ellipsoid_of_rev(Clear)
      //translate -1*(p0-nu*n0)
      translate -1*nu_center
    }
    no_shadow
    //translate p0-nu*n0
    translate nu_center
  }
};

//meridional plane in point p0, plane with maximum curvature
#declare rho = 0;
#declare rho_center = 0;
lon_curvature_radius(p0,rho,rho_center);

//#declare n0m = vcross(n0,vcross(n0,-z));
#declare n0m = vcross(rho_center-p0,vcross(rho_center-p0, -z));

#declare o_lon_plane = union {
  intersection {
    object { ellipsoid_of_rev(Clear) }
    object { tplane(color Blue transmit .7, 10, p0, n0m) }

  }
  union {
    sphere { // center
      0, 1.66e-2
      pigment { color .5*White+.5*Blue }
    }
    cylinder { // radius
      //-rho*n0, rho*n0, 5e-3
      -rho*vnormalize((rho_center-p0)),
      +rho*vnormalize((rho_center-p0)), 5.2e-3
      pigment { color Blue transmit -max(-1, min(0, 1.66*sin(slat))) }
    }
    torus { // circle
      rho, 5e-3
      pigment { color Blue }
  
      Reorient_Trans(y, n0m)
    }
    difference {
      intersection {
        object { sphere { 0, rho } }
        object { tplane(color .5*White+.5*Blue transmit .8, 10, 0, y) }

        Reorient_Trans(y, n0m)
        //translate p0-rho*n0
        translate rho_center
      }
      ellipsoid_of_rev(Clear)
      //translate -1*(p0-rho*n0)
      translate -1*rho_center
    }
    no_shadow
    //translate p0-rho*n0
    translate rho_center
  }
};

// text elements
#declare o_texts = union {
  union {
    text {
      ttf "timrom.ttf" concat("curvature radius nu: ", str(nu,8,5)) .0125, 0
      no_shadow pigment { Red }

      translate 0*y-12*x
    }
    text {
      ttf "timrom.ttf" concat("curvature 1/nu: ", str(1/nu,8,5)) .0125, 0
      no_shadow pigment { Red }

      translate -1*y-10*x
    }
    text {
      ttf "timrom.ttf" concat("a =", str(a,5,4), "  b = ", str(b,5,4)) .0125, 0
      no_shadow pigment { Gray10 }
  
      scale .8
      translate -2*y-7*x
    }
    text {
      ttf "timrom.ttf" concat("phi = ", str(degrees(p0_lat),8,4)) .0125, 0
      no_shadow pigment { Gray10 }
  
      scale .8
      translate -3*y-4.5*x
    }
    text {
      ttf "timrom.ttf" concat("frame = ", str(frame_number,-4,0)) .0125, 0
      no_shadow pigment { Gray10 }
  
      scale .8
      translate -4*y-4*x rotate 1*y
    }

    scale .1 rotate 90*x
    translate -1.25*z-.5*y
  }
  union {
    text {
      ttf "timrom.ttf" concat("curvature radius rho: ", str(rho,8,5)) .0125, 0
      no_shadow pigment { Blue }
  
      translate 0*y-6.66*x
    }
    text {
      ttf "timrom.ttf" concat("curvature 1/rho: ", str(1/rho,8,5)) .0125, 0
      no_shadow pigment { Blue }
  
      translate -1*y-4.65*x
    }
  
    scale .09-0.01*y rotate -90*x
    translate -1.5*z
    rotate -90*y rotate 90*x+5*z
  }
}

// scene
union {
  object { o_axes }
  object { o_ellps }
  object { o_tan_plane }
  object { o_lon_plane }
  object { o_lat_plane }
  
  rotate -16*sin(slat)*y
  translate .1*y + min(0, sin(slat))*.2*y
}

object { o_texts }