found in ian@chiark:things/moebius.old-before-cvs

diff --git a/.depend b/.depend
new file mode 100644 (file)
index 0000000..bda22ed
--- /dev/null
+++ b/.depend
@@ -0,0 +1,10 @@
+dualx11.o : dualx11.cc dualx11.hh output.hh library.hh parameter.hh 
+graphics.o : graphics.cc library.hh graphics.hh output.hh 
+library.o : library.cc library.hh transforms.hh 
+main.o : main.cc output.hh library.hh x11.hh dualx11.hh postscript.hh moebius.hh \
+  transforms.hh graphics.hh parameter.hh 
+moebius.o : moebius.cc library.hh moebius.hh 
+parameter.o : parameter.cc parameter.hh 
+postscript.o : postscript.cc postscript.hh output.hh library.hh 
+transforms.o : transforms.cc transforms.hh library.hh 
+x11.o : x11.cc x11.hh output.hh library.hh 

diff --git a/Makefile b/Makefile
new file mode 100644 (file)
index 0000000..80111df
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,18 @@
+C++FLAGS= -Wall -O2
+CPPFLAGS= -I/usr/include
+LDFLAGS= 
+LIBS= 
+
+OBJS=  dualx11.o x11.o main.o parameter.o graphics.o library.o \
+       transforms.o postscript.o moebius.o
+
+a.out: $(OBJS)
+       $(LINK.cc) -o a.out $(OBJS) -lm -L/usr/X11R6/lib -lX11 -lg++
+
+depend:
+       $(C++) $(CPPFLAGS) -E -MM *.cc >.depend
+
+clean:
+       rm -f $(OBJS) a.out
+
+include .depend

diff --git a/Makefile.bak b/Makefile.bak
new file mode 100644 (file)
index 0000000..ee52661
--- /dev/null
+++ b/Makefile.bak
@@ -0,0 +1,6 @@
+
+OBJS= main.o moebius.o output.o library.o
+
+a.out: $(OBJS)
+       g++ -o a.out $(OBJS)
+

diff --git a/dualx11.cc b/dualx11.cc
new file mode 100644 (file)
index 0000000..4bfc39c
--- /dev/null
+++ b/dualx11.cc
@@ -0,0 +1,79 @@
+/*
+ * X11 functions
+ */
+
+#include "dualx11.hh"
+#include "parameter.hh"
+
+DualX11Output::DualX11Output() {
+  display= XOpenDisplay(0);
+  Colormap cmap= DefaultColormap(display,DefaultScreen(display));
+
+  XAllocColorCells(display,cmap,
+                   False, // contig
+                   planemasks,
+                   2, // nplanes
+                   &pixelbase,
+                   1); // ncolors
+  XColor colours[4], *cp= colours;
+  for (int i=0; i<2; i++)
+    for (int j=0; j<2; j++, cp++) {
+      cp->pixel= pixelbase + (i? planemasks[0] :0) + (j? planemasks[1] :0);
+      cp->red= i?65535:0;  cp->green= j?/*65535*/32000:0;  cp->blue= 0;
+      cp->flags= DoRed|DoGreen|DoBlue;
+    }
+  XStoreColors(display,cmap,colours,4);
+  
+  window= XCreateSimpleWindow(display,
+                              DefaultRootWindow(display),
+                              0,0, 500,500, 0,0, pixelbase);
+
+  XGCValues gcvalues;
+  for (i=0; i<2; i++) {
+    gcvalues.plane_mask= planemasks[i];
+    // First, the fabric
+    gcvalues.function= GXclear;
+    fabric[i]= XCreateGC(display,window,
+                         GCFunction|GCPlaneMask,
+                         &gcvalues);
+    // Then, the mesh
+    gcvalues.function= GXset;
+    mesh[i]= XCreateGC(display,window,
+                       GCFunction|GCPlaneMask,
+                       &gcvalues);
+  }
+  XSelectInput(display,window,0);
+  XMapWindow(display,window);
+  XFlush(display);
+}
+
+void DualX11Output::startimage() {
+  XClearWindow(display,window);
+}
+
+void DualX11Output::endimage() {
+  XFlush(display);
+}
+
+DualX11Output::~DualX11Output() {
+  XCloseDisplay(display);
+}
+
+void DualX11Output::drawcell(const Point* list, int n) {
+  static Parameter<double> eyeseparation("eyesep",
+                                         "Distance from projection eye to origin", 
+                                         0.3, .1, 0., 100.);
+  for (int i=0; i<2; i++) {
+    Point::seteyex(eyeseparation*(i-0.5));
+    XPoint xp[n+1];
+    for (int j=0; j<n; j++) {
+      Onscreen here= Onscreen(list[j]);
+      xp[j].x= (int)((here.x+1.0)*250.0);
+      xp[j].y= (int)((-here.y+1.0)*250.0);
+    }
+    XFillPolygon(display,window,fabric[i],xp,n,Nonconvex,CoordModeOrigin);
+    xp[n]= xp[0];
+    XDrawLines(display,window,mesh[i],xp,n+1,CoordModeOrigin);
+  }
+  Point::seteyex(0);
+}

diff --git a/dualx11.hh b/dualx11.hh
new file mode 100644 (file)
index 0000000..042bc11
--- /dev/null
+++ b/dualx11.hh
@@ -0,0 +1,27 @@
+/*
+ * X11 functions
+ */
+
+#ifndef DUALX11_HH
+#define DUALX11_HH
+
+#include <stdlib.h>
+#include <X11/Xlib.h>
+
+#include "output.hh"
+
+struct DualX11Output : Output {
+  DualX11Output();
+  ~DualX11Output();
+  void drawcell(const Point*, int);
+  void startimage();
+  void endimage();
+private:
+  Display *display;
+  Window window;
+  GC fabric[2];
+  GC mesh[2];
+  unsigned long planemasks[2], pixelbase;
+};
+
+#endif

