--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+
+OBJS= main.o moebius.o output.o library.o
+
+a.out: $(OBJS)
+ g++ -o a.out $(OBJS)
+
--- /dev/null
+/*
+ * 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);
+}
--- /dev/null
+/*
+ * 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
--- /dev/null
+epsffit -c -s 50 450 550 750 t
--- /dev/null
+/*
+ * 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();
+}
--- /dev/null
+/*
+ * 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
--- /dev/null
+tear=1
+twist=-1
+theta=1.57
+bulge=1.5
+sheary=3
+un=20
+tn=35
+x11
--- /dev/null
+/*
+ * 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));
+}
--- /dev/null
+/*
+ * 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
--- /dev/null
+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 ./
--- /dev/null
+/*
+ * 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";
+ }
+ }
+ }
+}
+
--- /dev/null
+/*
+ * 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));
+ }
+}
--- /dev/null
+/*
+ * 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) );
+ }
+}
--- /dev/null
+/*
+ * 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
--- /dev/null
+/*
+ * 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
--- /dev/null
+/*
+ * 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;
--- /dev/null
+/*
+ * 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
--- /dev/null
+/*
+ * 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";
+}
--- /dev/null
+/*
+ * 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();
+};
--- /dev/null
+/*
+ * 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);
+}
--- /dev/null
+/*
+ * 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
--- /dev/null
+/*
+ * 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);
+}
--- /dev/null
+/*
+ * 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