X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=project.c;fp=project.c;h=0000000000000000000000000000000000000000;hp=1ce6db66d492867a869ed833b2a9e9c6931067a8;hb=512fa4d69efa3f71c61a581c137ee6da2a9bf2f9;hpb=47850b1d988483f5939b5bdb4a3bd60604f6930c

diff --git a/project.c b/project.c
deleted file mode 100644
index 1ce6db6..0000000
--- a/project.c
+++ /dev/null
@@ -1,516 +0,0 @@
-/*
- * Displays a conformation
- */
-
-#include <X11/Xlib.h>
-#include <X11/Xutil.h>
-
-#include "mgraph.h"
-
-#define MAXTRIS (N*2)
-
-typedef struct { double vertex[3][D3]; } Triangle;
-
-static Triangle trisbuffer[MAXTRIS], *displaylist[MAXTRIS];
-static int ntris;
-static Vertices conformation;
-
-static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
-static GSL_MATRIX(transform);
-
-const char *input_filename;
-
-static void read_input(void) {
-  FILE *f;
-  int r;
-  
-  f= fopen(input_filename, "rb");  if (!f) diee("input file");
-  errno= 0;
-  r= fread(&conformation,sizeof(conformation),1,f);  if (r!=1) diee("fread");
-  fclose(f);
-}
-
-static void transform_coordinates(void) {
-  double result[D3];
-  GSL_VECTOR(result);
-  gsl_vector input_gsl= { D3,1 };
-
-  int v, k;
-  
-  FOR_VERTEX(v) {
-    input_gsl.data= &conformation[v][0];
-    GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
-		       0.0, &result_gsl) );
-    K conformation[v][k]= result[k];
-  }
-}
-
-static void addtriangle(int va, int vb, int vc) {
-  Triangle *t= &trisbuffer[ntris];
-  int k;
-  
-  assert(ntris < MAXTRIS);
-  K {
-    t->vertex[0][k]= conformation[va][k];
-    t->vertex[1][k]= conformation[vb][k];
-    t->vertex[2][k]= conformation[vc][k];
-  }
-  displaylist[ntris++]= t;
-}
-
-static void generate_display_list(void) {
-  int vb, ve[3], e;
-
-  ntris= 0;
-  FOR_VERTEX(vb) {
-    /* We use the two triangles in the parallelogram vb, vb+e0, vb+e1, vb+e2.
-     * We go round each triangle clockwise (although our surface is non-
-     * orientable so it shouldn't matter).
-     */
-    for (e=0; e<3; e++) ve[e]= EDGE_END2(vb,e);
-    if (ve[1]>=0) {
-      if (ve[0]>=0) addtriangle(vb,ve[0],ve[1]);
-      if (ve[2]>=0) addtriangle(vb,ve[1],ve[2]);
-    }
-  }
-}    
-
-static int dl_compare(const void *tav, const void *tbv) {
-  int i;
-  const Triangle *const *tap= tav, *ta= *tap;
-  const Triangle *const *tbp= tbp, *tb= *tbp;
-  double za=0, zb=0;
-  for (i=0; i<3; i++) {
-    za += ta->vertex[i][2];
-    zb += tb->vertex[i][2];
-  }
-  return za > zb ? -1 :
-         za < zb ? +1 : 0;
-}
-
-static void sort_display_list(void) {
-  qsort(displaylist, ntris, sizeof(*displaylist), dl_compare);
-}
-
-/*---------- X stuff ----------*/
-
-#define WSZ 400
-
-typedef struct { GC fillgc, linegc; } DrawingMode;
-
-static Display *display;
-static Pixmap pixmap, doublebuffers[2];
-static Window window;
-
-static DrawingMode dmred, dmblue, dmwhite;
-static const DrawingMode *dmcurrent;
-static int wwidth=WSZ, wheight=WSZ, wmindim=WSZ, wmaxdim=WSZ;
-static int ncut, currentbuffer, x11depth, x11screen;
-XVisualInfo visinfo;
-
-static double sizeadj_scale= 0.3, eyes_apart, scale_wmindim;
-static double eye_z= -10, eye_x=0;
-static double cut_z= -9;
-
-static void drawtriangle(const Triangle *t) {
-  XPoint points[4];
-  int i;
-
-  for (i=0; i<3; i++) {
-    double *v= t->vertex[i];
-    double x= v[0];
-    double y= v[1];
-    double z= v[2];
-
-    if (z < cut_z) { ncut++; return; }
-    
-    double zezezp= eye_z / (eye_z - z);
-    points[i].x= scale_wmindim * (zezezp * (x - eye_x) + eye_x) + wwidth/2;
-    points[i].y= scale_wmindim * (zezezp *  y                 ) + wheight/2;
-  }
-  points[3]= points[0];
-
-  XA( XFillPolygon(display,pixmap, dmcurrent->fillgc,
-		   points,3,Convex,CoordModeOrigin) );
-  XA( XDrawLines(display,pixmap, dmcurrent->linegc,
-		 points, 4,CoordModeOrigin) );
-}
-
-static const unsigned long core_event_mask=
-  ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask;
-
-static void mkpixmaps(void) {
-  for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
-    XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,x11depth) );
-    doublebuffers[currentbuffer]= pixmap;
-  }
-  currentbuffer= 0;
-}
-
-static void mkgcs(DrawingMode *dm, unsigned long planes) {
-  XGCValues gcv;
-
-  gcv.function= GXcopy;
-  gcv.foreground= WhitePixel(display,x11screen);
-  gcv.plane_mask= planes;
-  dm->linegc= XCreateGC(display,pixmap,
-			GCFunction|GCForeground|GCPlaneMask,
-			&gcv);
-
-  gcv.function= GXclear;
-  dm->fillgc= XCreateGC(display,pixmap,
-			GCFunction|GCPlaneMask,
-			&gcv);
-}
-
-static void display_prepare(void) {
-  XSetWindowAttributes wa;
-  XSizeHints hints;
-  
-  XA( display= XOpenDisplay(0) );
-  x11screen= DefaultScreen(display);
-  x11depth= DefaultDepth(display,x11screen);
-  XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&visinfo) );
-  
-  wa.event_mask= core_event_mask;
-  XA( window= XCreateWindow(display, DefaultRootWindow(display),
-			    0,0, wwidth,wheight, 0,x11depth,
-			    InputOutput, visinfo.visual,
-			    CWEventMask, &wa) );
-
-  hints.flags= USPosition;
-  hints.x= 10;
-  hints.y= 10;
-  XSetWMNormalHints(display,window,&hints);
-
-  mkpixmaps();
-
-  mkgcs(&dmwhite, AllPlanes);
-  mkgcs(&dmblue, visinfo.blue_mask);
-  mkgcs(&dmred, visinfo.red_mask);
-}
-
-static void drawtriangles(const DrawingMode *dm) {
-  Triangle *const *t;
-  int i;
-
-  dmcurrent= dm;
-  for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
-    drawtriangle(*t);
-}
-
-static void display_conformation(void) {
-  pixmap= doublebuffers[currentbuffer];
-
-  XA( XFillRectangle(display,pixmap,dmwhite.fillgc,0,0,wwidth,wheight) );
-
-  if (eyes_apart > 0) {
-    const double preferred=0.05, beyond=0.07;
-
-    eye_x= eyes_apart < preferred ? eyes_apart :
-           eyes_apart < beyond ? preferred :
-           eyes_apart - (beyond - preferred);
-    eye_x /= sizeadj_scale;
-    drawtriangles(&dmblue);
-    eye_x= -eye_x;
-    drawtriangles(&dmred);
-  } else {
-    drawtriangles(&dmwhite);
-    printf("shown, %d/%d triangles cut\n", ncut, ntris);
-  }
-  
-  XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
-  XA( XClearWindow(display,window) );
-  currentbuffer= !currentbuffer;
-}
-
-static void show(void) {
-  scale_wmindim= sizeadj_scale * wmindim;
-  read_input();
-  transform_coordinates();
-  generate_display_list();
-  sort_display_list();
-  display_conformation();
-}
-
-typedef struct {
-  const char *name;
-  void (*start)(void);
-  void (*delta)(double dx, double dy);
-  void (*conclude)(void);
-  void (*abandon)(void);
-} Drag;
-
-#define DRAG(x)					\
-  static const Drag drag_##x= {			\
-    #x, drag_##x##_start, drag_##x##_delta,	\
-    drag_##x##_conclude, drag_##x##_abandon	\
-  }
-
-#define DRAG_SAVING(x, thing)				\
-  static typeof(thing) original_##thing;		\
-  static void drag_##x##_start(void) {			\
-    memcpy(&original_##thing, &thing, sizeof(thing));	\
-  }							\
-  static void drag_##x##_conclude(void) { }		\
-  static void drag_##x##_abandon(void) {		\
-    memcpy(&thing, &original_##thing, sizeof(thing));	\
-    show();						\
-  }							\
-  DRAG(x)
-
-static void drag_none_start(void) { }
-static void drag_none_delta(double dx, double dy) { }
-static void drag_none_conclude(void) { }
-static void drag_none_abandon(void) { }
-DRAG(none);
-
-static void pvectorcore(const char *n, double v[D3]) {
-  int k;
-  printf("%10s [ ",n);
-  K printf("%# 10.10f ",v[k]);
-  printf("]\n");
-}
-static void pvector(const char *n, double v[D3]) {
-  pvectorcore(n,v);
-  putchar('\n');
-}
-static void pmatrix(const char *n, double m[D3][D3]) {
-  int j;
-  for (j=0; j<D3; j++) { pvectorcore(n,m[j]); n=""; }
-  putchar('\n');
-}
-#define PMATRIX(x) pmatrix(#x,x);
-
-static void drag_rotate_delta(double dx, double dy) {
-  /* We multiple our transformation matrix by a matrix:
-   *
-   * If we just had y movement, we would rotate about x axis:
-   *  rotation X = [  1    0   0 ]
-   *               [  0   cy  sy ]
-   *               [  0  -sy  cy ]
-   *  where cy and sy are sin and cos of y rotation
-   *
-   * But we should pre-rotate this by a rotation about the z axis
-   * to get it to the right angle (to include x rotation).  So
-   * we make cy and sy be cos() and sin(hypot(x,y)) and use
-   * with cr,sr as cos() and sin(atan2(y,y)):
-   *
-   * Ie we would do  T' = R^T X R T   where
-   *             or  T' =    C    T   where  C = R^T X R  and
-   *
-   *  adjustment R = [  cr  sr  0 ]
-   *                 [ -sr  cr  0 ]
-   *                 [  0    0  1 ]
-   */
-
-  double rotx[D3][D3], adjr[D3][D3];
-  GSL_MATRIX(rotx);
-  GSL_MATRIX(adjr);
-
-  static double temp[D3][D3], change[D3][D3];
-  static GSL_MATRIX(temp);
-  static GSL_MATRIX(change);
-
-  double d= hypot(dx,dy);
-  if (d < 1e-6) return;
-
-  double ang= d*2.0;
-  
-  double cy= cos(ang);
-  double sy= sin(ang);
-  double cr= -dy / d;
-  double sr=  dx / d;
-  printf("\n d=%g cy,sy=%g,%g cr,sr=%g,%g\n\n", d,cy,sy,cr,sr);
-  
-  rotx[0][0]=   1;    rotx[0][1]=   0;     rotx[0][2]=  0;
-  rotx[1][0]=   0;    rotx[1][1]=  cy;     rotx[1][2]= sy;
-  rotx[2][0]=   0;    rotx[2][1]= -sy;     rotx[2][2]= cy;
-  PMATRIX(rotx);
-
-  adjr[0][0]=  cr;    adjr[0][1]=  sr;     adjr[0][2]=  0;
-  adjr[1][0]= -sr;    adjr[1][1]=  cr;     adjr[1][2]=  0;
-  adjr[2][0]=   0;    adjr[2][1]=   0;     adjr[2][2]=  1;
-  PMATRIX(adjr);
-
-  GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
-		     &rotx_gsl,&adjr_gsl,
-		     0.0, &temp_gsl) );
-  PMATRIX(temp);
-  
-  GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
-		     &adjr_gsl,&temp_gsl,
-		     0.0, &change_gsl) );
-  PMATRIX(change);
-
-  static double skew[D3][D3];
-  static GSL_MATRIX(skew);
-  
-  GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
-		     &change_gsl,&transform_gsl,
-		     0.0, &skew_gsl) );
-  PMATRIX(skew);
-
-  memcpy(&transform,&skew,sizeof(transform));
-  show();
-  return;
-
-  /* Now we want to normalise skew, the result becomes new transform */
-  double svd_v[D3][D3];
-  GSL_MATRIX(svd_v);
-
-  double sigma[D3], tau[D3];
-  GSL_VECTOR(sigma);
-  GSL_VECTOR(tau);
-  
-  /* We use notation from Wikipedia Polar_decomposition
-   *     Wikipedia's  W      is GSL's U
-   *     Wikipedia's  Sigma  is GSL's S
-   *     Wikipedia's  V      is GSL's V
-   *     Wikipedia's  U      is our desired result
-   * Wikipedia which says if the SVD is    A = W Sigma V*
-   * then the polar decomposition is       A = U P
-   *   where                               P = V Sigma V*
-   *   and                                 U = W V*
-   */
-  
-  GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
-  pmatrix("W",    skew);
-  pvector("Sigma",sigma);
-  pmatrix("V",    svd_v);
-
-  /* We only need U, not P. */
-  GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
-		     &skew_gsl,&svd_v_gsl,
-		     0.0,&transform_gsl) );
-
-  pmatrix("U", transform);
-
-  printf("drag_rotate_delta...\n");
-  show();
-}
-DRAG_SAVING(rotate, transform);
-
-static void drag_sizeadj_delta(double dx, double dy) {
-  sizeadj_scale *= pow(3.0, -dy);
-  show();
-}
-DRAG_SAVING(sizeadj, sizeadj_scale);
-
-static void drag_3d_delta(double dx, double dy) {
-  const double min_eyes_apart= -0.02;
-  eyes_apart += dx * 0.1;
-  if (eyes_apart < min_eyes_apart) eyes_apart= min_eyes_apart;
-  printf("sizeadj eyes_apart %g\n", eyes_apart);
-  show();
-}
-DRAG_SAVING(3d, eyes_apart);
-
-static const Drag *drag= &drag_none;
-
-static int drag_last_x, drag_last_y;
-
-static void drag_position(int x, int y) {
-  drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
-	      (y - drag_last_y) * 1.0 / wmaxdim);
-  drag_last_x= x;
-  drag_last_y= y;
-}
-
-static void event_button(XButtonEvent *e) {
-  if (e->window != window || !e->same_screen) return;
-  if (e->type == ButtonPress) {
-    if (e->state || drag != &drag_none) {
-      printf("drag=%s press state=0x%lx abandon\n",
-	     drag->name, (unsigned long)e->state);
-      drag->abandon();
-      drag= &drag_none;
-      return;
-    }
-    switch (e->button) {
-    case Button1: drag= &drag_rotate;  break;
-    case Button2: drag= &drag_sizeadj; break;
-    case Button3: drag= &drag_3d;      break;
-    default: printf("unknown drag start %d\n", e->button);
-    }
-    printf("drag=%s press button=%lu start %d,%d\n",
-	   drag->name, (unsigned long)e->button, e->x, e->y);
-    drag_last_x= e->x;
-    drag_last_y= e->y;
-    drag->start();
-  }
-  if (e->type == ButtonRelease) {
-    printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
-    drag_position(e->x, e->y);
-    drag->conclude();
-    drag= &drag_none;
-  }
-}
-
-static void event_motion(int x, int y) {
-  printf("drag=%s motion %d,%d\n", drag->name, x, y);
-  drag_position(x,y);
-}
-
-static void event_config(XConfigureEvent *e) {
-  if (e->width == wwidth && e->height == wheight)
-    return;
-
-  wwidth= e->width;  wheight= e->height;
-  wmaxdim= wwidth > wheight ? wwidth : wheight;
-  wmindim= wwidth < wheight ? wwidth : wheight;
-  
-  XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
-  for (currentbuffer=0; currentbuffer<2; currentbuffer++)
-    XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
-
-  mkpixmaps();
-  show();
-}
-
-int main(int argc, const char *const *argv) {
-  XEvent event;
-  int k;
-  int motion_deferred=0, motion_x=-1, motion_y=-1;
-  
-  if (argc != 2 || argv[1][0]=='-') {
-    fputs("need filename\n",stderr); exit(8);
-  }
-  input_filename= argv[1];
-
-  read_input();
-  K transform[k][k]= 1.0;
-  display_prepare();
-  show();
-
-  XMapWindow(display,window);
-  for (;;) {
-    if (motion_deferred) {
-      int r= XCheckMaskEvent(display,~0UL,&event);
-      if (!r) {
-	event_motion(motion_x, motion_y);
-	motion_deferred=0;
-	continue;
-      }
-    } else {
-      XNextEvent(display,&event);
-    }
-    switch (event.type) {
-
-    case ButtonPress:
-    case ButtonRelease:     event_button(&event.xbutton);     break;
-      
-    case ConfigureNotify:   event_config(&event.xconfigure);  break;
-
-    case MotionNotify:
-      motion_x= event.xmotion.x;
-      motion_y= event.xmotion.y;
-      motion_deferred= 1;
-      continue;
-      
-    default:
-      printf("unknown event type %u 0x%x\n", event.type,event.type);
-    }
-  }
-}
-