diff --git a/fit b/fit
new file mode 100644 (file)
index 0000000..f7bd11a
--- /dev/null
+++ b/fit
@@ -0,0 +1 @@
+epsffit -c -s 50 450 550 750 t

diff --git a/graphics.cc b/graphics.cc
new file mode 100644 (file)
index 0000000..5fdfb25
--- /dev/null
+++ b/graphics.cc
@@ -0,0 +1,112 @@
+/*
+ * Graphics library
+ */
+
+#include <iostream.h>
+
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "library.hh"
+#include "graphics.hh"
+#include "output.hh"
+
+static int checkneededline(Onscreen a, Onscreen b) {
+  double xd= a.x - b.x;
+  double yd= a.y - b.y;
+  int n= (int)ceil(sqrt(xd*xd+yd*yd)*100.0);
+  return n>0 ? n : 1;
+}
+
+void Cell::display(Output& o) {
+  int nn[4];
+  int totalnn= 0;
+  for (int i=0; i<4; i++) {
+    nn[i]= checkneededline(p[i],p[(i+1)&3]);
+    totalnn+= nn[i];
+  }
+  Point *array= new Point[totalnn];
+  Point *inarrayp= array;
+  for (i=0; i<4; i++) {
+    for (int a=0; a<nn[i]; a++) {
+      double fp= (double)a         / nn[i];
+      double fn= (double)(nn[i]-a) / nn[i];
+      *inarrayp++= p[i]*fn + p[(i+1)&3]*fp;
+    }
+  }
+  o.drawcell(array,totalnn);
+  delete[] array;
+}                          
+
+SortedCellList::Fragment* SortedCellList::insertfragment(int at) {
+  if (used==size) {
+    size+=5; size<<=1;
+    fragments= (Fragment**)realloc(fragments,sizeof(Fragment*)*size);
+    if (!fragments) { perror("realloc"); exit(1); }
+  }
+  memmove(fragments+at+1, fragments+at, sizeof(Fragment*)*(used-at));
+  used++;
+  Fragment *nf= new Fragment;
+  fragments[at]= nf;
+  return nf;
+}
+
+void SortedCellList::insert(Cell *it) {
+  double index= it->index();
+  for (int fragment= used-1;
+       fragment >= 0 && index < fragments[fragment]->entries[0].index;
+       fragment--);
+  if (fragment < 0) {
+    Fragment *nf= insertfragment(0);
+    nf->used= 1;
+    nf->entries[0].index= index;
+    nf->entries[0].cell= it;
+    return;
+  }
+  for (int entry= 0;
+       entry < fragments[fragment]->used &&
+       index > fragments[fragment]->entries[entry].index;
+       entry++);
+  if (fragments[fragment]->used >= fragmentsize) {
+    Fragment *nf= insertfragment(fragment+1);
+    nf->used= fragmentsize>>1;
+    fragments[fragment]->used -= nf->used;
+    memcpy(nf->entries,
+           fragments[fragment]->entries + fragments[fragment]->used,
+           nf->used*sizeof(Entry));
+    if (entry >= fragments[fragment]->used) {
+      entry-= fragments[fragment]->used;
+      fragment++;
+    }
+  }
+  memmove(fragments[fragment]->entries + entry+1,
+          fragments[fragment]->entries + entry,
+          sizeof(Entry) * (fragments[fragment]->used - entry));
+  fragments[fragment]->entries[entry].index= index;
+  fragments[fragment]->entries[entry].cell= it;
+  fragments[fragment]->used++;
+}
+
+SortedCellList::~SortedCellList() {
+  for (int fragment=0; fragment<used; fragment++) {
+    for (int entry=0; entry<fragments[fragment]->used; entry++) {
+      delete fragments[fragment]->entries[entry].cell;
+    }
+    delete fragments[fragment];
+  }
+  free(fragments);
+}
+
+void SortedCellList::display(Output& o) {
+  o.startimage();
+  for (int fragment=0; fragment<used; fragment++) {
+    for (int entry=0; entry<fragments[fragment]->used; entry++) {
+      if (Point::indexvisible(fragments[fragment]->entries[entry].index)) {
+        fragments[fragment]->entries[entry].cell->display(o);
+      }
+    }
+  }
+  o.endimage();
+}

diff --git a/graphics.hh b/graphics.hh
new file mode 100644 (file)
index 0000000..bf7c3d2
--- /dev/null
+++ b/graphics.hh
@@ -0,0 +1,40 @@
+/*
+ * Graphics library
+ */
+
+#ifndef GRAPHICS_HH
+#define GRAPHICS_HH
+
+#include "library.hh"
+#include "output.hh"
+
+class Cell {
+  void display(Output& o);
+  friend class SortedCellList;
+public:
+  Point p[4];
+  double index() { return p[0].index(); }
+};
+
+#define fragmentsize 10
+class SortedCellList {
+  struct Entry {
+    double index;
+    Cell *cell;
+  };
+  struct Fragment {
+    int used;
+    Entry entries[fragmentsize];
+  };
+  Fragment **fragments;
+  int size,used;
+  
+  Fragment *insertfragment(int at);
+public:
+  SortedCellList() { fragments=0; size=used=0; }
+  ~SortedCellList();
+  void insert(Cell*);
+  void display(Output& o);
+};
+
+#endif

