2 * Displays a conformation
12 typedef struct { double vertex[3][D3]; } Triangle;
14 static Triangle trisbuffer[MAXTRIS], *displaylist[MAXTRIS];
16 static Vertices conformation;
18 static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
19 static GSL_MATRIX(transform);
21 const char *input_filename;
23 static void read_input(void) {
27 f= fopen(input_filename, "rb"); if (!f) diee("input file");
29 r= fread(&conformation,sizeof(conformation),1,f); if (r!=1) diee("fread");
33 static void transform_coordinates(void) {
36 gsl_vector input_gsl= { D3,1 };
41 input_gsl.data= &conformation[v][0];
42 GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
44 K conformation[v][k]= result[k];
48 static void addtriangle(int va, int vb, int vc) {
49 Triangle *t= &trisbuffer[ntris];
52 assert(ntris < MAXTRIS);
54 t->vertex[0][k]= conformation[va][k];
55 t->vertex[1][k]= conformation[vb][k];
56 t->vertex[2][k]= conformation[vc][k];
58 displaylist[ntris++]= t;
61 static void generate_display_list(void) {
66 /* We use the two triangles in the parallelogram vb, vb+e1, vb+e0, vb+e2.
67 * We go round each triangle clockwise (although our surface is non-
68 * orientable so it shouldn't matter).
70 for (e=0; e<3; e++) ve[e]= EDGE_END2(vb,e);
72 if (ve[1]>=0) addtriangle(vb,ve[0],ve[1]);
73 if (ve[2]>=0) addtriangle(vb,ve[2],ve[0]);
78 static int dl_compare(const void *tav, const void *tbv) {
79 const Triangle *const *tap= tav, *ta= *tap;
80 const Triangle *const *tbp= tbp, *tb= *tbp;
81 double za= ta->vertex[0][2];
82 double zb= tb->vertex[0][2];
87 static void sort_display_list(void) {
88 qsort(displaylist, ntris, sizeof(*displaylist), dl_compare);
91 /*---------- X stuff ----------*/
95 static Display *display;
96 static Pixmap pixmap, doublebuffers[2];
98 static GC linegc, fillgc;
99 static int wwidth=WSZ, wheight=WSZ, wmindim=WSZ, wmaxdim=WSZ;
100 static int ncut, currentbuffer, x11depth, x11screen;
102 static double scale= 0.3;
103 static double eye_z= -10, eye_x= 0;
104 static double cut_z= -9;
106 static void drawtriangle(const Triangle *t) {
110 for (i=0; i<3; i++) {
111 double *v= t->vertex[i];
116 if (z < cut_z) { ncut++; return; }
118 double zezezp= eye_z / (eye_z - z);
119 points[i].x= scale * wmindim * (zezezp * (x - eye_x) + eye_x) + wwidth/2;
120 points[i].y= scale * wmindim * (zezezp * y ) + wheight/2;
122 points[3]= points[0];
124 XA( XFillPolygon(display,pixmap,fillgc, points,3,Convex,CoordModeOrigin) );
125 XA( XDrawLines(display,pixmap,linegc,points, 4,CoordModeOrigin) );
128 static const unsigned long core_event_mask=
129 ButtonPressMask|ButtonReleaseMask|StructureNotifyMask|ButtonMotionMask;
131 static void mkpixmaps(void) {
132 for (currentbuffer=0; currentbuffer<2; currentbuffer++) {
133 XA( pixmap= XCreatePixmap(display,window,wwidth,wheight,x11depth) );
134 doublebuffers[currentbuffer]= pixmap;
139 static void display_prepare(void) {
141 XSetWindowAttributes wa;
145 XA( display= XOpenDisplay(0) );
146 x11screen= DefaultScreen(display);
147 x11depth= DefaultDepth(display,x11screen);
148 XA( XMatchVisualInfo(display,x11screen,x11depth, TrueColor,&vinfo) );
150 wa.event_mask= core_event_mask;
151 XA( window= XCreateWindow(display, DefaultRootWindow(display),
152 0,0, wwidth,wheight, 0,x11depth,
153 InputOutput, vinfo.visual,
156 hints.flags= USPosition;
159 XSetWMNormalHints(display,window,&hints);
163 gcv.function= GXcopy;
164 gcv.plane_mask= ~0UL;
165 gcv.foreground= WhitePixel(display,x11screen);
166 linegc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask|GCForeground, &gcv);
168 gcv.function= GXclear;
169 gcv.plane_mask= ~0UL;
170 fillgc= XCreateGC(display,pixmap, GCFunction|GCPlaneMask, &gcv);
173 static void display_conformation(void) {
177 pixmap= doublebuffers[currentbuffer];
178 XA( XFillRectangle(display,pixmap,fillgc,0,0,wwidth,wheight) );
179 for (i=0, t=displaylist, ncut=0; i<ntris; i++, t++)
181 printf("shown, %d/%d triangles cut\n", ncut, ntris);
183 XA( XSetWindowBackgroundPixmap(display,window,pixmap) );
184 XA( XClearWindow(display,window) );
185 currentbuffer= !currentbuffer;
188 static void show(void) {
190 transform_coordinates();
191 generate_display_list();
193 display_conformation();
199 void (*delta)(double dx, double dy);
200 void (*conclude)(void);
201 void (*abandon)(void);
205 static const Drag drag_##x= { \
206 #x, drag_##x##_start, drag_##x##_delta, \
207 drag_##x##_conclude, drag_##x##_abandon \
210 #define DRAG_SAVING(x, thing) \
211 static typeof(thing) original_##thing; \
212 static void drag_##x##_start(void) { \
213 memcpy(&original_##thing, &thing, sizeof(thing)); \
215 static void drag_##x##_conclude(void) { } \
216 static void drag_##x##_abandon(void) { \
217 memcpy(&thing, &original_##thing, sizeof(thing)); \
222 static void drag_none_start(void) { }
223 static void drag_none_delta(double dx, double dy) { }
224 static void drag_none_conclude(void) { }
225 static void drag_none_abandon(void) { }
228 static void pmatrix(const char *n, double m[D3][D3]) {
230 for (j=0; j<D3; j++) {
232 K printf("%# f ",m[j][k]);
238 #define PMATRIX(x) pmatrix(#x,x);
240 static void drag_rotate_delta(double dx, double dy) {
241 /* We multiple our transformation matrix by this matrix:
245 * and then renormalise.
248 static double rotateby[D3][D3]= {{1,0,0},{0,1,0},{-20,-20,1}};
249 static GSL_MATRIX(rotateby);
252 gsl_matrix qr_gsl= { D3,D3,D3,&qr[0][0] };
262 GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
263 &rotateby_gsl,&transform_gsl, 0.0,&qr_gsl) );
265 pmatrix("input", qr);
267 GA( gsl_linalg_QR_decomp(&qr_gsl, &tau_gsl) );
269 pmatrix("mangled", qr);
271 GA( gsl_linalg_QR_unpack(&qr_gsl, &tau_gsl,
272 &transform_gsl, &rotateby_gsl /*dummy*/) );
274 pmatrix("Q", transform);
275 pmatrix("R", rotateby);
277 printf("drag_rotate_delta...\n");
280 DRAG_SAVING(rotate, transform);
282 static void drag_scale_delta(double dx, double dy) {
283 scale *= pow(3.0, dy);
286 DRAG_SAVING(scale, scale);
288 static const Drag *drag= &drag_none;
290 static int drag_last_x, drag_last_y;
292 static void drag_position(int x, int y) {
293 drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
294 (y - drag_last_y) * 1.0 / wmaxdim);
299 static void event_button(XButtonEvent *e) {
300 if (e->window != window || !e->same_screen) return;
301 if (e->type == ButtonPress) {
302 if (e->state || drag != &drag_none) {
303 printf("drag=%s press state=0x%lx abandon\n",
304 drag->name, (unsigned long)e->state);
310 case Button1: drag= &drag_rotate; break;
311 case Button2: drag= &drag_scale; break;
312 default: printf("unknown drag start %d\n", e->button);
314 printf("drag=%s press button=%lu start %d,%d\n",
315 drag->name, (unsigned long)e->button, e->x, e->y);
320 if (e->type == ButtonRelease) {
321 printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
322 drag_position(e->x, e->y);
328 static void event_motion(int x, int y) {
329 printf("drag=%s motion %d,%d\n", drag->name, x, y);
333 static void event_config(XConfigureEvent *e) {
334 if (e->width == wwidth && e->height == wheight)
337 wwidth= e->width; wheight= e->height;
338 wmaxdim= wwidth > wheight ? wwidth : wheight;
339 wmindim= wwidth < wheight ? wwidth : wheight;
341 XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
342 for (currentbuffer=0; currentbuffer<2; currentbuffer++)
343 XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
349 int main(int argc, const char *const *argv) {
352 int motion_deferred=0, motion_x=-1, motion_y=-1;
354 if (argc != 2 || argv[1][0]=='-') {
355 fputs("need filename\n",stderr); exit(8);
357 input_filename= argv[1];
360 K transform[k][k]= 1.0;
364 XMapWindow(display,window);
366 if (motion_deferred) {
367 int r= XCheckMaskEvent(display,~0UL,&event);
369 event_motion(motion_x, motion_y);
374 XNextEvent(display,&event);
376 switch (event.type) {
379 case ButtonRelease: event_button(&event.xbutton); break;
381 case ConfigureNotify: event_config(&event.xconfigure); break;
385 motion_x= event.xmotion.x;
386 motion_y= event.xmotion.y;
392 printf("unknown event type %u 0x%x\n", event.type,event.type);