-
Notifications
You must be signed in to change notification settings - Fork 1
/
idkit.c
1954 lines (1786 loc) · 64.8 KB
/
idkit.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*
* idkit.c: fit Identakit N-body model to interacting galaxy XYV data.
* Copyright (C) 2012 by Joshua Edward Barnes, Honolulu, Hawaii. This
* is free software, distributed under the terms of GNU General Public
* License. See <http://www.gnu.org/licenses/> for details.
*/
#include "stdinc.h"
#include "mathfns.h"
#include "getparam.h"
#include "vectmath.h"
#include "filestruct.h"
#include "phatbody.h"
#if defined(MACOSX)
#include <GLUT/glut.h>
#else
#include <GL/glut.h>
#endif
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
// ______________________________________________________
// Command line parameters and their description strings.
string defv[] = { ";Fit Identakit model to XYVz data.",
";Drag mouse to adjust view; command keys:",
";? <sp> <del> n p x y z r R l m 1 2 c s o v",
";t Q q D L j J V M w W d C Z S b f $ <esc>",
"jhat=???", ";Spin direction and weight data.",
";Used to select disk particles to display.",
"sims=???", ";Patterns for identakit simulation files.",
";Embedded printf directive uses current",
";value of step to generate filenames.",
"zrot=0.0", ";Rotate simulation about initial z axis",
"isim=0", ";Initially display this simulation",
"step=0x200", ";Initially display this step number",
"incr=0x020", ";Increment/decrement for step number",
"ncenter=2048", ";Number of bodies used to track centers",
"maxview=32768", ";Maximum number of bodies to display",
"panelsize=400", ";Size of panels; entire window is 2x2",
"xyimage=", ";Image for XY projection",
"vyimage=", ";Image for VY projection",
"xvimage=", ";Image for XV projection",
"datacube=", ";XYV data in particle form",
"initview=", ";Initial view file",
"saveview=idk%03d.dat", ";Output pattern for view files",
"saveimage=idk%03d.ppm", ";Output pattern for screen dump",
"spindepth=4", ";Tessellation depth for spin ball",
"viewdepth=2", ";Tessellation depth for view search",
"xwidth=250.0", ";Image X,Y width in kpc",
"vwidth=500.0", ";Image V width in km/sec",
"script=", ";String of keyboard commands",
"VERSION=2.2mod" VARIANT, ";Joshua Barnes 11 June 2012.",
";(modified George Privon, 05 October 2016).",
";This program is free software; you can",
";redistribute it under certain conditions.",
";This program comes with NO WARRANTY.",
NULL,
};
// __________________________________________________________
// idkstate: structure for current and past identikit states.
#define NUMPARS 10 // num. of parameter-pair variables
typedef struct {
real par[NUMPARS][2]; // values of variable parameters
matrix pmat; // preliminary transformation mat.
int step; // current value of step number
int isim; // index into list of simulations
} idkstate;
// _____________________________________________________
// XYROT, ..., JVIEW: indexes into parameter-pair array.
#define XYROT 0 // rotation about X and Y axes
#define SCALE 1 // X,Y scale and rotation about Z
#define DISK1 2 // incl. and peri. arg. for disk 1
#define DISK2 3 // incl. and peri. arg. for disk 2
#define XYOFF 4 // offsets in X and Y directions
#define VELOS 5 // V offset and scale factor
#define ITRAN 6 // image contrast and zeropoint
#define BOXXY 7 // selection box X and Y coords.
#define BOXVR 8 // selection box V coord and size
#define JVIEW 9 // rotate view of spin ball
// ________________________________________________
// box: structure for phase-space selection region.
typedef struct {
vector XYbox; // screen X, Y coords. of box
real Vbox; // screen V coord. of box
real Rbox; // size of box
bool velflag; // select V as well as X and Y
int disk; // which disk box selects
} box;
#define MAXBOX 32 // max. number of phase-space boxes
// ___________________________________________________________________
// facet: structure for triangular face of solid approximating a ball.
typedef struct {
vector vert1, vert2, vert3; // coordinates of corners
vector midp; // unit vector through midpoint
real spin[2], prod[2], view[2]; // values on ball, indexed by disk
} facet;
#define NUMICOS 20 // num. of faces on icosahedron
// ______________________________________________________________________
// XZSHOW, ..., VWSHOW: content of lower-right panel, selected by lrshow.
#define XZSHOW 0 // show (X,Z) projection
#define BXSHOW (MAXBOX) // show spin ball (1:BXSHOW)
#define PRSHOW (MAXBOX+1) // show product ball
#define VWSHOW (MAXBOX+2) // show view ball
// ______________________________________
// Global variables and input parameters.
idkstate idk; // current identikit state
bool parset[NUMPARS]; // flag changed parameter-pairs
#define MAXHIST 1024 // max. number of saved states
idkstate old[MAXHIST]; // past identikit history
int ihist = 0; // index into history array
int nhist = 0; // num. of items in history array
box boxtab[MAXBOX+1]; // storage for phase-space boxes
int nbox = 0; // box indices run from 1 to nbox
facet *ball; // tessellation of icosahedron
int nfacet; // number of facets in tessellation
bodyptr btab = NULL; // identikit array of bodies
int nbody = 0; // number of bodies in btab[]
real tnow; // current time from input file
bodyptr vtab = NULL; // array of visible bodies
int maxview, nview1, nview2; // max size, number from each disk
bodyptr ctab = NULL; // data cube body array
int ncube = 0; // number of bodies in data cube
string *sims; // patterns for input filenames
string *jhat; // names of jhat files for above
string *zrot; // values for initial z rot'n
int nsims; // number of names in these lists
GLsizei pansize, imgsize = -1; // size of panels and bkgnd images
GLubyte *xyimage = NULL, *xyimage1; // image behind XY projection
GLubyte *vyimage = NULL, *vyimage1; // image behind VY projection
GLubyte *xvimage = NULL, *xvimage1; // image behind XV projection
int butbind[] = { XYROT,DISK1,DISK2 }; // map mouse buttons to parameters
int actpars = -1, actlast = XYROT; // current and last param adjusted
int xlast, ylast; // previous position of mouse
int ncenter; // bodies used to find centers
real disktol = 0.01; // sets number of bodies visible
vector cpos1, cpos2; // screen positions of centers
real cvel1[2], cvel2[2]; // screen vel. ranges of centers
bool cposflag = FALSE; // TRUE if cpos vectors are set
bool cvelflag = FALSE; // TRUE if cvel ranges are set
int centflag = 0; // if > 0, lock centers at cpos
int lrshow = XZSHOW; // what to show in LR panel
int parflag = 1; // display parameter values
bool cubeflag = FALSE; // if TRUE, show data cube
int zoomflag = 0; // if >0, zoom on one panel
int ndirect; // number of directions to scan
char *script; // command keys to be executed
bool animate = FALSE; // TRUE if animation in progress
char msgbuf[256], titlebuf[256]; // buffers for messages and title
// __________________________________________________________
// Function and procedure prototypes, in order of appearance.
void display(void);
void seldisks(void);
void lockcent(void);
void mapdisks(void);
void adjustimage(GLubyte *img1, GLubyte *img0, real scale, real bias);
void displayimage(GLubyte *image);
void displaycube(int xind, int yind, real xfield, real vfield);
void displaydisks(int xind, int yind);
void displayball(void);
bool loadball(matrix vmat, real xbox, real ybox, real vbox,
real rbox, real hbox, int disk);
facet *showball(matrix tmat, real *rgb1, real *rgb2, int ind);
void showspins(matrix jmat, vector spin1, vector spin2);
void displaypars(void);
void keyboard(unsigned char key, int x, int y);
void printhelp(void);
void storehist(void);
void scanboxes(void);
void scanviews(int arg);
void findmax(void);
void printscale(void);
void mousebut(int but, int state, int x, int y);
void mousemove(int x, int y);
void special(int key, int x, int y);
void getjhat(string jhat, string zrot, bool firstcall);
void getdata(string spat, string zrot, bool firstcall);
void getcube(string cube);
void getview(void);
void putview(void);
void getimage(GLubyte **image, GLubyte **imcpy, string ifile);
string inputline(stream istr);
void putscreen(void);
void clrpars(bool all);
void getpars(matrix vmat, real *xfield, vector spin1, vector spin2,
vector xyoff, real *vfield, real *voff, real *scale, real *bias,
vector XYbox, real *Vbox, real *Rbox, matrix jmat);
void setpars(real *xyzrot, real *xfield, vector spin1, vector spin2,
vector xyoff, vector XYbox, real *Vbox, real *Rbox);
void showtext(string str, int x, int y);
void showcross(real x, real y, real size);
void showbox(real x, real y, real size);
void rotmatrix(matrix rmat, real xrot, real yrot, real zrot);
string *padlist(string *old, int len, string name);
void findcenters(void);
bool centbbox(real xmin[2], real xmax[2], real ymin[2], real ymax[2],
GLubyte *img, int size);
void buildsphere(int lev);
void buildtriangle(vector v1, vector v2, vector v3, int lev);
// __________________
// Macro definitions.
#define rsinD(x) (((x) == 0 || (x) == 180) ? 0 : rsin((x) * (PI / 180)))
#define rcosD(x) (((x) == 90 || (x) == 270) ? 0 : rcos((x) * (PI / 180)))
#define rasinD(x) (rasin(x) * (180 / PI))
#define racosD(x) (racos(x) * (180 / PI))
#define ratan2D(x,y) ((x<0 ? 360 : 0) + ratan2(x,y) * (180 / PI))
#define PERDIST 4.0 // viewing distance for ball
#define PerTrans(x,z) (1.125 * (x) / (1.0 - (z) / PERDIST))
// _________________________________________________________________
// main: read initial data, initialize display, and start main loop.
int main(int argc, string argv[])
{
string bodytags[] = { PosTag, VelTag, AuxTag, AuxVecTag, NULL, };
glutInit(&argc, argv);
initparam(argv, defv);
layout_body(bodytags, Precision, NDIM);
sims = burststring(getparam("sims"), ",; ");
nsims = xstrlen(sims, sizeof(string)) - 1;
if (nsims < 1)
error("%s: input file list list empty\n", getargv0());
jhat = padlist(burststring(getparam("jhat"), ",; "), nsims, "jhat");
zrot = padlist(burststring(getparam("zrot"), ",; "), nsims, "zrot");
idk.isim = getiparam("isim");
if (idk.isim < 0 || idk.isim >= nsims)
error("%s: isim out of range\n", getargv0());
idk.step = getiparam("step");
getjhat(jhat[idk.isim], zrot[idk.isim], TRUE);
getdata(sims[idk.isim], zrot[idk.isim], TRUE);
ncenter = getiparam("ncenter");
maxview = getiparam("maxview");
vtab = (bodyptr) allocate(maxview * SizeofBody);
pansize = getiparam("panelsize");
if (! strnull(getparam("xyimage")))
getimage(&xyimage, &xyimage1, getparam("xyimage"));
if (! strnull(getparam("vyimage")))
getimage(&vyimage, &vyimage1, getparam("vyimage"));
if (! strnull(getparam("xvimage")))
getimage(&xvimage, &xvimage1, getparam("xvimage"));
findcenters();
if (! strnull(getparam("datacube")))
getcube(getparam("datacube"));
clrpars(TRUE);
SETMI(idk.pmat);
buildsphere(getiparam("spindepth"));
if (! strnull(getparam("initview")))
getview();
storehist();
ndirect = MIN(NUMICOS * (1 << (2 * getiparam("viewdepth"))), nfacet);
script = getparam("script");
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
glutInitWindowSize(2 * pansize, 2 * pansize);
glutCreateWindow(getargv0());
glClearColor(0.0, 0.0, 0.0, 0.0);
glutDisplayFunc(display);
glutMouseFunc(mousebut);
glutMotionFunc(mousemove);
glutKeyboardFunc(keyboard);
glutSpecialFunc(special);
#if defined(BOLD_GRAPHICS)
glLineWidth((GLfloat) BOLD_GRAPHICS);
glPointSize((GLfloat) BOLD_GRAPHICS);
#endif
glutMainLoop();
return (0);
}
// ________________________________________________
// display: main routine called to update graphics.
void display(void)
{
int box, i;
real xfield, vfield, scale, bias;
if (parset[DISK1] || parset[DISK2])
seldisks();
if (centflag > 0)
lockcent();
box = (lrshow <= BXSHOW ? lrshow : 0); // box == 0 has no effect
getpars(NULL, &xfield, NULL, NULL, NULL, &vfield, NULL, &scale, &bias,
boxtab[box].XYbox, &boxtab[box].Vbox, &boxtab[box].Rbox, NULL);
mapdisks();
if (parset[ITRAN] && xyimage != NULL)
adjustimage(xyimage1, xyimage, scale, bias);
if (parset[ITRAN] && xvimage != NULL)
adjustimage(xvimage1, xvimage, scale, bias);
if (parset[ITRAN] && vyimage != NULL)
adjustimage(vyimage1, vyimage, scale, bias);
glClear(GL_COLOR_BUFFER_BIT);
glViewport(0, 0, 2*pansize, 2*pansize);
if (zoomflag == 0 || zoomflag == 1) {
if (zoomflag == 0)
glViewport(0, pansize, pansize, pansize);
displayimage(xyimage1);
displaycube(0, 1, xfield, vfield);
displaydisks(0, 1);
}
if (zoomflag == 0 || zoomflag == 2) {
if (zoomflag == 0)
glViewport(0, 0, pansize, pansize);
displayimage(xvimage1);
displaycube(0, 3, xfield, vfield);
displaydisks(0, 3);
}
if (zoomflag == 0 || zoomflag == 3) {
if (zoomflag == 0)
glViewport(pansize, pansize, pansize, pansize);
displayimage(vyimage1);
displaycube(3, 1, xfield, vfield);
displaydisks(3, 1);
}
if (zoomflag == 0 || zoomflag == 4) {
if (zoomflag == 0)
glViewport(pansize, 0, pansize, pansize);
if (lrshow == XZSHOW)
displaydisks(0, 2); // show XZ projection
else
displayball(); // show function on ball
}
displaypars(); // show pars on last panel
glFlush();
glutSwapBuffers();
glutSetWindowTitle(titlebuf);
for (i = 0; i < NUMPARS; i++)
parset[i] = FALSE;
if (*script != (long) NULL && !animate) // if script is playing
keyboard(*script++, 0, 0); // handle keys one-by-one
}
// _______________________________________________________________
// seldisks: select bodies in current disks, and locate centroids.
void seldisks(void)
{
vector spin1, spin2, cmpos, cmvel;
matrix zmat;
int block, nc;
bodyptr vp, bp, cp;
getpars(NULL, NULL, spin1, spin2, NULL, NULL, NULL,
NULL, NULL, NULL, NULL, NULL, NULL);
vp = NthBody(vtab, 2);
nview1 = nview2 = 0;
block = 0; // step thru body blocks
for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp))
if (AuxVec(bp)[0] == 0 && AuxVec(bp)[1] == 0 && AuxVec(bp)[2] == 0) {
#ifndef BULGELESS
if (block % 2 == 0) { // start of sphr block?
#else
if (block % 2 == 1) { // start of sphr block?
#endif
block++;
#ifndef BULGELESS
cp = NthBody(vtab, block == 1 ? 0 : 1);
#else
cp = NthBody(vtab, block == 2 ? 0 : 1);
#endif
CLRV(cmpos);
CLRV(cmvel);
nc = 0;
}
if (nc < ncenter) {
nc++;
ADDV(cmpos, cmpos, Pos(bp));
ADDV(cmvel, cmvel, Vel(bp));
DIVVS(Pos(cp), cmpos, nc);
DIVVS(Vel(cp), cmvel, nc);
}
} else {
#ifndef BULGELESS
if (block % 2 == 1) // start of disk block?
#else
if (block % 2 == 0) // start of disk block?
#endif
block++;
#ifndef BULGELESS
if (dotvp(AuxVec(bp), block==2 ? spin1 : spin2) > 1 - disktol*Aux(bp)) {
#else
if (dotvp(AuxVec(bp), block==1 ? spin1 : spin2) > 1 - disktol*Aux(bp)) {
#endif
if (vp >= NthBody(vtab, maxview)) {
sprintf(msgbuf, "#rtoo many bodies to view");
break;
}
SETV(Pos(vp), Pos(bp));
SETV(Vel(vp), Vel(bp));
vp = NextBody(vp);
#ifndef BULGELESS
(* (block==2 ? &nview1 : &nview2))++;
#else
(* (block==1 ? &nview1 : &nview2))++;
#endif
}
}
}
// __________________________________________________________________
// lockcent: adjust scale, z rotation, offsets to keep centers fixed.
void lockcent(void)
{
real dxref, dyref, xfield, dxmod, dymod, sfact, dxoff, dyoff;
matrix vmat;
vector xyoff;
bodyptr bp, c1 = NthBody(vtab, 0), c2 = NthBody(vtab, 1);
dxref = cpos2[0] - cpos1[0]; // work in screen coords.
dyref = cpos2[1] - cpos1[1];
getpars(vmat, &xfield, NULL, NULL, xyoff, NULL, NULL,
NULL, NULL, NULL, NULL, NULL, NULL);
for (bp = vtab; bp < NthBody(vtab, 2); bp = NextBody(bp)) {
MULMV(AuxVec(bp), vmat, Pos(bp));
ADDV(AuxVec(bp), AuxVec(bp), xyoff);
DIVVS(AuxVec(bp), AuxVec(bp), xfield);
}
dxmod = AuxVec(c2)[0] - AuxVec(c1)[0];
dymod = AuxVec(c2)[1] - AuxVec(c1)[1];
sfact = rsqrt((dxmod*dxmod + dymod*dymod) / (dxref*dxref + dyref*dyref));
/********** convert following to call to setpars() **********/
idk.par[SCALE][0] += (2/PI) * (ratan2(dxref, dyref) - ratan2(dxmod, dymod));
idk.par[SCALE][1] += rlog10(sfact);
parset[SCALE] = TRUE;
/********** convert previous to call to setpars() **********/
getpars(vmat, &xfield, NULL, NULL, xyoff, NULL, NULL,
NULL, NULL, NULL, NULL, NULL, NULL);
for (bp = vtab; bp < NthBody(vtab, 2); bp = NextBody(bp)) {
MULMV(AuxVec(bp), vmat, Pos(bp));
ADDV(AuxVec(bp), AuxVec(bp), xyoff);
DIVVS(AuxVec(bp), AuxVec(bp), xfield);
}
dxoff = (cpos1[0] + cpos2[0])/2 - (AuxVec(c1)[0] + AuxVec(c2)[0])/2;
dyoff = (cpos1[1] + cpos2[1])/2 - (AuxVec(c1)[1] + AuxVec(c2)[1])/2;
xyoff[0] += xfield * dxoff;
xyoff[1] += xfield * dyoff;
setpars(NULL, NULL, NULL, NULL, xyoff, NULL, NULL, NULL);
}
// _________________________________________________________________
// mapdisks: apply viewing transformation to sampled disk particles.
void mapdisks(void)
{
matrix vmat;
real xfield, vfield, voff;
vector xyoff, tmpv;
bodyptr bp;
getpars(vmat, &xfield, NULL, NULL, xyoff, &vfield, &voff,
NULL, NULL, NULL, NULL, NULL, NULL);
for (bp = vtab; bp < NthBody(vtab, 2+nview1+nview2); bp = NextBody(bp)) {
MULMV(AuxVec(bp), vmat, Pos(bp));
ADDV(AuxVec(bp), AuxVec(bp), xyoff);
DIVVS(AuxVec(bp), AuxVec(bp), xfield);
MULMV(tmpv, vmat, Vel(bp));
Aux(bp) = (tmpv[2] + voff) / vfield;
}
}
// ____________________________________________________
// adjustimage: perform linear transformation on image.
void adjustimage(GLubyte *img1, GLubyte *img0, real scale, real bias)
{
int i, imgval;
for (i = 0; i < 3 * imgsize * imgsize; i++) {
imgval = 256.0 * (bias + scale * img0[i] / 256.0);
img1[i] = MIN(255, MAX(0, imgval));
}
}
// ________________________________________
// displayimage: fill panel with RGB image.
void displayimage(GLubyte *image)
{
int effsize = (zoomflag > 0 ? 2 * pansize : pansize);
if (image != NULL) {
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(0, effsize, effsize, 0, -1.0, 1.0);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glPushMatrix();
glRasterPos2f((GLfloat) 0, (GLfloat) 0);
glPixelZoom((effsize / (real) imgsize), -(effsize / (real) imgsize));
glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
glDrawPixels(imgsize, imgsize, GL_RGB, GL_UNSIGNED_BYTE, image);
glPopMatrix();
}
}
// _________________________________________________________
// displaycube: display data cube, highlighting current box.
void displaycube(int xind, int yind, real xfield, real vfield)
{
bodyptr cp;
real XYV[3];
if (cubeflag && ctab != NULL) {
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(-0.5, 0.5, -0.5, 0.5, -100.0, 100.0);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glPushMatrix();
glBegin(GL_POINTS);
for (cp = ctab; cp < NthBody(ctab, ncube); cp = NextBody(cp)) {
XYV[0] = Pos(cp)[0] / xfield;
XYV[1] = Pos(cp)[1] / xfield;
XYV[2] = Pos(cp)[2] / vfield;
if ((1 <= lrshow && lrshow <= BXSHOW) && // lrshow is box index
ABS(XYV[0] - boxtab[lrshow].XYbox[0]) < boxtab[lrshow].Rbox &&
ABS(XYV[1] - boxtab[lrshow].XYbox[1]) < boxtab[lrshow].Rbox &&
(ABS(XYV[2] - boxtab[lrshow].Vbox) < boxtab[lrshow].Rbox ||
! boxtab[lrshow].velflag))
glColor3f(0.0, 1.0, 0.0);
else
glColor3f(0.0, 0.5, 0.0);
glVertex2f(xind<2 ? XYV[xind] : XYV[2],
yind<2 ? XYV[yind] : XYV[2]);
}
glEnd(/* GL_POINTS */);
glPopMatrix();
}
}
// _____________________________________________________________
// displaydisks: display specified projection of disk particles.
#define SelCoord(bp,ind) ((ind) < 3 ? AuxVec(bp)[ind] : Aux(bp))
void displaydisks(int xind, int yind)
{
int j, i;
bodyptr bp;
vector XYbox;
real Rbox, Vbox, Xbox, Ybox, boxlev;
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(-0.5, 0.5, -0.5, 0.5, -100.0, 100.0);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glPushMatrix();
glBegin(GL_POINTS);
for (j = 0; j < 2 * MAX(nview1, nview2); j++)
if ((j & 1) == 0 && j/2 < nview1) {
glColor3f(0.0, 1.0, 1.0);
i = j/2 + 2;
glVertex2f(SelCoord(NthBody(vtab, i), xind),
SelCoord(NthBody(vtab, i), yind));
} else if (j/2 < nview2) {
glColor3f(1.0, 0.0, 1.0);
i = j/2 + 2 + nview1;
glVertex2f(SelCoord(NthBody(vtab, i), xind),
SelCoord(NthBody(vtab, i), yind));
}
glEnd(/* GL_POINTS */);
if (centflag > 0 && xind == 0 && yind == 1) {
glColor3f(0.0, 1.0, 0.0);
glRectf(cpos1[0]-0.01, cpos1[1]-0.01, cpos1[0]+0.01, cpos1[1]+0.01);
glRectf(cpos2[0]-0.01, cpos2[1]-0.01, cpos2[0]+0.01, cpos2[1]+0.01);
}
if (centflag >= 0 || (centflag == -1 && xind == 0 && yind == 1)) {
glColor3f(0.0, 0.0, 1.0);
for (bp = vtab; bp < NthBody(vtab, 2); bp = NextBody(bp))
showcross(SelCoord(bp, xind), SelCoord(bp, yind), 0.05);
}
for (i = 1; i <= nbox; i++)
if ((xind == 0 && yind == 1) ||
(boxtab[i].velflag && (xind == 3 || yind == 3))) {
Xbox = (xind < 3 ? boxtab[i].XYbox[xind] : boxtab[i].Vbox);
Ybox = (yind < 3 ? boxtab[i].XYbox[yind] : boxtab[i].Vbox);
boxlev = (i == lrshow ? 1.0 : 0.7); // highlight current box
glColor3f((boxtab[i].disk == 1 ? 0.0 : boxlev),
(boxtab[i].disk == 2 ? 0.0 : boxlev),
(boxtab[i].disk == 0 ? 0.0 : boxlev));
showbox(Xbox, Ybox, 2 * boxtab[i].Rbox);
}
glPopMatrix();
}
// _________________________________________________
// displayball: display spin, product, or view ball.
void displayball(void)
{
matrix vmat, jmat;
real xfield, vfield, voff, xbox, ybox, vbox, rbox, hbox;
real rgb1[3] = { 0.0, 0.8, 0.8 }, rgb2[3] = { 0.8, 0.0, 0.8 };
vector spin1, spin2, xyoff, jvec, jorbit = { 0.0, 0.0, 1.0 };
facet *fp;
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(-1.5, 1.5, -1.5, 1.5, -100.0, 100.0);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glPushMatrix();
getpars(vmat, &xfield, spin1, spin2, xyoff, &vfield, &voff,
NULL, NULL, NULL, NULL, NULL, jmat);
if (1 <= lrshow && lrshow <= BXSHOW) { // lrshow is box index
xbox = boxtab[lrshow].XYbox[0] * xfield - xyoff[0];
ybox = boxtab[lrshow].XYbox[1] * xfield - xyoff[1];
vbox = boxtab[lrshow].Vbox * vfield - voff;
rbox = boxtab[lrshow].Rbox * xfield;
hbox = (boxtab[lrshow].velflag ? boxtab[lrshow].Rbox : 100.0) * vfield;
(void) loadball(vmat, xbox, ybox, vbox, rbox, hbox, boxtab[lrshow].disk);
}
if (lrshow != VWSHOW) { // show spin or prod ball
(void) showball(jmat, rgb1, rgb2, lrshow <= BXSHOW ? 0 : 1);
showspins(jmat, spin1, spin2);
MULMV(jvec, jmat, jorbit); // rotate orbital (z) axis
glColor3f(1.0, 1.0, 1.0);
glRectf(PerTrans(jvec[0] - 0.02, jvec[2]), // mark z direction
PerTrans(jvec[1] - 0.02, jvec[2]),
PerTrans(jvec[0] + 0.02, jvec[2]),
PerTrans(jvec[1] + 0.02, jvec[2]));
} else { // show view ball
fp = showball(vmat, rgb1, rgb2, 2); // get nearest facet
sprintf(msgbuf, "#yscores = %6.2f #c%6.2f #m%6.2f",
fp->view[0] * fp->view[1] > 0 ?
rlog10(fp->view[0] * fp->view[1]) / 2 : -99.0,
fp->view[0] > 0 ? rlog10(fp->view[0]) : -99.0,
fp->view[1] > 0 ? rlog10(fp->view[1]) : -99.0);
if (centflag >= 0) {
glColor3f(1.0, 1.0, 1.0);
showbox(0.0, 0.0, 0.04);
}
}
glPopMatrix();
}
// _____________________________________________________________
// loadball: find particles in box and load them onto spin ball.
bool loadball(matrix vmat, real xbox, real ybox, real vbox,
real rbox, real hbox, int disk)
{
real rhosmth, rhoface, dotsmth, dotcrit, q, max1, max2;
facet *sp;
int block, npart[2], icf;
bodyptr bp;
rhosmth = racos(1 - rsqrt(2.0) * disktol);
rhoface = racos(ball[0].midp[2]);
dotsmth = rcos(rhosmth);
dotcrit = rcos(rhosmth + rhoface);
for (sp = ball; sp < &ball[nfacet]; sp++)
sp->spin[0] = sp->spin[1] = 0.0;
block = 0; // count particle blocks
npart[0] = npart[1] = 0;
for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp))
if (AuxVec(bp)[0] != 0 || AuxVec(bp)[1] != 0 || AuxVec(bp)[2] != 0) {
block = (block % 2 == 1 ? block + 1 : block);
if ((2*disk == block || disk == 0) &&
rabs(dotvp(vmat[0], Pos(bp)) - xbox) < rbox &&
rabs(dotvp(vmat[1], Pos(bp)) - ybox) < rbox &&
rabs(dotvp(vmat[2], Vel(bp)) - vbox) < hbox) {
npart[block/2 - 1]++; // count particles in box
for (icf = 0; icf < NUMICOS; icf++) { // loop over primary facets
sp = &ball[icf * (nfacet / NUMICOS)]; // point to middle facet
if (dotvp(sp->midp, AuxVec(bp)) > dotcrit)
for (sp = sp; sp < &ball[(icf + 1) * (nfacet / NUMICOS)]; sp++)
if (dotvp(sp->midp, AuxVec(bp)) > dotsmth) {
q = racos(dotvp(sp->midp, AuxVec(bp))) / (0.5 * rhosmth);
#if defined(NO_WEIGHTING)
sp->spin[block/2 - 1] +=
(q < 1 ? 1 - 1.5*rsqr(q) + 0.75*rqbe(q) : rqbe(2 - q) / 4);
#else
sp->spin[block/2 - 1] += Aux(bp) *
(q < 1 ? 1 - 1.5*rsqr(q) + 0.75*rqbe(q) : rqbe(2 - q) / 4);
#endif
}
}
}
} else
block = (block % 2 == 0 ? block + 1 : block);
#if !defined(NO_NORMALIZE)
max1 = max2 = 0;
for (sp = ball; sp < &ball[nfacet]; sp++) {
max1 = MAX(max1, sp->spin[0]);
max2 = MAX(max2, sp->spin[1]);
}
for (sp = ball; sp < &ball[nfacet]; sp++) {
sp->spin[0] = (max1 > 0 ? sp->spin[0] / max1 : 0);
sp->spin[1] = (max2 > 0 ? sp->spin[1] / max2 : 0);
}
#endif
return (npart[0] + npart[1] > 0);
}
// _______________________________________________
// showball: display functions on surface of ball.
facet *showball(matrix tmat, real *rgb1, real *rgb2, int ind)
{
facet *fp, *fnear;
real *func, max1 = 0, max2 = 0, znear = 0, f, g;
vector mp, vp1, vp2, vp3;
for (fp = ball; fp < &ball[nfacet]; fp++) {
func = (ind == 0 ? fp->spin : ind == 1 ? fp->prod : fp->view);
max1 = MAX(max1, func[0]);
max2 = MAX(max2, func[1]);
}
for (fp = ball; fp < &ball[nfacet]; fp++) {
MULMV(mp, tmat, fp->midp);
if (mp[2] > 1.0 / PERDIST) { // is midp visible?
func = (ind == 0 ? fp->spin : ind == 1 ? fp->prod : fp->view);
f = rcbrt(max1 > 0 ? func[0] / max1 : 0);
g = rcbrt(max2 > 0 ? func[1] / max2 : 0);
glColor3f(rgb1[0] * f + rgb2[0] * g,
rgb1[1] * f + rgb2[1] * g,
rgb1[2] * f + rgb2[2] * g);
MULMV(vp1, tmat, fp->vert1);
MULMV(vp2, tmat, fp->vert2);
MULMV(vp3, tmat, fp->vert3);
glBegin(GL_TRIANGLES);
glVertex2f(PerTrans(vp1[0], vp1[2]), PerTrans(vp1[1], vp1[2]));
glVertex2f(PerTrans(vp2[0], vp2[2]), PerTrans(vp2[1], vp2[2]));
glVertex2f(PerTrans(vp3[0], vp3[2]), PerTrans(vp3[1], vp3[2]));
glEnd(/* GL_TRIANGLES */);
glColor3f(0.25, 0.25, 0.25);
glBegin(GL_POINTS);
glVertex2f(PerTrans(mp[0], mp[2]), PerTrans(mp[1], mp[2]));
glEnd(/* GL_POINTS */);
fnear = (mp[2] > znear ? fp : fnear);
znear = (mp[2] > znear ? mp[2] : znear);
}
}
return (fnear); // return nearest facet
}
// ________________________________________________
// showspins: plot spins of particles in each disk.
void showspins(matrix jmat, vector spin1, vector spin2)
{
int block = 0;
bodyptr bp;
vector jvec;
glBegin(GL_POINTS);
for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp))
if (AuxVec(bp)[0] != 0 || AuxVec(bp)[1] != 0 || AuxVec(bp)[2] != 0) {
block = (block % 2 == 1 ? block + 1 : block);
if (dotvp(AuxVec(bp), block==2 ? spin1 : spin2) > 1 - disktol) {
MULMV(jvec, jmat, AuxVec(bp));
glColor3f(block==2 ? 0.0 : 0.5, block==2 ? 0.5 : 0.0, 0.5);
glVertex2f(PerTrans(jvec[0], jvec[2]), PerTrans(jvec[1], jvec[2]));
}
} else
block = (block % 2 == 0 ? block + 1 : block);
glEnd(/* GL_POINTS */);
}
// ________________________________________________________
// displaypars: display parameters describing current view.
void displaypars(void)
{
int effsize = (zoomflag > 0 ? 2 * pansize : pansize);
matrix vmat;
real xfield, vfield, voff;
vector spin1, spin2, xyoff;
char buf[64];
if (parflag > 0) {
getpars(vmat, &xfield, spin1, spin2, xyoff, &vfield, &voff,
NULL, NULL, NULL, NULL, NULL, NULL);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(0, effsize, effsize, 0, -1.0, 1.0);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glPushMatrix();
sprintf(buf, "#wview: %.1f %.1f %.1f", 90 * idk.par[XYROT][1],
90 * idk.par[XYROT][0], 90 * idk.par[SCALE][0]);
showtext(buf, 10, 15);
sprintf(buf, "#wtime: %.2f", tnow);
showtext(buf, effsize - (8 * strlen(buf) + 10), 15);
glColor3f(1.0, 1.0, 0.3);
showtext(msgbuf, 10, 30);
if (parflag > 1) {
sprintf(buf, "#wxyoff: %.3f %.3f", xyoff[0], xyoff[1]);
showtext(buf, 10, effsize - 40);
sprintf(buf, "#wvoff: %.3f", voff);
showtext(buf, effsize - (8 * strlen(buf) + 10), effsize - 40);
}
sprintf(buf, "#cdisk1: %.1f %.1f (%.1f)",
racosD(spin1[2]), ratan2D(spin1[0], spin1[1]),
racosD(dotvp(spin1, vmat[2])));
showtext(buf, 10, effsize - 25);
sprintf(buf, "#wxfield: %.2f", xfield);
showtext(buf, effsize - (8 * strlen(buf) + 10), effsize - 25);
sprintf(buf, "#mdisk2: %.1f %.1f (%.1f)",
racosD(spin2[2]), ratan2D(spin2[0], spin2[1]),
racosD(dotvp(spin2, vmat[2])));
showtext(buf, 10, effsize - 10);
sprintf(buf, "#wvfield: %.2f", vfield);
showtext(buf, effsize - (8 * strlen(buf) + 10), effsize - 10);
glPopMatrix();
}
}
// _______________________________
// keyboard: process command keys.
void keyboard(unsigned char key, int x, int y)
{
int i;
vector XYbox;
switch (key) {
case '?':
printhelp();
break;
case 040:
case 010:
case 0177:
storehist();
idk.step += (key == 040 ? 1 : -1) * getiparam("incr");
getdata(sims[idk.isim], zrot[idk.isim], FALSE);
break;
case 'n':
case 'p':
storehist();
idk.isim = (key == 'n' ? (idk.isim < nsims - 1 ? idk.isim + 1 : 0) :
(idk.isim > 0 ? idk.isim - 1 : nsims - 1));
getjhat(jhat[idk.isim], zrot[idk.isim], FALSE);
getdata(sims[idk.isim], zrot[idk.isim], FALSE);
break;
case 'x':
case 'y':
case 'z':
storehist();
rotmatrix(idk.pmat, key == 'y' ? 90 : 0, key == 'x' ? 90 : 0, 0);
break;
case 'r':
case 'R':
storehist();
clrpars(key == 'R');
break;
case 'l':
case 'm':
disktol = disktol * rsqrt(rsqrt(key == 'l' ? 0.5 : 2.0));
seldisks();
sprintf(msgbuf, "#ydisk tol: %.4f nview: %d,%d",
disktol, nview1, nview2);
break;
case '1':
case '2':
if (1 <= lrshow && lrshow <= BXSHOW) { // current box = lrshow
boxtab[lrshow].disk = (key == '1' ? 1 : 2);
sprintf(msgbuf, "#ycurrent box selects disk %c", key);
} else { // no current box?
butbind[key == '1' ? 1 : 2] = actlast = (key == '1' ? DISK1 : DISK2);
sprintf(msgbuf, "#y%s button adjusts disk %c",
key == '1' ? "2nd (opt)" : "3rd (cmd)", key);
}
break;
case 'c':
centflag = (centflag == 1 ? -2 : centflag + 1);
if (centflag > 0 && cposflag == FALSE) {
SETV(cpos1, AuxVec(NthBody(vtab, 0))); // use current positions
SETV(cpos2, AuxVec(NthBody(vtab, 1)));
}
sprintf(msgbuf, "#ycenters %s", centflag > 0 ? "locked" : "unlocked");
break;
case 's':
butbind[1] = actlast = SCALE;
sprintf(msgbuf, "#y2rd (opt) button sets scale & rotation");
break;
case 'o':
butbind[2] = actlast = XYOFF;
sprintf(msgbuf, "#y3rd (cmd) button sets XY offsets");
break;
case 'v':
butbind[1] = actlast = VELOS;
sprintf(msgbuf, "#y2nd (opt) button sets velocities");
break;
case 't':
butbind[2] = actlast = ITRAN;
sprintf(msgbuf, "#y3rd (cmd) button sets contrast");
break;
case 'Q':
case 'q':
if (nbox < MAXBOX) {
lrshow = ++nbox; // make new box current
boxtab[lrshow].velflag = (key == 'Q');
boxtab[lrshow].disk = 0;
getpars(NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
XYbox, NULL, NULL, NULL);
XYbox[0] += 3 / ((real) pansize);
XYbox[1] -= 3 / ((real) pansize);
setpars(NULL, NULL, NULL, NULL, NULL, XYbox, NULL, NULL);
butbind[1] = BOXXY;
butbind[2] = BOXVR;
sprintf(msgbuf, "#y2nd and 3rd buttons set box %d", lrshow);
} else
sprintf(msgbuf, "#rtoo many boxes (max = %d)", MAXBOX);
break;
case 'D':
if (nbox > 0 && 1 <= lrshow && lrshow <= BXSHOW) {
for (i = lrshow; i < nbox; i++)
boxtab[i] = boxtab[i + 1];
nbox--;
lrshow = MIN(lrshow, nbox);
setpars(NULL, NULL, NULL, NULL, NULL, boxtab[lrshow].XYbox,
&boxtab[lrshow].Vbox, &boxtab[lrshow].Rbox);
sprintf(msgbuf, nbox > 0 ? "#y%d boxes left; box %d is current" :
"#yno boxes left", nbox, lrshow);
} else
sprintf(msgbuf, "#rno box to delete");
break;
case 'L':
printf("\n%3s %6s %6s %6s %6s\n",
"box", "Xbox", "Ybox", "Vbox", "Rbox");
for (i = 1; i <= nbox; i++)
printf("%2d%1s %6.3f %6.3f %6.3f %6.3f\n", i,
i == lrshow ? "*" : " ", boxtab[i].XYbox[0],
boxtab[i].XYbox[1], boxtab[i].Vbox, boxtab[i].Rbox);
break;
case 'j':
lrshow = (lrshow < nbox ? lrshow + 1 : lrshow == nbox ? PRSHOW :
lrshow == PRSHOW ? VWSHOW : XZSHOW);
butbind[0] = (lrshow == XZSHOW || lrshow == VWSHOW ? XYROT : JVIEW);
if (1 <= lrshow && lrshow <= BXSHOW) { // is there a current box?
butbind[1] = BOXXY;
butbind[2] = BOXVR;
setpars(NULL, NULL, NULL, NULL, NULL, boxtab[lrshow].XYbox,
&boxtab[lrshow].Vbox, &boxtab[lrshow].Rbox);
}
sprintf(msgbuf, "#y1st button adjusts %s",
lrshow == VWSHOW ? "view ball" :
lrshow == PRSHOW ? "product ball" :
lrshow == XZSHOW ? "viewing direction" : "spin ball");
break;
case 'J':
if (nbox > 0)
scanboxes();
else
sprintf(msgbuf, "#rno boxes to scan");
break;
case 'V':
if (nbox > 0 && centflag > 0) {
animate = TRUE;
glutTimerFunc(10, scanviews, 0);
} else
sprintf(msgbuf,
nbox == 0 ? "#rno boxes to scan" : "#rcenters not locked");
break;
case 'M':
findmax();