diff --git a/ins b/ins
new file mode 100644 (file)
index 0000000..70ba461
--- /dev/null
+++ b/ins
@@ -0,0 +1,8 @@
+tear=1
+twist=-1
+theta=1.57
+bulge=1.5
+sheary=3
+un=20
+tn=35
+x11

diff --git a/library.cc b/library.cc
new file mode 100644 (file)
index 0000000..60f8446
--- /dev/null
+++ b/library.cc
@@ -0,0 +1,50 @@
+/*
+ * Points library
+ */
+
+#include <stdlib.h>
+
+#include "library.hh"
+#include "transforms.hh"
+
+double Point::planedistance= 1.0;
+double Point::eyedistance= 1.0;
+double Point::cutoffdistance= 10.0;
+double Point::eyex= 0.0;
+
+TransformList Point::usertransforms;
+TransformList Point::povtransform;
+
+Point TransformList::operator()(Point it) {
+  for (int i=0; i<used; i++) { it= (*a[i])(it); }
+  return it;
+}
+
+void TransformList::append(Transform *it) {
+  if (used >= size) {
+    size+=5; size<<=1;
+    a= (Transform**)realloc(a,sizeof(Transform*)*size);
+    if (!a) { perror("realloc"); exit(1); }
+  }
+  a[used++]= it;
+}
+
+void TransformList::clearcontents() {
+  for (int i=0; i<used; i++) delete a[i];
+}
+  
+Point::operator Onscreen() const {
+  Point it= transformed();
+  double factor= eyedistance / (eyedistance + planedistance - it[2]);
+  return Onscreen(it[0] * factor + eyex * (1.0-factor), it[1] * factor);
+}
+
+void Point::setobserver(double theta, double eta, double planedist, double eyedist,
+                        double cutoffdist) {
+  planedistance= planedist;
+  eyedistance= eyedist;
+  cutoffdistance= cutoffdist;
+  povtransform.reset();
+  if (theta != 0.0) povtransform.append(new XZRotationTransform(theta));
+  if (eta != 0.0) povtransform.append(new YZRotationTransform(eta));
+}

diff --git a/library.hh b/library.hh
new file mode 100644 (file)
index 0000000..bb226dd
--- /dev/null
+++ b/library.hh
@@ -0,0 +1,72 @@
+/*
+ * Points library
+ */
+
+#ifndef POINTS_HH
+#define POINTS_HH
+
+#include <stdlib.h>
+
+/*
+ * Coordinate system:
+ *
+ *          Y |
+ *            |
+ *            |
+ *            /-------- X
+ *           /
+ *          /
+ *       Z (coming out of the monitor)
+ */
+
+struct Onscreen {
+  double x, y;
+  Onscreen(){};
+  Onscreen(double ix, double iy) { x=ix; y=iy; }
+};
+
+class Point;
+
+class Transform {
+public:
+  virtual Point operator()(Point)= 0;
+};
+
+class TransformList {
+  Transform **a;
+  int size,used;
+  void clearcontents();
+ public:
+  TransformList() { size=used=0; a=0; }
+  ~TransformList() { clearcontents(); free(a); }
+  void reset() { clearcontents(); used=0; }
+  Point operator()(Point);
+  void append(Transform*);
+};
+
+class Point {
+  static TransformList usertransforms;
+  static TransformList povtransform;
+  static double planedistance, eyedistance, cutoffdistance, eyex;
+  double xyz[3];
+public:
+  Point(){}
+  Point(double x, double y, double z) { xyz[0]=x; xyz[1]=y; xyz[2]=z; }
+  Point operator+(Point r) { return Point(xyz[0]+r[0], xyz[1]+r[1], xyz[2]+r[2]); }
+  Point operator*(double f) { return Point(xyz[0]*f, xyz[1]*f, xyz[2]*f); }
+  double& operator[](int i) { return xyz[i]; }
+
+  operator Onscreen() const;
+  Point transformed() const { return povtransform(usertransforms(*this)); }
+  double index() const { return transformed()[2]; }
+
+  static int indexvisible(double index) { return index < cutoffdistance; }
+  static void appendtransform(Transform *t) { usertransforms.append(t); }
+  static void cleartransform() { usertransforms.reset(); }
+  static void setobserver(double theta, double eta=0.0,
+                          double planedist=1.0, double eyedist=1.0,
+                          double cutoff=10.0);
+  static void seteyex(double n) { eyex=n; }
+};
+
+#endif

