Add SHEETWARP case.

src/warp.c

@@ -562,6 +562,32 @@ void windup_o2d( struct warpstuff *ws, Point *out, CONST Point *in ) {
     dvtfmpoint( out, &tout, ws->d2o );
 }
 
+void sheetwarp( struct warpstuff *ws, Point *xyzout, CONST Point *xyz )
+{
+  /*
+   * v(r) =  (vmax/rcore)*r  r in 0..rcore
+   *         vmax             r > rcore
+   * period(r) = 2pi/(vmax/rcore)  r<rcore  [Call this T0.]
+   *             2pi/(vmax/r) = (r/rcore)*T0
+   * theta(t/T0,r) = (t/T0)*T0*(r < rcore ? 1 : rcore/r)
+   * So rotate in the XZ plane about (0,0,0) by theta(t/T0,r)
+   */
+   float sx = xyz->x[0] / ws->rcore;
+   float sz = xyz->x[2] / ws->routercore;
+
+   xyzout->x[0] = xyz->x[0];
+   xyzout->x[1] = xyz->x[1] + sx * ws->sheetampl * expf( -sx*sx ) * expf( -sz*sz );
+   xyzout->x[2] = xyz->x[2];
+}
+
+void sheetwarp_o2d( struct warpstuff *ws, Point *out, CONST Point *in ) {
+    Point tin, tout;
+
+    dvtfmpoint( &tin, in, ws->o2d );
+    sheetwarp( ws, &tout, &tin );
+    dvtfmpoint( out, &tout, ws->d2o );
+}
+
 void extrap( struct warpstuff *ws, Point *out, CONST Point *in, CONST struct speck *sp, int deg )
 {
     int i;
@@ -754,6 +780,17 @@ void warpspecks( struct warpstuff *ws,
     }
     break;
 
+  case SHEETWARP:
+    for(i = 0; i < n; i++) {
+	if(ws->has_o2d)
+	    sheetwarp_o2d( ws, &sp->p, &osp->p );
+	else
+	    sheetwarp( ws, &sp->p, &osp->p );
+	osp = NextSpeck(osp, osl, 1);
+	sp = NextSpeck(sp, sl, 1);
+    }
+    break;
+
   case GALAXYRIDE:
     if(ws->coef0 >= 0 && ws->coef0 + 8 <= SPECKMAXVAL(osl)) {
 	for(i = 0; i < n; i++) {