diff --git a/ll -tr b/ll -tr
new file mode 100644 (file)
index 0000000..338407e
--- /dev/null
+++ b/ll -tr
@@ -0,0 +1,29 @@
+total 112
+-rw-rw-r--  1 ian ian   82 Jul 17  1993 Makefile.bak
+-rw-rw-r--  1 ian ian  317 Jul 17  1993 main.cc%
+-rw-rw-r--  1 ian ian  645 Jul 17  1993 graphics.hh
+-rw-rw-r--  1 ian ian 1988 Jul 18  1993 transforms.hh
+-rw-rw-r--  1 ian ian   62 Jul 18  1993 ins
+-rw-rw-r--  1 ian ian 2119 Jul 18  1993 transforms.cc
+-rw-rw-r--  1 ian ian 1559 Jul 19  1993 moebius.hh
+-rw-rw-r--  1 ian ian 2079 Jul 20  1993 moebius.cc
+-rw-rw-r--  1 ian ian   31 Jul 20  1993 fit
+-rw-rw-r--  1 ian ian  245 Sep 11  1993 output.hh
+-rw-rw-r--  1 ian ian  337 Sep 11  1993 x11.hh
+-rw-rw-r--  1 ian ian  302 Sep 11  1993 postscript.hh
+-rw-rw-r--  1 ian ian 1104 Sep 11  1993 x11.cc
+-rw-rw-r--  1 ian ian 3068 Sep 11  1993 graphics.cc
+-rw-rw-r--  1 ian ian 1102 Sep 11  1993 postscript.cc
+-rw-rw-r--  1 ian ian 1788 Sep 11  1993 library.hh
+-rw-rw-r--  1 ian ian 1291 Sep 11  1993 library.cc
+-rw-rw-r--  1 ian ian 5488 Sep 11  1993 main.cc
+-rw-rw-r--  1 ian ian  405 Sep 11  1993 dualx11.hh
+-rw-rw-r--  1 ian ian  553 Sep 11  1993 .depend
+-rw-rw-r--  1 ian ian 1280 Sep 11  1993 parameter.hh
+-rw-rw-r--  1 ian ian  526 Sep 11  1993 parameter.cc
+-rw-rw-r--  1 ian ian 2299 Sep 11  1993 dualx11.cc
+-rw-rw-r--  1 ian ian  338 Aug 31  1996 Makefile
+drwxrwsr-x 60 ian ian 4096 Nov 16 11:06 ../
+drwxrwsr-x  7 ian ian 4096 Mar 16 23:50 .git/
+-rw-rw-r--  1 ian ian    0 Mar 16 23:50 ll -tr
+drwxrwsr-x  3 ian ian 4096 Mar 16 23:50 ./

diff --git a/main.cc b/main.cc
new file mode 100644 (file)
index 0000000..7007194
--- /dev/null
+++ b/main.cc
@@ -0,0 +1,152 @@
+/*
+ * Moebius main program.
+ */
+
+#include <math.h>
+#include <iostream.h>
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "output.hh"
+#include "x11.hh"
+#include "dualx11.hh"
+#include "postscript.hh"
+#include "library.hh"
+#include "moebius.hh"
+#include "transforms.hh"
+#include "graphics.hh"
+#include "parameter.hh"
+
+Parameter<int> tn("tn", "Number of sections around `strip'", 30, 1, 1, 500);
+Parameter<int> un("un", "Number of sections across `strip'", 15, 1, 1, 500);
+Parameter<double> theta("theta", "Angle of model rotation",
+                        /*-PI/2.*/0, PI/8., -PI*2., PI*2.);
+Parameter<double> eta("eta", "Angle of model elevation", 0., PI/8., -PI*2., PI*2.);
+Parameter<double> planedistance("plane", "Distance from projection plane to origin", 
+                                4., .25, 0., 100.);
+Parameter<double> eyedistance("eye", "Distance from projection eye to origin", 
+                              1.5, .25, 0., 100.);
+Parameter<double> cutoffdistance("cutoff", "Distance from projection cutoff to origin", 
+                                 10., .25, -10., 100.);
+Parameter<double> width("width", "Width of the `strip' before transformation", 
+                        .4, .1, 0., 1.);
+Parameter<double> thickness("thick", "Thickness of the `pancake'", 1., .1, 0., 1.);
+Parameter<double> bottomportion("bp", "Proportion in the bottom half", 
+                                .7, .1, 0., 1.);
+//Parameter<double> tear("tear", "Angle factor in tear-open", 0/*1.1*/, .1, -5., 5.);
+//Parameter<double> teary("teary", "Axis y coord in tear-open", 1.3, .1, -10., 10.);
+//Parameter<double> twist("twist", "Angle per length in twist", 0/*.8*/, .25, -5., 5.);
+//Parameter<double> twisty("twisty", "Axis y coord in twist", .85, .1, -10., 10.);
+//Parameter<double> bulge("bulge", "Amount of bulge", 0/*1.*/, .25, -15., 15.);
+//Parameter<double> bulgesd("bulgesd", "S.D. of bulge", 1., .2, .000001, 20.);
+//Parameter<double> bulgey("bulgey", "Axis y coord in bulge", .85, .1, -10., 10.);
+//Parameter<double> shearx("shearx", "Amount of shear (x)", 0., .1, -15., 15.);
+//Parameter<double> sheary("sheary", "Amount of shear (y)", 0/*1.*/, .1, -15., 15.);
+//Parameter<double> shearsdx("shearsdx", "S.D. of shear (x)", .5, .1, -15., 15.);
+//Parameter<double> shearsdy("shearsdy", "S.D. of shear (y)", .5, .1, -15., 15.);
+//Parameter<double> shearey("shearey", "Centre of shear (y)", 1.3, .1, -10., 10.);
+
+SortedCellList *list=0;
+X11Output *x11= 0;
+DualX11Output *dualx11= 0;
+
+void generate() {
+  if (list) delete list;
+  // MoebiusStrip strip(width);
+  MoebiusEnfoldment strip(thickness,bottomportion);
+  list= new SortedCellList;
+  double t0= 0.0;
+  for (int ti=0; ti<tn; ti++) {
+    double t1= (double)(ti+1)/tn;
+    double u0= 0.0;
+    for (int ui=0; ui<un; ui++) {
+      double u1= (double)(ui+1)/un;
+      Cell *cell= new Cell;
+//u0=u1=tear;
+      cell->p[0]= strip.middlepoint(t0,u0);
+//fprintf(stderr,"t %7.4f u %7.4f  %7.4f,%7.4f,%7.4f\n",
+//        t0,u0, cell->p[0][0], cell->p[0][1], cell->p[0][2]);
+      cell->p[1]= strip.middlepoint(t0,u1);
+      cell->p[2]= strip.middlepoint(t1,u1);
+      cell->p[3]= strip.middlepoint(t1,u0);
+      list->insert(cell);
+      u0= u1;
+    }
+    t0= t1;
+  }
+}
+
+int main(int argc, char **argv) {
+  char buf[100];
+
+  for (;;) {
+    Point::setobserver(theta,eta,planedistance,eyedistance,cutoffdistance);
+    Point::cleartransform();
+//    if (tear != 0.0) Point::appendtransform(new TearOpenTransform(tear,teary));
+//    if (twist != 0.0) Point::appendtransform(new TwistTransform(twist,twisty));
+//    if (bulge != 0.0) Point::appendtransform(new BulgeTransform(bulge,bulgesd,bulgey));
+//    if (shearx != 0.0 || sheary != 0.0)
+//      Point::appendtransform(new ShearTransform(Point(0.0,shearey,0.0),
+//                                                Point(shearx,sheary,0.0),
+//                                                Point(shearsdx,shearsdy,1e20)));
+    generate();
+    if (x11) list->display(*x11);
+    if (dualx11) list->display(*dualx11);
+    for (;;) {
+      cin.get(buf,100,'\n');
+      char c;
+      if (!cin.get(c) || c!='\n') {
+        cerr << "error reading command input, giving up\n";
+        exit(1);
+      }
+      char *equals= strchr(buf,'=');
+      if (equals) {
+        *equals++= 0;
+        AnyParameter *param= AnyParameter::find(buf);
+        if (param) {
+          *param= atof(equals);
+          break;
+        } else {
+          cerr << "parameter `" << buf << "' not found\n";
+        }
+      } else if (!strncmp(buf,"postscript ",11)) {
+        const char *filename= buf+11;
+        ofstream file(filename);
+        if (file) {
+          PostScriptOutput psout(file);
+          list->display(psout);
+        } else {
+          cerr << "Failed to open `" << filename << "'\n";
+          continue;
+        }
+        cerr << "PostScript written to `" << filename <<"'\n";
+      } else if (!strcmp(buf,"x11")) {
+        if (x11) {
+          delete x11;
+          x11= 0;
+        } else {
+          x11= new X11Output;
+        }
+        break;
+      } else if (!strcmp(buf,"dualx11")) {
+        if (dualx11) {
+          delete dualx11;
+          dualx11= 0;
+        } else {
+          dualx11= new DualX11Output;
+        }
+        break;
+      } else if (!strcmp(buf,"quit")) {
+        exit(0);
+      } else if (!strcmp(buf,"list")) {
+        AnyParameter::list();
+      } else if (!*buf) {
+        break;
+      } else {
+        cerr << "command not understood\n";
+      }
+    }
+  }
+}
+

diff --git a/main.cc% b/main.cc%
new file mode 100644 (file)
index 0000000..8324362
--- /dev/null
+++ b/main.cc%
@@ -0,0 +1,19 @@
+/*
+ * Moebius main program.
+ */
+
+#include <math.h>
+
+#include "output.hh"
+#include "library.hh"
+#include "graphics.hh"
+#include "moebius.hh"
+
+Moebius strip;
+PostScriptOutput output;
+
+int main() {
+  for (double t=-PI*2.0; t<= PI*2.0; t+=0.1) {
+    output.drawedgeline(strip.edgepoint(t), strip.edgepoint(t+0.1));
+  }
+}

diff --git a/moebius.cc b/moebius.cc
new file mode 100644 (file)
index 0000000..d68f4af
--- /dev/null
+++ b/moebius.cc
@@ -0,0 +1,76 @@
+/*
+ * Equation for a Moebius strip
+ */
+
+#include <math.h>
+#include <stdio.h>
+
+#include "library.hh"
+#include "moebius.hh"
+
+Point MoebiusStrip::edgepoint(double t) {
+  double theta= (t-0.5)*4.0*PI;
+  double r= 1.0 - halfgap*(1.0 - sin(theta/2.0));
+  return Point(r*sin(theta),
+               -r*cos(theta),
+               halfbreadth*cos(theta/2.0));
+}
+
+Point MoebiusStrip::middlepoint(double t, double u) {
+  return edgepoint(t*0.5)*u + edgepoint((t+1)*0.5)*(1.0-u);
+}
+
+class Bezier {
+  double e,f,g,h;
+public:
+  Bezier(double x0, double x1, double dx0, double dx1);
+  double operator()(double t) { return h + t*(g + t*(f + t*e)); }
+  void debug();
+};
+
+Bezier::Bezier(double x0, double x1, double dx0, double dx1) {
+//  e= 2*x0 - 2*x1 + dx0 + dx1;
+//  f= x1 - dx0 - x0 - e;
+  h= x0;
+  g= dx0;
+  e= g + 2*h + dx1 - 2*x1;
+  f= x1 - e - g - h;
+}
+
+void Bezier::debug() {
+  fprintf(stderr,"bz e %7.4f f %7.4f g %7.4f h %7.4f\n",e,f,g,h);
+}
+
+// The first end is at [sin(theta/2),-cos(theta/2),0]
+// The second end is at [-theta/pi,0,sin(theta)]
+// The first end is oriented towards [0,cos(theta),sin(theta)]
+// The second end is oriented towards [0,-1,0]
+
+Point MoebiusEnfoldment::edgepoint(double t) {
+  double theta= t*2.0*PI;
+  return Point(sin(theta),cos(theta),0);
+}
+
+Point MoebiusEnfoldment::middlepoint(double t, double u) {
+  if (t > bottomportion) {
+    t -= bottomportion;
+    t /= (1.0 - bottomportion);
+    double sizehere= sqrt(1-t*t);
+    return Point((u*2.0-1.0) * sizehere,
+                 t,
+                 sizehere * thickness * sin(u*2.0*PI));
+  } else {
+    t /= bottomportion;
+    double theta= (.5-u)*2*PI;
+    Bezier bx(sin(theta*.5), -theta/PI, 0, 0);
+    double ypushiness= (1-cos(theta))*2.5+1;
+//        double fudge= (PI*sin(theta*.5)*cos(theta*.5))*(.5-u)*(.5-u)*4;
+    double fudge= (.5-u)*(.5-u)*4*cos(theta*.5);
+    Bezier by(-cos(theta*.5), 0,
+              cos(theta)*ypushiness + fudge*ypushiness,
+              ypushiness);
+//by.debug();
+    Bezier bz(0, sin(theta), sin(theta), 0);
+    return Point( bx(t), by(t), thickness * bz(t) );
+  }
+}    

diff --git a/moebius.hh b/moebius.hh
new file mode 100644 (file)
index 0000000..e834d3b
--- /dev/null
+++ b/moebius.hh
@@ -0,0 +1,56 @@
+/*
+ * Moebius strip
+ */
+
+#ifndef MOEBIUS_HH
+#define MOEBIUS_HH
+
+class Surface {
+public:
+  // t and u vary from 0 to 1.0
+  virtual Point middlepoint(double t, double u) =0;
+};
+
+class Edge {
+public:
+  // t varies from 0 to 1.0
+  virtual Point edgepoint(double t) =0;
+};
+
+class Moebius : public Surface, Edge {
+public:
+  virtual Point edgepoint(double t) =0;
+  virtual Point middlepoint(double t, double u) =0;
+};
+
+class MoebiusStrip : public Moebius {
+  // Strip is `lying' in the XY plane. The crossing is at the
+  // maximal Y value, where the edge which has X*Z positive
+  // has smaller Y at the point where X=Z=0 than the `other' edge,
+  // which has X*Z negative.  The flat part of the strip
+  // is at minimal Y value.  Here one edge has Z positive and
+  // the other Z negative.
+  // X and Y range from -1.0 to 1.0;
+  // Z ranges from -halfbreadth to +halfbreadth.
+  double halfgap;     // 1/2 the width of the strip at the crossing point
+  double halfbreadth; // 1/2 the width of the strip at the flat part
+public:
+  MoebiusStrip() { halfgap= 0.15; halfbreadth= 0.15; }
+  MoebiusStrip(double b) { halfgap= halfbreadth= b/2.0; }
+  MoebiusStrip(double hg, double hb) { halfgap= hg/2.0; halfbreadth= hb/2.0; }
+    
+  Point edgepoint(double t);
+  Point middlepoint(double t, double u);
+};
+
+class MoebiusEnfoldment : public Moebius {
+  double thickness;
+  double bottomportion;
+public:
+  MoebiusEnfoldment(double t=.35, double bp=.5) { thickness= t; bottomportion= bp; }
+    
+  Point edgepoint(double t);
+  Point middlepoint(double t, double u);
+};
+
+#endif

diff --git a/output.hh b/output.hh
new file mode 100644 (file)
index 0000000..92cc33b
--- /dev/null
+++ b/output.hh
@@ -0,0 +1,17 @@
+/*
+ * Output functions
+ */
+
+#ifndef OUTPUT_HH
+#define OUTPUT_HH
+
+#include "library.hh"
+
+struct Output {
+  virtual ~Output(){};
+  virtual void drawcell(const Point*, int) =0;
+  virtual void startimage(){};
+  virtual void endimage(){};
+};
+
+#endif

diff --git a/parameter.cc b/parameter.cc
new file mode 100644 (file)
index 0000000..0eb6df6
--- /dev/null
+++ b/parameter.cc
@@ -0,0 +1,27 @@
+/*
+ * Parameters library
+ */
+
+#include "parameter.hh"
+
+AnyParameter::AnyParameter(const char *n, const char *d) {
+  name= n; description= d;
+  next= first;
+  first= this;
+}
+
+AnyParameter* AnyParameter::find(const char *n) {
+  for (AnyParameter* search= first;
+       search && strcmp(search->name,n);
+       search= search->next);
+  return search;
+}
+
+void AnyParameter::list() {
+  for (AnyParameter* search= first;
+       search;
+       search= search->next)
+    search->rangecheck();
+}
+
+AnyParameter* AnyParameter::first= 0;

diff --git a/parameter.hh b/parameter.hh
new file mode 100644 (file)
index 0000000..07b0137
--- /dev/null
+++ b/parameter.hh
@@ -0,0 +1,53 @@
+/*
+ * Parameters library
+ */
+
+#ifndef PARAMETER_HH
+#define PARAMETER_HH
+
+#include <iostream.h>
+
+class AnyParameter {
+  AnyParameter *next;
+  static AnyParameter *first;
+protected:
+  const char *name;
+  const char *description;
+  virtual void rangecheck() =0;
+public:
+  AnyParameter(const char *n, const char *d);
+  virtual void operator ++() =0;
+  virtual void operator --() =0;
+  virtual void operator =(double) =0;
+  static AnyParameter *find(const char *n);
+  static void list();
+};
+
+template<class T> class Parameter : AnyParameter {
+  T value, delta, min, max;
+  void rangecheck();
+public:
+  Parameter(const char *n, const char *d, T i, T de, T mi, T ma);
+  operator T (){ return value; }
+  void operator ++(){ value+= delta; rangecheck(); }
+  void operator --(){ value-= delta; rangecheck(); }
+  void operator =(double v) { value= (T)v; rangecheck(); }
+};
+
+template<class T> Parameter<T>::Parameter
+    (const char *n, const char *d, T i, T de, T mi, T ma)
+: AnyParameter(n,d) {
+  value=i; delta=de; min=mi; max=ma;
+}
+
+template<class T> void Parameter<T>::rangecheck() {
+  cerr << name << " ";
+  if (value<min) {
+    value= min;  cerr << "underflowed; ";
+  } else if (value>max) {
+    value= max;  cerr << "overflowed; ";
+  }
+  cerr << "set to " << value << "\n";
+}
+
+#endif

diff --git a/postscript.cc b/postscript.cc
new file mode 100644 (file)
index 0000000..3fbe4e7
--- /dev/null
+++ b/postscript.cc
@@ -0,0 +1,47 @@
+/*
+ * PostScript output
+ */
+
+#include <iostream.h>
+#include "postscript.hh"
+
+PostScriptOutput::PostScriptOutput(ofstream& f)
+: file(f) {
+
+  file <<
+    "%!PS-Adobe-2.0 EPSF-2.0\n"
+    "%%Creator: moebius\n"
+    "%%BoundingBox: -500 -500 500 500\n"
+    "%%Pages: 1\n"
+    "%%DocumentFonts: \n"
+    "%%EndComments\n"
+    "    1 setlinecap 1 setlinejoin 0.002 setlinewidth\n"
+    "    500 500 scale\n"
+    "%%EndProlog\n"
+    "%%Page: 1 0\n"
+    "    save\n";
+}
+
+PostScriptOutput::~PostScriptOutput() {
+  file <<
+    "    restore\n"
+    "%%Trailer\n"
+    "%%EOF\n";
+}
+
+void PostScriptOutput::drawcell(const Point* list, int n) {
+  Onscreen p[4];
+  for (int i=0; i<4; i++) p[i]= Onscreen(list[i]);
+  docell(p,n,"1 setgray fill");
+  docell(p,n,"0 setgray stroke");
+}
+
+void PostScriptOutput::docell(const Onscreen* list, int n, const char *what) {
+  file << "    newpath\n";
+  file << "        " << list->x << " " << list->y << " moveto\n";
+  list++;
+  for (int i=1; i<n; i++, list++) {
+    file << "        " << list->x << " " << list->y << " lineto\n";
+  }
+  file << "      closepath  " << what << "\n";
+}

diff --git a/postscript.hh b/postscript.hh
new file mode 100644 (file)
index 0000000..d385513
--- /dev/null
+++ b/postscript.hh
@@ -0,0 +1,17 @@
+/*
+ * PostScript output
+ */
+
+#include <fstream.h>
+
+#include "output.hh"
+
+struct PostScriptOutput : Output {
+  PostScriptOutput(ofstream&);
+  ~PostScriptOutput();
+  void drawcell(const Point*, int);
+private:
+  ofstream& file;
+  void docell(const Onscreen*, int, const char* what);
+  void preamble();
+};

diff --git a/transforms.cc b/transforms.cc
new file mode 100644 (file)
index 0000000..d5b6cd6
--- /dev/null
+++ b/transforms.cc
@@ -0,0 +1,81 @@
+/*
+ * Transformations
+ */
+
+#include <math.h>
+#include <stdio.h>
+
+#include "transforms.hh"
+
+RotationTransform::RotationTransform(double theta) {
+  s= sin(theta); c= cos(theta);
+}
+
+Point XZRotationTransform::operator()(Point in) {
+  return Point(in[0]*c + in[2]*s,
+               in[1],
+               in[2]*c - in[0]*s);
+}
+
+Point YZRotationTransform::operator()(Point in) {
+  return Point(in[0],
+               in[1]*c + in[2]*s,
+               in[2]*c - in[1]*s);
+}
+
+TearOpenTransform::TearOpenTransform(double amt, double apy) {
+  factor= pow(10.0,amt);
+  apexy= apy;
+}
+
+Point TearOpenTransform::operator()(Point in) {
+  double w= apexy - in[1];
+  double r= sqrt(in[2]*in[2] + w*w);
+  double theta= r>0.0 ? factor*atan2(in[2],w) : 0.0;
+  double ny= apexy - r*cos(theta);
+  double nz= r*sin(theta);
+  return Point(in[0],ny,nz);
+}
+
+BulgeTransform::BulgeTransform(double amt, double stddev, double apy) {
+  amount=amt; apexy=apy;
+  variance= stddev*stddev;
+}
+
+Point BulgeTransform::operator()(Point in) {
+  double w= in[1] - apexy;
+  double r2= w*w + in[2]*in[2];
+  double nror= 1.0 + amount*exp(-r2/variance);
+  double ny= w*nror + apexy;
+  return Point(in[0],ny,in[2]);
+}
+
+Point ShearTransform::operator()(Point in) {
+  double dx= (in[0] - epicentre[0]) / deviation[0];
+  double dy= (in[1] - epicentre[1]) / deviation[1];
+  double dz= (in[2] - epicentre[2]) / deviation[2];
+  double r2= dx*dx + dy*dy + dz*dz;
+  double amount= exp(-r2);
+  return in + magnitude*amount;
+}
+
+Point TwistTransform::operator()(Point in) {
+  double w= in[1] - apexy;
+  double theta= factor*in[2];
+  double c= cos(theta), s= sin(theta);
+  double nx= in[0]*c - w*s;
+  double ny= w*c + in[0]*s + apexy;
+  return Point(nx,ny,in[2]);
+}
+
+Point MagicTransform::operator()(Point in) {
+  double nz= in[2]*zfactor;
+  double r= sqrt(in[0]*in[0] + in[1]*in[1]);
+  double nx= xfactor * (r - xkonstant);
+  double nr= yoffset + yfactor*in[1];
+  double ny2=  nr*nr- nz*nz;
+//fprintf(stderr,"M %7.4f,%7.4f,%7.4f r=%7.4f nx=%7.4f nz=%7.4f nr=%7.4f ny2=%7.4f\n",
+//        in[0],in[1],in[2], r,nx,nz,nr,ny2);
+  double ny= sqrt(ny2);
+  return Point(nx,ny,nz);
+}

diff --git a/transforms.hh b/transforms.hh
new file mode 100644 (file)
index 0000000..e398d27
--- /dev/null
+++ b/transforms.hh
@@ -0,0 +1,91 @@
+/*
+ * Transformations
+ */
+
+#ifndef TRANSFORMS_HH
+#define TRANSFORMS_HH
+
+#include "library.hh"
+
+class RotationTransform: public Transform {
+protected:
+  double s,c;
+public:
+  RotationTransform(double theta);
+};
+
+class XZRotationTransform : public RotationTransform {
+  /*
+   *          Y |
+   *            |    <\ (positive theta)
+   *            |     /
+   *            /-------- X
+   *           /
+   *        Z / \--->
+   */
+public:
+  XZRotationTransform(double theta) : RotationTransform(theta) { }
+  Point operator()(Point);
+};
+
+class YZRotationTransform : public RotationTransform {
+  /*
+   *          Y | _
+   *            |/ \ (positive theta)
+   *          ^ |  V
+   *          | /-------- X
+   *          \/
+   *        Z /
+   */
+public:
+  YZRotationTransform(double theta) : RotationTransform(theta) { }
+  Point operator()(Point);
+};
+
+class TearOpenTransform : public Transform {
+  double factor;
+  double apexy;
+public:
+  TearOpenTransform(double amt, double apy=1.3);
+  Point operator()(Point in);
+};
+
+class BulgeTransform : public Transform {
+  double amount;
+  double variance;
+  double apexy;
+public:
+  BulgeTransform(double amt, double stddev, double apy=1.1);
+  Point operator()(Point in);
+};
+
+class TwistTransform : public Transform {
+  double factor;
+  double apexy;
+public:
+  TwistTransform(double amt, double apy=1.2) { factor=amt; apexy=apy; }
+  Point operator()(Point in);
+};
+
+class ShearTransform : public Transform {
+  Point epicentre;
+  Point magnitude;
+  Point deviation;
+public:
+  ShearTransform(Point ec, Point mag, Point dev) {
+    epicentre=ec; magnitude=mag; deviation=dev;
+  }
+  Point operator()(Point in);
+};
+
+class MagicTransform : public Transform {
+  double xkonstant, xfactor, yoffset, yfactor, zfactor;
+public:
+  MagicTransform(double xk, double xf, double zf,
+                 double yfudge=0.25, double yo=0.5, double yf=-0.5) {
+    xkonstant=xk; xfactor=xf; yoffset=yo+yfudge; yfactor=yf; zfactor=zf;
+  }
+  Point operator()(Point in);
+};
+
+#endif

diff --git a/x11.cc b/x11.cc
new file mode 100644 (file)
index 0000000..b724594
--- /dev/null
+++ b/x11.cc
@@ -0,0 +1,46 @@
+/*
+ * X11 functions
+ */
+
+
+#include "x11.hh"
+
+X11Output::X11Output() {
+  XGCValues gcvalues;
+  display= XOpenDisplay(0);
+  window= XCreateSimpleWindow(display,
+                              DefaultRootWindow(display),
+                              0,0, 500,500, 0,0,0);
+  gcvalues.background= 0;
+  fabric= XCreateGC(display,window,GCBackground,&gcvalues);
+  gcvalues.foreground= 1;
+  gcvalues.background= 1;
+  mesh= XCreateGC(display,window,GCForeground|GCBackground,&gcvalues);
+  XSelectInput(display,window,0);
+  XMapWindow(display,window);
+  XFlush(display);
+}
+
+void X11Output::startimage() {
+  XClearWindow(display,window);
+}
+
+void X11Output::endimage() {
+  XFlush(display);
+}
+
+X11Output::~X11Output() {
+  XCloseDisplay(display);
+}
+
+void X11Output::drawcell(const Point* list, int n) {
+  XPoint xp[n+1];
+  for (int i=0; i<n; i++) {
+    Onscreen here= Onscreen(list[i]);
+    xp[i].x= (int)((here.x+1.0)*250.0);
+    xp[i].y= (int)((-here.y+1.0)*250.0);
+  }
+  XFillPolygon(display,window,fabric,xp,n,Nonconvex,CoordModeOrigin);
+  xp[n]= xp[0];
+  XDrawLines(display,window,mesh,xp,n+1,CoordModeOrigin);
+}

diff --git a/x11.hh b/x11.hh
new file mode 100644 (file)
index 0000000..d711741
--- /dev/null
+++ b/x11.hh
@@ -0,0 +1,26 @@
+/*
+ * X11 functions
+ */
+
+#ifndef X11_HH
+#define X11_HH
+
+#include <stdlib.h>
+#include <X11/Xlib.h>
+
+#include "output.hh"
+
+struct X11Output : Output {
+  X11Output();
+  ~X11Output();
+  void drawcell(const Point*, int);
+  void startimage();
+  void endimage();
+private:
+  Display *display;
+  Window window;
+  GC fabric;
+  GC mesh;
+};
+
+#endif