1 /* graphics.c
2
3 MolScript v2.1.2
4
5 Graphics: construct the geometries and call the output procedures.
6
7 Copyright (C) 1997-1998 Per Kraulis
8 6-Dec-1996 first attempts
9 10-Oct-1997 fairly finished
10 */
11
12 #include <assert.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <ctype.h>
16 #include <math.h>
17
18 #include "clib/str_utils.h"
19 #include "clib/angle.h"
20 #include "clib/extent3d.h"
21 #include "clib/hermite_curve.h"
22 #include "clib/matrix3.h"
23
24 #include "graphics.h"
25 #include "global.h"
26 #include "segment.h"
27 #include "lex.h"
28 #include "select.h"
29 #include "xform.h"
30
31 /* empirically determined factors */
32 #define HELIX_HERMITE_FACTOR 4.7
33 #define STRAND_HERMITE_FACTOR 0.5
34 #define HELIX_ALPHA (to_radians(32.0))
35 #define HELIX_BETA (to_radians(-11.0))
36
37
38 /*============================================================*/
39 boolean frame;
40 double area [4];
41 colour background_colour;
42 double window;
43 double slab;
44 boolean headlight;
45 boolean shadows;
46 double fog;
47
48 double aspect_ratio;
49 double aspect_window_x, aspect_window_y;
50
51
52 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
53 void (*output_first_plot) (void);
54 void (*output_finish_output) (void);
55 void (*output_start_plot) (void);
56 void (*output_finish_plot) (void);
57
58 void (*set_area) (void);
59 void (*set_background) (void);
60 void (*anchor_start) (char *str);
61 void (*anchor_description) (char *str);
62 void (*anchor_parameter) (char *str);
63 void (*anchor_start_geometry) (void);
64 void (*anchor_finish) (void);
65 void (*lod_start) (void);
66 void (*lod_finish) (void);
67 void (*lod_start_group) (void);
68 void (*lod_finish_group) (void);
69 void (*viewpoint_start) (char *str);
70 void (*viewpoint_output) (void);
71 void (*output_directionallight) (void);
72 void (*output_pointlight) (void);
73 void (*output_spotlight) (void);
74 void (*output_comment) (char *str);
75
76 void (*output_coil) (void);
77 void (*output_cylinder) (vector3 *v1, vector3 *v2);
78 void (*output_helix) (void);
79 void (*output_label) (vector3 *p, char *label, colour *c);
80 void (*output_line) (boolean polylines);
81 void (*output_sphere) (at3d *at, double radius);
82 void (*output_stick) (vector3 *v1, vector3 *v2,
83 double r1, double r2, colour *c);
84 void (*output_strand) (void);
85
86 void (*output_start_object) (void);
87 void (*output_object) (int code, vector3 *triplets, int count);
88 void (*output_finish_object) (void);
89
90 void (*output_pickable) (at3d *atom);
91
92
93 /*------------------------------------------------------------*/
94 static void
msg_chain(char * type,mol3d_chain * ch)95 msg_chain (char *type, mol3d_chain *ch)
96 {
97 assert (type);
98 assert (ch);
99
100 if (message_mode)
101 fprintf (stderr, "%s from %s to %s (%i residues)\n", type,
102 (ch->residues[0])->name, (ch->residues[ch->length-1])->name,
103 ch->length);
104 }
105
106
107 /*------------------------------------------------------------*/
108 static void
set_aspect_ratio(void)109 set_aspect_ratio (void)
110 {
111 aspect_ratio = (area[3] - area[1]) / (area[2] - area[0]);
112 if (aspect_ratio >= 1.0) {
113 aspect_window_x = window;
114 aspect_window_y = aspect_ratio * window;
115 } else {
116 aspect_window_x = window / aspect_ratio;
117 aspect_window_y = window;
118 }
119 }
120
121
122 /*------------------------------------------------------------*/
123 void
graphics_plot_init(void)124 graphics_plot_init (void)
125 {
126 frame = TRUE;
127 window = -1.0;
128 slab = -1.0;
129 headlight = TRUE;
130 shadows = FALSE;
131 fog = 0.0;
132 ext3d_initialize();
133 }
134
135
136 /*------------------------------------------------------------*/
137 void
set_area_values(double xlo,double ylo,double xhi,double yhi)138 set_area_values (double xlo, double ylo, double xhi, double yhi)
139 {
140 area[0] = xlo;
141 area[1] = ylo;
142 area[2] = xhi;
143 area[3] = yhi;
144 set_aspect_ratio();
145 }
146
147
148 /*------------------------------------------------------------*/
149 void
set_window(void)150 set_window (void)
151 {
152 assert (dstack_size == 1);
153
154 if (dstack[0] <= 0.0) {
155 yyerror ("invalid window value");
156 return;
157 }
158
159 window = dstack[0] / 2.0;
160 clear_dstack();
161
162 set_aspect_ratio();
163
164 assert (dstack_size == 0);
165 }
166
167
168 /*------------------------------------------------------------*/
169 void
set_slab(void)170 set_slab (void)
171 {
172 assert (dstack_size == 1);
173
174 if (dstack[0] <= 0.0) {
175 yyerror ("invalid slab value");
176 return;
177 }
178
179 slab = dstack[0] / 2.0;
180 clear_dstack();
181
182 assert (dstack_size == 0);
183 }
184
185
186 /*------------------------------------------------------------*/
187 void
set_fog(void)188 set_fog (void)
189 {
190 assert (dstack_size == 1);
191
192 if (dstack[0] < 0.0) {
193 yyerror ("invalid fog value");
194 return;
195 }
196
197 fog = dstack[0];
198 clear_dstack();
199
200 assert (dstack_size == 0);
201 }
202
203
204 /*------------------------------------------------------------*/
205 void
set_extent(void)206 set_extent (void)
207 {
208 vector3 center, size;
209 double d;
210
211 ext3d_get_center_size (¢er, &size);
212 v3_scale (&size, 0.5);
213
214 if (window < 0.0) {
215 if (aspect_ratio < 1.0) {
216 size.x *= aspect_ratio;
217 } else {
218 size.y /= aspect_ratio;
219 }
220 window = fabs (center.x - size.x);
221 d = fabs (center.x + size.x);
222 if (d > window) window = d;
223 d = fabs (center.y - size.y);
224 if (d > window) window = d;
225 d = fabs (center.y + size.y);
226 if (d > window) window = d;
227 window += 2.0; /* make fit less snug */
228 if (message_mode)
229 fprintf (stderr, "setting window to %.2f\n", 2.0 * window);
230 }
231
232 if (slab < 0.0) {
233 slab = fabs (center.z - size.z);
234 d = fabs (center.z + size.z);
235 if (d > slab) slab = d;
236 if (message_mode)
237 fprintf (stderr, "setting slab to %.2f\n", 2.0 * slab);
238 }
239 }
240
241
242 /*------------------------------------------------------------*/
243 int
outside_extent_radius(vector3 * v,double radius)244 outside_extent_radius (vector3 *v, double radius)
245 {
246 assert (v);
247 assert (radius >= 0.0);
248
249 if (window > 0.0) {
250 if (v->x + radius < -aspect_window_x) return TRUE;
251 if (v->x - radius > aspect_window_x) return TRUE;
252 if (v->y + radius < -aspect_window_y) return TRUE;
253 if (v->y - radius > aspect_window_y) return TRUE;
254 }
255
256 if (slab > 0.0) {
257 if (v->z + radius < -slab) return TRUE;
258 if (v->z - radius > slab) return TRUE;
259 }
260
261 return FALSE;
262 }
263
264
265 /*------------------------------------------------------------*/
266 int
outside_extent_2v(vector3 * v1,vector3 * v2)267 outside_extent_2v (vector3 *v1, vector3 *v2)
268 {
269 assert (v1);
270 assert (v2);
271
272 if (window > 0.0) {
273 if ((v1->x < -aspect_window_x) && (v2->x < -aspect_window_x)) return TRUE;
274 if ((v1->x > aspect_window_x) && (v2->x > aspect_window_x)) return TRUE;
275 if ((v1->y < -aspect_window_y) && (v2->y < -aspect_window_y)) return TRUE;
276 if ((v1->y > aspect_window_y) && (v2->y > aspect_window_y)) return TRUE;
277 }
278
279 if (slab > 0.0) {
280 if ((v1->z < -slab) && (v2->z < -slab)) return TRUE;
281 if ((v1->z > slab) && (v2->z > slab)) return TRUE;
282 }
283
284 return FALSE;
285 }
286
287
288 /*------------------------------------------------------------*/
289 double
depthcue(double depth,state * st)290 depthcue (double depth, state *st)
291 {
292 double factor;
293
294 assert (st);
295
296 if (slab <= 0.0) return 1.0;
297
298 factor = depth / slab;
299 if (factor < -1.0) {
300 factor = -1.0;
301 } else if (factor > 1.0) {
302 factor = 1.0;
303 }
304 factor = (factor / 2.0 - 0.5) * st->depthcue + 1.0;
305
306 return factor;
307 }
308
309
310 /*------------------------------------------------------------*/
311 void
ball_and_stick(int single_selection)312 ball_and_stick (int single_selection)
313 {
314 at3d **atoms1, **atoms2;
315 int atom_count1, atom_count2, slot1, slot2, start;
316 vector3 *v1, *v2;
317 double radius, dist;
318
319 if (single_selection) {
320
321 assert (count_atom_selections() == 1);
322
323 atoms1 = select_atom_list (&atom_count1);
324
325 if (message_mode)
326 fprintf (stderr, "%i atoms selected for ball-and-stick\n", atom_count1);
327
328 if (atoms1 == NULL) return;
329
330 for (slot1 = 0; slot1 < atom_count1; slot1++) { /* balls output */
331 radius = 0.25 * atoms1[slot1]->radius;
332 if (radius > 0.0) {
333 output_sphere (atoms1[slot1], radius);
334 ext3d_update (&(atoms1[slot1]->xyz), radius);
335 }
336 }
337
338 atoms2 = atoms1;
339 atom_count2 = atom_count1;
340
341 } else {
342
343 assert (count_atom_selections() == 2);
344
345 atoms1 = select_atom_list (&atom_count1);
346 atoms2 = select_atom_list (&atom_count2);
347
348 if (message_mode) {
349 fprintf (stderr, "%i atoms in first set and %i atoms in second set for ball-and-stick\n", atom_count1, atom_count2);
350 }
351
352 if (atoms1 == NULL) {
353 if (atoms2 != NULL) free (atoms2);
354 return;
355 } else if (atoms2 == NULL) {
356 free (atoms1);
357 return;
358 }
359 }
360
361 assert (count_atom_selections() == 0);
362
363 if (current_state->colourparts) { /* sticks output, atom colour */
364
365 vector3 middle;
366
367 for (slot1 = 0; slot1 < atom_count1; slot1++) {
368 v1 = &(atoms1[slot1]->xyz);
369
370 if (single_selection) {
371 start = slot1 + 1;
372 } else {
373 start = 0;
374 }
375
376 for (slot2 = start; slot2 < atom_count2; slot2++) {
377 if (atoms1[slot1] == atoms2[slot2]) continue;
378
379 v2 = &(atoms2[slot2]->xyz);
380
381 dist = v3_distance (v1, v2);
382 if ((dist > current_state->bonddistance) || (dist < 0.001)) continue;
383
384 if (colour_unequal (&(atoms1[slot1]->colour),
385 &(atoms2[slot2]->colour))) {
386 v3_middle (&middle, v1, v2);
387 output_stick (v1, &middle, 0.25 * atoms1[slot1]->radius, -1.0,
388 &(atoms1[slot1]->colour));
389 output_stick (&middle, v2, -1.0, 0.25 * atoms2[slot2]->radius,
390 &(atoms2[slot2]->colour));
391 } else {
392 output_stick (v1, v2,
393 0.25 * atoms1[slot1]->radius,
394 0.25 * atoms2[slot2]->radius,
395 &(atoms1[slot1]->colour));
396 }
397 ext3d_update (v1, current_state->stickradius);
398 ext3d_update (v2, current_state->stickradius);
399 }
400 }
401
402 } else { /* sticks output, overall colour */
403
404 for (slot1 = 0; slot1 < atom_count1; slot1++) {
405 v1 = &(atoms1[slot1]->xyz);
406
407 if (single_selection) {
408 start = slot1 + 1;
409 } else {
410 start = 0;
411 }
412
413 for (slot2 = start; slot2 < atom_count2; slot2++) {
414 if (atoms1[slot1] == atoms2[slot2]) continue;
415
416 v2 = &(atoms2[slot2]->xyz);
417
418 dist = v3_distance (v1, v2);
419 if ((dist > current_state->bonddistance) || (dist < 0.001)) continue;
420
421 output_stick (v1, v2,
422 0.25 * atoms1[slot1]->radius,
423 0.25 * atoms2[slot2]->radius,
424 NULL);
425 ext3d_update (v1, current_state->stickradius);
426 ext3d_update (v2, current_state->stickradius);
427 }
428 }
429 }
430
431 if (single_selection) {
432 free (atoms1);
433 } else {
434 free (atoms2);
435 free (atoms1);
436 }
437
438 assert (count_atom_selections() == 0);
439 }
440
441
442 /*------------------------------------------------------------*/
443 void
bonds(int single_selection)444 bonds (int single_selection)
445 {
446 at3d **atoms1, **atoms2;
447 int atom_count1, atom_count2, slot1, slot2, start;
448 vector3 *v1, *v2;
449 colour *col1;
450 line_segment *ls;
451 double dist;
452
453 if (single_selection) {
454
455 assert (count_atom_selections() == 1);
456 atoms1 = select_atom_list (&atom_count1);
457 assert (count_atom_selections() == 0);
458
459 if (message_mode)
460 fprintf (stderr, "%i atoms selected for bonds\n", atom_count1);
461
462 if (atoms1 == NULL) return;
463
464 atoms2 = atoms1;
465 atom_count2 = atom_count1;
466
467 for (slot1 = 0; slot1 < atom_count1; slot1++) atoms1[slot1]->rval = 0.0;
468
469 } else {
470
471 assert (count_atom_selections() == 2);
472 atoms1 = select_atom_list (&atom_count1);
473 atoms2 = select_atom_list (&atom_count2);
474 assert (count_atom_selections() == 0);
475
476 if (message_mode) {
477 fprintf (stderr,
478 "%i atoms in first set and %i atoms in second set for bonds\n",
479 atom_count1, atom_count2);
480 }
481
482 if (atoms1 == NULL) {
483 if (atoms2 != NULL) free (atoms2);
484 return;
485 } else if (atoms2 == NULL) {
486 free (atoms1);
487 return;
488 }
489 }
490
491 line_segment_init();
492
493 if (current_state->colourparts) {
494 vector3 middle;
495
496 for (slot1 = 0; slot1 < atom_count1; slot1++) {
497 v1 = &(atoms1[slot1]->xyz);
498 col1 = &(atoms1[slot1]->colour);
499
500 if (single_selection) {
501 start = slot1 + 1;
502 } else {
503 start = 0;
504 }
505
506 for (slot2 = start; slot2 < atom_count2; slot2++) {
507 v2 = &(atoms2[slot2]->xyz);
508 if (v2 == v1) continue;
509
510 dist = v3_distance (v1, v2);
511 if ((dist > current_state->bonddistance) || (dist < 0.001)) continue;
512
513 if (colour_unequal (col1, &(atoms2[slot2]->colour))) {
514 v3_middle (&middle, v1, v2);
515
516 ls = line_segment_next();
517 ls->p = *v1;
518 ls->c = *col1;
519 ls = line_segment_next();
520 ls->p = middle;
521
522 ls = line_segment_next();
523 ls->p = middle;
524 ls->c = atoms2[slot2]->colour;
525 ls = line_segment_next();
526 ls->p = *v2;
527
528 } else {
529 ls = line_segment_next();
530 ls->p = *v1;
531 ls->c = *col1;
532 ls = line_segment_next();
533 ls->p = *v2;
534 }
535
536 ext3d_update (v1, 0.0);
537 ext3d_update (v2, 0.0);
538
539 atoms1[slot1]->rval++;
540 atoms2[slot2]->rval++;
541 }
542 }
543
544 } else {
545
546 for (slot1 = 0; slot1 < atom_count1; slot1++) {
547 v1 = &(atoms1[slot1]->xyz);
548
549 if (single_selection) {
550 start = slot1 + 1;
551 } else {
552 start = 0;
553 }
554
555 for (slot2 = start; slot2 < atom_count2; slot2++) {
556 v2 = &(atoms2[slot2]->xyz);
557 if (v2 == v1) continue;
558
559 dist = v3_distance (v1, v2);
560 if ((dist > current_state->bonddistance) || (dist < 0.001)) continue;
561
562 ls = line_segment_next();
563 ls->p = *v1;
564 ls = line_segment_next();
565 ls->p = *v2;
566
567 ext3d_update (v1, 0.0);
568 ext3d_update (v2, 0.0);
569 atoms1[slot1]->rval++;
570 atoms2[slot2]->rval++;
571 }
572
573 if (line_segment_count > 0) line_segments->c = current_state->linecolour;
574 }
575 }
576
577 if (single_selection && (current_state->bondcross != 0.0)) {
578 double radius = 0.5 * current_state->bondcross;
579
580 if (! current_state->colourparts) col1 = NULL;
581
582 for (slot1 = 0; slot1 < atom_count1; slot1++) {
583 if (atoms1[slot1]->rval != 0.0) continue;
584
585 atoms1[slot1]->rval = 1.0;
586
587 v1 = &(atoms1[slot1]->xyz);
588 if (current_state->colourparts) col1 = &(atoms1[slot1]->colour);
589
590 ls = line_segment_next();
591 v3_sum_scaled (&(ls->p), v1, radius, &xaxis);
592 if (col1) ls->c = *col1;
593 ls = line_segment_next();
594 v3_sum_scaled (&(ls->p), v1, -radius, &xaxis);
595
596 ls = line_segment_next();
597 v3_sum_scaled (&(ls->p), v1, radius, &yaxis);
598 if (col1) ls->c = *col1;
599 ls = line_segment_next();
600 v3_sum_scaled (&(ls->p), v1, -radius, &yaxis);
601
602 ls = line_segment_next();
603 v3_sum_scaled (&(ls->p), v1, radius, &zaxis);
604 if (col1) ls->c = *col1;
605 ls = line_segment_next();
606 v3_sum_scaled (&(ls->p), v1, -radius, &zaxis);
607
608 ext3d_update (v1, 0.0);
609 }
610 }
611
612 output_line (FALSE);
613
614 if (output_pickable) {
615 if (single_selection) {
616 for (slot1 = 0; slot1 < atom_count1; slot1++) {
617 if (atoms1[slot1]->rval != 0.0) {
618 output_pickable (atoms1[slot1]);
619 }
620 }
621 } else {
622 }
623 }
624
625 if (single_selection) {
626 free (atoms1);
627 } else {
628 free (atoms2);
629 free (atoms1);
630 }
631
632 assert (count_atom_selections() == 0);
633 }
634
635
636 /*------------------------------------------------------------*/
637 static vector3 *
get_atom_positions(mol3d_chain * ch)638 get_atom_positions (mol3d_chain *ch)
639 {
640 int slot;
641 vector3 *points;
642
643 assert (ch);
644
645 points = malloc (ch->length * sizeof (vector3));
646 for (slot = 0; slot < ch->length; slot++) {
647 *(points + slot) = ch->atoms[slot]->xyz;
648 }
649
650 return points;
651 }
652
653
654 /*------------------------------------------------------------*/
655 static void
priestle_smoothing(vector3 * points,int length,int steps)656 priestle_smoothing (vector3 *points, int length, int steps)
657 {
658 int sm, slot;
659 vector3 *ptmp = malloc (length * sizeof (vector3));
660
661 assert (points);
662
663 for (sm = 0; sm < steps; sm++) {
664
665 for (slot = 1; slot < length-1; slot++) {
666 v3_middle (ptmp + slot, points + slot - 1, points + slot + 1);
667 v3_middle (ptmp + slot, ptmp + slot, points + slot);
668 }
669
670 for (slot = 1; slot < length-1; slot++) {
671 *(points + slot) = *(ptmp + slot);
672 }
673 }
674
675 free (ptmp);
676 }
677
678
679 /*------------------------------------------------------------*/
680 void
coil(int is_peptide_chain,int smoothing)681 coil (int is_peptide_chain, int smoothing)
682 {
683 char *atomname, *coilname;
684 double chain_distance, t;
685 mol3d_chain *first_ch, *ch;
686 vector3 *points;
687 res3d *res;
688 at3d *first, *last;
689 vector3 vec1, vec2;
690 colour *col;
691 int slot, segment;
692 double hermite_factor = 0.5 * current_state->splinefactor;
693 int segments = current_state->segments;
694
695 assert (count_residue_selections() == 1);
696
697 /* colourparts requires even number of seg's */
698 if (current_state->colourparts && (segments % 2)) segments++;
699
700 if (is_peptide_chain) {
701 atomname = PEPTIDE_CHAIN_ATOMNAME;
702 chain_distance = PEPTIDE_DISTANCE;
703 first_ch = get_peptide_chains();
704 if (smoothing) {
705 coilname = "coil";
706 } else {
707 coilname = "turn";
708 }
709
710 } else {
711 atomname = NUCLEOTIDE_CHAIN_ATOMNAME;
712 chain_distance = NUCLEOTIDE_DISTANCE;
713 coilname = "double-helix";
714 first_ch = get_nucleotide_chains();
715 }
716
717 for (ch = first_ch; ch; ch = ch->next) {
718 if (ch->length < 2) continue;
719
720 msg_chain (coilname, ch);
721 points = get_atom_positions (ch);
722
723 first = NULL; /* find atom before or in first residue */
724 res = ch->residues[0]->prev;
725 if (res) {
726 first = at3d_lookup (res, atomname);
727 if (first) {
728 if (v3_distance (&(first->xyz), points) >=
729 chain_distance) first = NULL;
730 }
731 }
732 if (first == NULL) first = ch->atoms[0];
733
734 last = NULL; /* find atom after or in last residue */
735 res = ch->residues[ch->length-1]->next;
736 if (res) {
737 last = at3d_lookup (res, atomname);
738 if (last) {
739 if (v3_distance (&(last->xyz), points + ch->length-1) >=
740 chain_distance) last = NULL;
741 }
742 }
743 if (last == NULL) last = ch->atoms[ch->length-1];
744
745 if (smoothing)
746 priestle_smoothing (points, ch->length, current_state->smoothsteps);
747
748 if (current_state->colourparts) {
749 col = &(ch->residues[0]->colour);
750 } else {
751 col = NULL;
752 }
753
754 if (current_state->coilradius < 0.01) { /* coil rendered as line */
755 line_segment *ls;
756
757 line_segment_init();
758 ls = line_segment_next();
759 ls->new = TRUE;
760 ls->p = *points;
761 if (col) {
762 ls->c = *col;
763 } else {
764 ls->c = current_state->planecolour;
765 }
766 ext3d_update (points, 0.0);
767
768 v3_difference (&vec2, points + 1, &(first->xyz));
769 if (first == ch->atoms[0]) v3_scale (&vec2, 2.0);
770 v3_scale (&vec2, hermite_factor);
771
772 for (slot = 0; slot < ch->length - 1; slot++) {
773
774 vec1 = vec2;
775 if (slot == ch->length - 2) {
776 v3_difference (&vec2, &(last->xyz), points + ch->length - 2);
777 if (last == ch->atoms[ch->length-1]) v3_scale (&vec2, 2.0);
778 } else {
779 v3_difference (&vec2, points + slot + 2, points + slot);
780 }
781 v3_scale (&vec2, hermite_factor);
782
783 hermite_set (points + slot, points + slot + 1, &vec1, &vec2);
784 for (segment = 1; segment < segments; segment++) {
785 t = ((double) segment) / ((double) segments);
786 ls = line_segment_next();
787 hermite_get (&(ls->p), t);
788 if (current_state->colourparts &&
789 (segment == segments / 2) &&
790 colour_unequal (col, &(ch->residues[slot+1]->colour))) {
791 ls->c = *col;
792 ls = line_segment_next();
793 ls->new = TRUE;
794 ls->p = (ls - 1)->p;
795 col = &(ch->residues[slot+1]->colour);
796 }
797 if (col) ls->c = *col;
798 ext3d_update (&(ls->p), 0.0);
799 }
800
801 ls = line_segment_next();
802 ls->p = *(points + slot + 1);
803 if (col) ls->c = *col;
804 ext3d_update (points + slot + 1, 0.0);
805 }
806
807 output_line (TRUE);
808
809 } else { /* coil rendered as solid */
810 coil_segment *cs;
811 vector3 side = {1.0, 0.0, 0.0};
812 vector3 normal, dir;
813 double radius = current_state->coilradius / sqrt (2.0);
814
815 coil_segment_init(); /* coil path and colour*/
816 cs = coil_segment_next();
817 cs->p = *points;
818 if (col) {
819 cs->c = *col;
820 } else {
821 cs->c = current_state->planecolour;
822 }
823 ext3d_update (&(cs->p), radius);
824
825 v3_difference (&vec2, points + 1, &(first->xyz));
826 if (first == ch->atoms[0]) v3_scale (&vec2, 2.0);
827 v3_scale (&vec2, hermite_factor);
828
829 for (slot = 0; slot < ch->length - 1; slot++) {
830
831 vec1 = vec2;
832 if (slot == ch->length - 2) {
833 v3_difference (&vec2, &(last->xyz), points + ch->length - 2);
834 if (last == ch->atoms[ch->length-1]) v3_scale (&vec2, 2.0);
835 } else {
836 v3_difference (&vec2, points + slot + 2, points + slot);
837 }
838 v3_scale (&vec2, hermite_factor);
839
840 hermite_set (points + slot, points + slot + 1, &vec1, &vec2);
841 for (segment = 1; segment < segments; segment++) {
842 t = ((double) segment) / ((double) segments);
843 cs = coil_segment_next();
844 hermite_get (&(cs->p), t);
845 if (current_state->colourparts &&
846 (segment == segments / 2)) col = &(ch->residues[slot+1]->colour);
847 if (col) cs->c = *col;
848 ext3d_update (&(cs->p), radius);
849 }
850
851 cs = coil_segment_next();
852 cs->p = *(points +slot + 1);
853 if (col) cs->c = *col;
854 ext3d_update (points + slot + 1, radius);
855 }
856 /* coil plane coordinates */
857 v3_difference (&dir, &(coil_segments[1].p), &(coil_segments[0].p));
858 v3_normalize (&dir);
859 v3_cross_product (&normal, &side, &dir);
860 if (v3_length (&normal) < 1.0e-6) {
861 v3_initialize (&side, 0.0, 1.0, 0.0);
862 v3_cross_product (&normal, &side, &dir);
863 if (v3_length (&normal) < 1.0e-6) {
864 v3_initialize (&side, 0.0, 0.0, 1.0);
865 v3_cross_product (&normal, &side, &dir);
866 }
867 }
868 v3_normalize (&normal);
869 v3_cross_product (&side, &dir, &normal);
870 v3_normalize (&side);
871
872 cs = coil_segments;
873 v3_sum_scaled (&(cs->p1), &(cs->p), radius, &normal);
874 v3_add_scaled (&(cs->p1), radius, &side);
875 v3_sum_scaled (&(cs->p2), &(cs->p), -radius, &normal);
876 v3_add_scaled (&(cs->p2), radius, &side);
877 v3_sum_scaled (&(cs->p3), &(cs->p), -radius, &normal);
878 v3_add_scaled (&(cs->p3), -radius, &side);
879 v3_sum_scaled (&(cs->p4), &(cs->p), radius, &normal);
880 v3_add_scaled (&(cs->p4), -radius, &side);
881
882 for (slot = 1; slot < coil_segment_count - 1; slot++) {
883 v3_difference (&dir, &(coil_segments[slot+1].p),
884 &(coil_segments[slot-1].p));
885 v3_cross_product (&side, &dir, &normal);
886 v3_normalize (&side);
887 v3_cross_product (&normal, &side, &dir);
888 v3_normalize (&normal);
889
890 cs = coil_segments + slot;
891 v3_sum_scaled (&(cs->p1), &(cs->p), radius, &normal);
892 v3_add_scaled (&(cs->p1), radius, &side);
893 v3_sum_scaled (&(cs->p2), &(cs->p), -radius, &normal);
894 v3_add_scaled (&(cs->p2), radius, &side);
895 v3_sum_scaled (&(cs->p3), &(cs->p), -radius, &normal);
896 v3_add_scaled (&(cs->p3), -radius, &side);
897 v3_sum_scaled (&(cs->p4), &(cs->p), radius, &normal);
898 v3_add_scaled (&(cs->p4), -radius, &side);
899 }
900
901 v3_difference (&dir, &(coil_segments[coil_segment_count-1].p),
902 &(coil_segments[coil_segment_count-2].p));
903 v3_cross_product (&side, &dir, &normal);
904 v3_normalize (&side);
905 v3_cross_product (&normal, &side, &dir);
906 v3_normalize (&normal);
907
908 cs = coil_segments + coil_segment_count - 1;
909 v3_sum_scaled (&(cs->p1), &(cs->p), radius, &normal);
910 v3_add_scaled (&(cs->p1), radius, &side);
911 v3_sum_scaled (&(cs->p2), &(cs->p), -radius, &normal);
912 v3_add_scaled (&(cs->p2), radius, &side);
913 v3_sum_scaled (&(cs->p3), &(cs->p), -radius, &normal);
914 v3_add_scaled (&(cs->p3), -radius, &side);
915 v3_sum_scaled (&(cs->p4), &(cs->p), radius, &normal);
916 v3_add_scaled (&(cs->p4), -radius, &side);
917
918 /* coil plane normals */
919 for (slot = 0; slot < coil_segment_count; slot++) {
920 cs = coil_segments + slot;
921 v3_difference (&(cs->n1), &(cs->p1), &(cs->p));
922 v3_normalize (&(cs->n1));
923 v3_difference (&(cs->n2), &(cs->p2), &(cs->p));
924 v3_normalize (&(cs->n2));
925 v3_difference (&(cs->n3), &(cs->p3), &(cs->p));
926 v3_normalize (&(cs->n3));
927 v3_difference (&(cs->n4), &(cs->p4), &(cs->p));
928 v3_normalize (&(cs->n4));
929 }
930
931 output_coil();
932 }
933
934 free (points);
935 }
936
937 if (first_ch) mol3d_chain_delete (first_ch);
938
939 assert (count_residue_selections() == 0);
940 }
941
942
943 /*------------------------------------------------------------*/
944 void
cpk(void)945 cpk (void)
946 {
947 at3d **atoms;
948 int atom_count, slot;
949 double radius;
950
951 assert (count_atom_selections() == 1);
952
953 atoms = select_atom_list (&atom_count);
954
955 assert (count_atom_selections() == 0);
956
957 if (message_mode)
958 fprintf (stderr, "%i atoms selected for cpk\n", atom_count);
959 if (atoms == NULL) return;
960
961 for (slot = 0; slot < atom_count; slot++) {
962 radius = atoms[slot]->radius;
963 if (radius > 0.0) {
964 output_sphere (atoms[slot], radius);
965 ext3d_update (&(atoms[slot]->xyz), radius);
966 }
967 }
968
969 free (atoms);
970 }
971
972
973 /*------------------------------------------------------------*/
974 void
cylinder(void)975 cylinder (void)
976 {
977 mol3d_chain *first_ch, *ch;
978 vector3 start, finish, dir, vec, top;
979
980 assert (count_residue_selections() == 1);
981
982 first_ch = get_peptide_chains();
983
984 for (ch = first_ch; ch; ch = ch->next) {
985 if (ch->length < 3) continue;
986
987 msg_chain ("cylinder", ch);
988
989 v3_middle (&start, &(ch->atoms[0]->xyz), &(ch->atoms[2]->xyz));
990 v3_middle (&finish, &(ch->atoms[ch->length-3]->xyz),
991 &(ch->atoms[ch->length-1]->xyz));
992
993 if (ch->length == 3) { /* compute helix axis as for ordinary helix */
994 vector3 vec1, vec2, rvec;
995
996 v3_difference (&vec, &(ch->atoms[2]->xyz), &(ch->atoms[0]->xyz));
997 v3_normalize (&vec);
998
999 v3_difference (&vec1, &(ch->atoms[1]->xyz), &(ch->atoms[0]->xyz));
1000 v3_difference (&vec2, &(ch->atoms[2]->xyz), &(ch->atoms[1]->xyz));
1001 v3_cross_product (&rvec, &vec1, &vec2);
1002 v3_normalize (&rvec);
1003
1004 v3_scaled (&dir, cos (HELIX_ALPHA), &rvec);
1005 v3_add_scaled (&dir, sin (HELIX_ALPHA), &vec);
1006
1007 } else { /* compute helix axis from terminii coords */
1008 v3_difference (&dir, &finish, &start);
1009 }
1010
1011 v3_normalize (&dir);
1012
1013 v3_difference (&vec, &(ch->atoms[2]->xyz), &(ch->atoms[0]->xyz));
1014 v3_scaled (&top, - 0.5 * v3_length (&vec) *
1015 cos (v3_angle (&vec, &dir)), &dir);
1016 v3_add (&start, &top);
1017
1018 v3_difference (&vec, &(ch->atoms[ch->length-1]->xyz),
1019 &(ch->atoms[ch->length-3]->xyz));
1020 v3_scaled (&top, 0.5 * v3_length (&vec) *
1021 cos (v3_angle (&vec, &dir)), &dir);
1022 v3_add (&finish, &top);
1023
1024 output_cylinder (&start, &finish);
1025 ext3d_update (&start, current_state->cylinderradius);
1026 ext3d_update (&finish, current_state->cylinderradius);
1027 }
1028
1029 if (first_ch) mol3d_chain_delete (first_ch);
1030
1031 assert (count_residue_selections() == 0);
1032 }
1033
1034
1035 /*------------------------------------------------------------*/
1036 void
helix(void)1037 helix (void)
1038 {
1039 mol3d_chain *first_ch, *ch;
1040 vector3 *points, *axes, *tangents;
1041 vector3 cvec, rvec, vec1, vec2, pos, dir;
1042 res3d *res;
1043 at3d *ca_first, *ca_last, *at;
1044 int slot, segment;
1045 double halfwidth, t;
1046 colour *col;
1047 helix_segment *hs;
1048 double coilradius = current_state->coilradius;
1049 int segments = current_state->segments;
1050
1051 assert (count_residue_selections() == 1);
1052
1053 first_ch = get_peptide_chains();
1054
1055 /* colourparts requires even number of seg's */
1056 if (current_state->colourparts && (segments % 2)) segments++;
1057
1058 if (coilradius < 0.01) coilradius = 0.01; /* not too thin */
1059
1060 for (ch = first_ch; ch; ch = ch->next) {
1061 if (ch->length < 3) continue;
1062
1063 msg_chain ("helix", ch);
1064
1065 points = get_atom_positions (ch); /* helix axis and tangent vectors */
1066 axes = malloc (ch->length * sizeof (vector3));
1067 tangents = malloc (ch->length * sizeof (vector3));
1068
1069 for (slot = 1; slot < ch->length - 1; slot++) {
1070 /* helix direction vector at (i) */
1071 v3_difference (&cvec, points + slot + 1, points + slot - 1);
1072 v3_normalize (&cvec);
1073 /* normal vector for plane (i-1),(i),(i+1) */
1074 v3_difference (&vec1, points + slot, points + slot - 1);
1075 v3_difference (&vec2, points + slot + 1, points + slot);
1076 v3_cross_product (&rvec, &vec1, &vec2);
1077 v3_normalize (&rvec);
1078 /* helix axis at (i) */
1079 v3_scaled (&vec1, cos (HELIX_ALPHA), &rvec);
1080 v3_scaled (&vec2, sin (HELIX_ALPHA), &cvec);
1081 v3_sum (axes + slot, &vec1, &vec2);
1082 /* helix tangent at (i) */
1083 v3_scaled (&vec1, cos (HELIX_BETA), &cvec);
1084 v3_scaled (&vec2, sin (HELIX_BETA), &rvec);
1085 v3_sum (tangents + slot, &vec1, &vec2);
1086 v3_scale (tangents + slot, HELIX_HERMITE_FACTOR);
1087 }
1088
1089 ca_first = ch->atoms[0]; /* find CA before and after chain, if any */
1090 res = ch->residues[0]->prev; /* used for helix terminii tangents */
1091 if (res) {
1092 at = at3d_lookup (res, PEPTIDE_CHAIN_ATOMNAME);
1093 if (at &&
1094 v3_distance (points, &(at->xyz)) <= PEPTIDE_DISTANCE) ca_first = at;
1095 }
1096
1097 ca_last = ch->atoms[ch->length-1];
1098 res = ch->residues[ch->length-1]->next;
1099 if (res) {
1100 at = at3d_lookup (res, PEPTIDE_CHAIN_ATOMNAME);
1101 if (at &&
1102 v3_distance (points + ch->length - 1,
1103 &(at->xyz)) <= PEPTIDE_DISTANCE) ca_last = at;
1104 }
1105 /* helix terminii axes */
1106 *(axes) = *(axes + 1);
1107 *(axes + ch->length - 1) = *(axes + ch->length - 2);
1108 /* helix terminii tangents */
1109 v3_difference (&vec1, points + 1, &(ca_first->xyz));
1110 v3_normalize (&vec1);
1111 v3_scaled (tangents, HELIX_HERMITE_FACTOR, &vec1);
1112 v3_difference (&vec1, &(ca_last->xyz), points + ch->length - 2);
1113 v3_normalize (&vec1);
1114 v3_scaled (tangents + ch->length - 1, HELIX_HERMITE_FACTOR, &vec1);
1115
1116 if (current_state->colourparts) {
1117 col = &(ch->residues[0]->colour);
1118 } else {
1119 col = NULL;
1120 }
1121
1122 helix_segment_init();
1123 hs = helix_segment_next();
1124 v3_sum_scaled (&(hs->p1), points, coilradius, axes);
1125 v3_sum_scaled (&(hs->p2), points, -coilradius, axes);
1126 hs->a = *axes;
1127 v3_cross_product (&(hs->n), tangents, axes);
1128 v3_normalize (&(hs->n));
1129 if (col) hs->c = *col;
1130 ext3d_update (points, coilradius);
1131
1132 /* helix start segments */
1133 pos = *points;
1134 hermite_set (points, points + 1, tangents, tangents + 1);
1135
1136 for (segment = 1; segment < segments; segment++) {
1137 t = ((double) segment) / ((double) segments);
1138 hermite_get (&pos, t);
1139 halfwidth = coilradius +
1140 (0.5 * current_state->helixwidth - coilradius) *
1141 0.5 * (- cos (ANGLE_PI * t) + 1.0);
1142 hs = helix_segment_next();
1143 v3_sum_scaled (&(hs->p1), &pos, halfwidth, axes);
1144 v3_sum_scaled (&(hs->p2), &pos, -halfwidth, axes);
1145 hs->a = *axes;
1146 hermite_get_tangent (&dir, t);
1147 v3_cross_product (&(hs->n), &dir, axes);
1148 v3_normalize (&(hs->n));
1149 if (current_state->colourparts &&
1150 (segment == segments / 2)) col = &(ch->residues[1]->colour);
1151 if (col) hs->c = *col;
1152 ext3d_update (&pos, halfwidth);
1153 }
1154
1155 halfwidth = 0.5 * current_state->helixwidth;
1156 hs = helix_segment_next();
1157 v3_sum_scaled (&(hs->p1), points + 1, halfwidth, axes + 1);
1158 v3_sum_scaled (&(hs->p2), points + 1, -halfwidth, axes + 1);
1159 hs->a = *(axes + 1);
1160 v3_cross_product (&(hs->n), tangents + 1, axes + 1);
1161 v3_normalize (&(hs->n));
1162 if (col) hs->c = *col;
1163 ext3d_update (points + 1, halfwidth);
1164
1165 /* helix main segments */
1166 for (slot = 1; slot < ch->length - 2; slot++) {
1167
1168 hermite_set (points + slot, points + slot + 1,
1169 tangents + slot, tangents + slot + 1);
1170
1171 for (segment = 1; segment < segments; segment++) {
1172 t = ((double) segment) / ((double) segments);
1173 hs = helix_segment_next();
1174 v3_scaled (&(hs->a), 1.0 - t, axes + slot);
1175 v3_add_scaled (&(hs->a), t, axes + slot + 1);
1176 v3_normalize (&(hs->a));
1177 hermite_get (&pos, t);
1178 v3_sum_scaled (&(hs->p1), &pos, halfwidth, &(hs->a));
1179 v3_sum_scaled (&(hs->p2), &pos, -halfwidth, &(hs->a));
1180 hermite_get_tangent (&dir, t);
1181 v3_cross_product (&(hs->n), &dir, &(hs->a));
1182 v3_normalize (&(hs->n));
1183 if (current_state->colourparts && (segment == segments / 2))
1184 col = &(ch->residues[slot+1]->colour);
1185 if (col) hs->c = *col;
1186 ext3d_update (&pos, halfwidth);
1187 }
1188
1189 hs = helix_segment_next();
1190 v3_sum_scaled (&(hs->p1), points + slot + 1, halfwidth, axes + slot + 1);
1191 v3_sum_scaled (&(hs->p2), points + slot + 1, -halfwidth, axes + slot +1);
1192 hs->a = *(axes + slot + 1);
1193 v3_cross_product (&(hs->n), tangents + slot + 1, axes + slot + 1);
1194 v3_normalize (&(hs->n));
1195 if (col) hs->c = *col;
1196 ext3d_update (points + slot + 1, halfwidth);
1197 }
1198 /* helix finish segments */
1199 pos = *(points + ch->length - 2);
1200 hermite_set (points + ch->length - 2, points + ch->length - 1,
1201 tangents + ch->length - 2, tangents + ch->length - 1);
1202
1203 for (segment = 1; segment < segments; segment++) {
1204 t = ((double) segment) / ((double) segments);
1205 hermite_get (&pos, t);
1206 halfwidth = coilradius +
1207 (0.5 * current_state->helixwidth - coilradius) *
1208 0.5 * (cos (ANGLE_PI * t) + 1.0);
1209 hs = helix_segment_next();
1210 v3_sum_scaled (&(hs->p1), &pos, halfwidth, axes + ch->length - 1);
1211 v3_sum_scaled (&(hs->p2), &pos, -halfwidth, axes + ch->length - 1);
1212 hs->a = *(axes + ch->length - 1);
1213 hermite_get_tangent (&dir, t);
1214 v3_cross_product (&(hs->n), &dir, axes + ch->length - 1);
1215 v3_normalize (&(hs->n));
1216 if (current_state->colourparts && (segment == segments / 2))
1217 col = &(ch->residues[ch->length-2]->colour);
1218 if (col) hs->c = *col;
1219 ext3d_update (&pos, halfwidth);
1220 }
1221
1222 hs = helix_segment_next();
1223 v3_sum_scaled (&(hs->p1), points + ch->length - 1,
1224 coilradius, axes + ch->length - 1);
1225 v3_sum_scaled (&(hs->p2), points + ch->length - 1,
1226 -coilradius, axes + ch->length - 1);
1227 hs->a = *(axes + ch->length - 1);
1228 v3_cross_product (&(hs->n),
1229 tangents + ch->length - 1, axes + ch->length - 1);
1230 v3_normalize (&(hs->n));
1231 if (col) hs->c = *col;
1232 ext3d_update (points + ch->length - 1, coilradius);
1233
1234 output_helix();
1235
1236 free (points);
1237 free (axes);
1238 free (tangents);
1239 }
1240
1241 if (first_ch) mol3d_chain_delete (first_ch);
1242
1243 assert (count_residue_selections() == 0);
1244 }
1245
1246
1247 /*------------------------------------------------------------*/
1248 void
line_start(void)1249 line_start (void)
1250 {
1251 assert (dstack_size == 3);
1252
1253 line_segment_init();
1254 line_next();
1255 line_segments->new = TRUE;
1256
1257 assert (dstack_size == 0);
1258 }
1259
1260
1261 /*------------------------------------------------------------*/
1262 void
line_next(void)1263 line_next (void)
1264 {
1265 line_segment *ls;
1266
1267 assert (dstack_size == 3);
1268
1269 ls = line_segment_next();
1270 ls->p.x = dstack[0];
1271 ls->p.y = dstack[1];
1272 ls->p.z = dstack[2];
1273 clear_dstack();
1274 ls->c = current_state->linecolour;
1275
1276 ext3d_update (&(ls->p), 0.0);
1277
1278 assert (dstack_size == 0);
1279 }
1280
1281
1282 /*------------------------------------------------------------*/
1283 static void
apply_labelmask_case(char * str,int len)1284 apply_labelmask_case (char *str, int len)
1285 {
1286 int slot;
1287 int last = (len > current_state->labelmasklength) ?
1288 current_state->labelmasklength : len;
1289 int *mask = current_state->labelmask;
1290
1291 assert (str);
1292 assert (len > 0);
1293
1294 for (slot = 0; slot < last; slot++) {
1295 switch (mask[slot] % 3) {
1296 case 1:
1297 *str = tolower (*str);
1298 break;
1299 case 2:
1300 *str = toupper (*str);
1301 break;
1302 default:
1303 break;
1304 }
1305 str++;
1306 }
1307 }
1308
1309
1310 /*------------------------------------------------------------*/
1311 void
label_atoms(char * label)1312 label_atoms (char *label)
1313 {
1314 at3d *at;
1315 at3d **atoms;
1316 char *orig_str, *copy_str;
1317 int slot, length, rpos, tpos, cpos, apos, atom_count;
1318
1319 assert (count_atom_selections() == 1);
1320 assert (label);
1321 assert (*label);
1322
1323 length = strlen (label) + 4 + 2 - 1 + 2 + 1; /* max possible length */
1324 orig_str = malloc (length * sizeof (char));
1325 copy_str = malloc (length * sizeof (char));
1326
1327 rpos = -1;
1328 tpos = -1;
1329 cpos = -1;
1330 apos = -1;
1331 length = 0;
1332 slot = 0;
1333 while (slot < strlen (label)) {
1334
1335 if (label[slot] == '%') {
1336 slot++;
1337 if (slot >= strlen (label)) goto skip;
1338
1339 switch (label[slot]) {
1340 case 'r':
1341 if (rpos >= 0) goto skip;
1342 rpos = length;
1343 length += 6;
1344 break;
1345 case 't':
1346 if (tpos >= 0) goto skip;
1347 tpos = length;
1348 length += 4;
1349 break;
1350 case 'c':
1351 if (cpos >= 0) goto skip;
1352 cpos = length;
1353 length += 1;
1354 break;
1355 case 'a':
1356 if (apos >= 0) goto skip;
1357 apos = length;
1358 length += 4;
1359 break;
1360 default:
1361 goto skip;
1362 }
1363 } else {
1364 orig_str[length] = label[slot];
1365 length++;
1366 }
1367
1368 slot++;
1369 }
1370
1371 atoms = select_atom_list (&atom_count);
1372
1373 assert (count_atom_selections() == 0);
1374
1375 if (message_mode)
1376 fprintf (stderr, "%i atoms selected for label\n", atom_count);
1377
1378 for (slot = 0; slot < atom_count; slot++) {
1379
1380 at = atoms[slot];
1381
1382 if (rpos >= 0) strncpy (orig_str + rpos, at->res->name, 6);
1383 if (tpos >= 0) strncpy (orig_str + tpos, at->res->type, 4);
1384 if (cpos >= 0) *(orig_str + cpos) = at->res->code;
1385 if (apos >= 0) strncpy (orig_str + apos, at->name, 4);
1386
1387 memcpy (copy_str, orig_str, length * sizeof (char));
1388 apply_labelmask_case (copy_str, length);
1389 str_remove_nulls (copy_str, length);
1390
1391 if (current_state->colourparts) {
1392 output_label (&(at->xyz), copy_str, &(at->colour));
1393 } else {
1394 output_label (&(at->xyz), copy_str, NULL);
1395 }
1396 ext3d_update (&(at->xyz), 0.0);
1397 }
1398
1399 free (copy_str);
1400 free (orig_str);
1401 free (atoms);
1402
1403 return;
1404
1405 skip:
1406 yyerror ("invalid label format string");
1407 }
1408
1409
1410 /*------------------------------------------------------------*/
1411 void
label_position(char * label)1412 label_position (char *label)
1413 {
1414 vector3 pos;
1415
1416 assert (dstack_size == 3);
1417 assert (label);
1418 assert (*label);
1419
1420 pos.x = dstack[0];
1421 pos.y = dstack[1];
1422 pos.z = dstack[2];
1423 clear_dstack();
1424
1425 apply_labelmask_case (label, strlen (label));
1426 output_label (&pos, label, NULL);
1427 ext3d_update (&pos, 0.0);
1428 }
1429
1430
1431 /*------------------------------------------------------------*/
1432 void
object(char * filename)1433 object (char *filename)
1434 {
1435 FILE *file;
1436 char object_code[4];
1437 int close_file, code, alloc, slot, count, number;
1438 vector3 *triplets;
1439 int total = 0;
1440
1441 if (filename) {
1442 file = fopen (filename, "r");
1443 if (file == NULL) {
1444 yyerror ("could not open the object file");
1445 return;
1446 }
1447 close_file = TRUE;
1448 } else {
1449 file = lex_input_file();
1450 close_file = FALSE;
1451 }
1452
1453 alloc = 1024;
1454 triplets = malloc (alloc * sizeof (vector3));
1455
1456 output_start_object();
1457
1458 while (fscanf (file, "%3s", object_code) != EOF) {
1459 if (str_eq (object_code, "P")) {
1460 code = OBJ_POINTS;
1461 } else if (str_eq (object_code, "PC")) {
1462 code = OBJ_POINTS_COLOURS;
1463 } else if (str_eq (object_code, "L")) {
1464 code = OBJ_LINES;
1465 } else if (str_eq (object_code, "LC")) {
1466 code = OBJ_LINES_COLOURS;
1467 } else if (str_eq (object_code, "T")) {
1468 code = OBJ_TRIANGLES;
1469 } else if (str_eq (object_code, "TC")) {
1470 code = OBJ_TRIANGLES_COLOURS;
1471 } else if (str_eq (object_code, "TN")) {
1472 code = OBJ_TRIANGLES_NORMALS;
1473 } else if (str_eq (object_code, "TNC")) {
1474 code = OBJ_TRIANGLES_NORMALS_COLOURS;
1475 } else if (str_eq (object_code, "S")) {
1476 code = OBJ_STRIP;
1477 } else if (str_eq (object_code, "SC")) {
1478 code = OBJ_STRIP_COLOURS;
1479 } else if (str_eq (object_code, "SN")) {
1480 code = OBJ_STRIP_NORMALS;
1481 } else if (str_eq (object_code, "SNC")) {
1482 code = OBJ_STRIP_NORMALS_COLOURS;
1483 } else if (str_eq (object_code, "Q")) {
1484 goto finish;
1485 } else {
1486 goto format_error;
1487 }
1488
1489 if (fscanf (file, "%i", &count) != 1) goto format_error;
1490
1491 switch (code) {
1492 case OBJ_POINTS:
1493 if (count < 1) goto format_error;
1494 number = 1;
1495 break;
1496 case OBJ_POINTS_COLOURS:
1497 if (count < 1) goto format_error;
1498 number = 2;
1499 break;
1500 case OBJ_LINES:
1501 if (count < 2) goto format_error;
1502 number = 1;
1503 break;
1504 case OBJ_LINES_COLOURS:
1505 if (count < 2) goto format_error;
1506 number = 2;
1507 break;
1508 case OBJ_TRIANGLES:
1509 if ((count < 3) || (count % 3 != 0)) goto format_error;
1510 number = 1;
1511 break;
1512 case OBJ_TRIANGLES_COLOURS:
1513 case OBJ_TRIANGLES_NORMALS:
1514 if ((count < 3) || (count % 3 != 0)) goto format_error;
1515 number = 2;
1516 break;
1517 case OBJ_TRIANGLES_NORMALS_COLOURS:
1518 if ((count < 3) || (count % 3 != 0)) goto format_error;
1519 number = 3;
1520 break;
1521 case OBJ_STRIP:
1522 if (count < 3) goto format_error;
1523 number = 1;
1524 break;
1525 case OBJ_STRIP_COLOURS:
1526 case OBJ_STRIP_NORMALS:
1527 if (count < 3) goto format_error;
1528 number = 2;
1529 break;
1530 case OBJ_STRIP_NORMALS_COLOURS:
1531 if (count < 3) goto format_error;
1532 number = 3;
1533 break;
1534 default:
1535 goto format_error;
1536 }
1537
1538 count *= number;
1539
1540 if (count > alloc) {
1541 alloc = count;
1542 triplets = realloc (triplets, count * sizeof (vector3));
1543 }
1544
1545 for (slot = 0; slot < count; slot++) {
1546 if (fscanf (file, "%lg %lg %lg",
1547 &(triplets[slot].x),
1548 &(triplets[slot].y),
1549 &(triplets[slot].z)) != 3) goto format_error;
1550 }
1551
1552 if (current_state->objecttransform) {
1553 for (slot = 0; slot < count; slot += number) { /* coordinate xform */
1554 matrix3_transform (triplets + slot, xform);
1555 }
1556 if (code == OBJ_TRIANGLES_NORMALS || /* normals rotate, not xform */
1557 code == OBJ_TRIANGLES_NORMALS_COLOURS ||
1558 code == OBJ_STRIP_NORMALS ||
1559 code == OBJ_STRIP_NORMALS_COLOURS) {
1560 for (slot = 1; slot < count; slot += number) {
1561 matrix3_rotate (triplets + slot, xform);
1562 }
1563 }
1564 }
1565
1566 for (slot = 0; slot < count; slot += number) {
1567 ext3d_update (triplets + slot, 0.0);
1568 }
1569
1570 output_object (code, triplets, count);
1571 total += count;
1572 }
1573
1574 finish:
1575 free (triplets);
1576 output_finish_object();
1577 if (close_file) fclose (file);
1578 if (message_mode)
1579 fprintf (stderr, "%i data triplets read from object file\n", total);
1580 return;
1581
1582 format_error:
1583 free (triplets);
1584 output_finish_object();
1585 if (close_file) fclose (file);
1586 yyerror ("invalid format or content in object file");
1587 }
1588
1589
1590 /*------------------------------------------------------------*/
1591 void
strand(void)1592 strand (void)
1593 {
1594 mol3d_chain *first_ch, *ch;
1595 vector3 *points, *normals, *smoothed;
1596 vector3 dir, dir1, dir2, pos, side, vec1, vec2, normal;
1597 colour *col;
1598 strand_segment *ss;
1599 int slot, segment, segments;
1600 double thickness, radius, t;
1601 double width = current_state->strandwidth / 2.0;
1602
1603 assert (count_residue_selections() == 1);
1604
1605 if (current_state->colourparts) { /* strand is less curved; */
1606 segments = current_state->segments / 2; /* colourparts requires */
1607 if (segments % 2 != 0) segments++; /* even number of segments */
1608 } else {
1609 segments = current_state->segments / 2 + 1;
1610 }
1611
1612 if (current_state->strandthickness > current_state->strandwidth) {
1613 radius = current_state->strandthickness / 2.0;
1614 } else {
1615 radius = current_state->strandwidth / 2.0;
1616 }
1617
1618 if (current_state->strandthickness < 0.01) {
1619 thickness = 0.0;
1620 } else {
1621 thickness = current_state->strandthickness / 2.0;
1622 }
1623
1624 first_ch = get_peptide_chains();
1625
1626 for (ch = first_ch; ch; ch = ch->next) {
1627 if (ch->length < 3) continue;
1628
1629 msg_chain ("strand", ch);
1630
1631 points = get_atom_positions (ch);
1632 /* normals for the strand */
1633 normals = malloc (ch->length * sizeof (vector3));
1634 for (slot = 1; slot < ch->length - 1; slot++) {
1635 v3_middle (&pos, points + slot - 1, points + slot + 1);
1636 v3_difference (normals + slot, points + slot, &pos);
1637 v3_normalize (normals + slot);
1638 }
1639 /* just copy the normals for the terminii */
1640 *(normals) = *(normals + 1);
1641 *(normals + ch->length - 1) = *(normals + ch->length - 2);
1642
1643 /* smooth CA positions */
1644 priestle_smoothing (points, ch->length, current_state->smoothsteps);
1645
1646 /* normals must point the same way */
1647 for (slot = 0; slot < ch->length - 1; slot++) {
1648 if (v3_dot_product (normals + slot, normals + slot + 1) < 0.0) {
1649 v3_reverse (normals + slot + 1);
1650 }
1651 }
1652 /* smooth normals, one iteration */
1653 smoothed = malloc (ch->length * sizeof (vector3));
1654 for (slot = 1; slot < ch->length - 1; slot++) {
1655 v3_sum (&dir, normals + slot - 1, normals + slot);
1656 v3_add (&dir, normals + slot + 1);
1657 v3_normalize (&dir);
1658 *(smoothed + slot) = dir;
1659 }
1660 for (slot = 1; slot < ch->length - 1; slot++) {
1661 *(normals + slot) = *(smoothed + slot);
1662 }
1663 free (smoothed);
1664
1665 /* normals exactly perpendicular to strand */
1666 v3_difference (&dir, points + 1, points);
1667 v3_cross_product (&side, &dir, normals);
1668 v3_cross_product (normals, &side, &dir);
1669 v3_normalize (normals);
1670 for (slot = 1; slot < ch->length - 1; slot++) {
1671 v3_difference (&dir, points + slot + 1, points + slot - 1);
1672 v3_cross_product (&side, &dir, normals + slot);
1673 v3_cross_product (normals + slot, &side, &dir);
1674 v3_normalize (normals + slot);
1675 }
1676 v3_difference (&dir, points + ch->length - 1, points + ch->length - 2);
1677 v3_cross_product (&side, &dir, normals + ch->length - 1);
1678 v3_cross_product (normals + ch->length - 1, &side, &dir);
1679 v3_normalize (normals + ch->length - 1);
1680
1681 if (current_state->colourparts) {
1682 col = &(ch->residues[0]->colour);
1683 } else {
1684 col = NULL;
1685 }
1686
1687 strand_segment_init();
1688 /* strand body */
1689 v3_difference (&dir2, points + 1, points);
1690 v3_normalize (&dir2);
1691 v3_difference (&vec2, points + 1, points);
1692 v3_scale (&vec2, STRAND_HERMITE_FACTOR);
1693
1694 for (slot = 0; slot < ch->length - 2; slot++) {
1695
1696 dir1 = dir2; /* direction vectors */
1697 v3_difference (&dir2, points + slot + 2, points + slot);
1698 v3_normalize (&dir2);
1699 vec1 = vec2; /* Hermite spline vectors */
1700 v3_difference (&vec2, points + slot + 2, points + slot);
1701 v3_scale (&vec2, STRAND_HERMITE_FACTOR);
1702 hermite_set (points + slot, points + slot + 1, &vec1, &vec2);
1703
1704 for (segment = 0; segment < segments; segment++) {
1705 t = ((double) segment) / ((double) segments);
1706 hermite_get (&pos, t);
1707 v3_scaled (&dir, 1.0 - t, &dir1);
1708 v3_add_scaled (&dir, t, &dir2);
1709 v3_scaled (&normal, 1.0 - t, normals + slot);
1710 v3_add_scaled (&normal, t, normals + slot + 1);
1711 v3_cross_product (&side, &normal, &dir);
1712 v3_normalize (&side);
1713 if (current_state->colourparts &&
1714 (segment == segments / 2)) col = &(ch->residues[slot+1]->colour);
1715 ss = strand_segment_next();
1716 v3_sum_scaled (&(ss->p1), &pos, width, &side);
1717 v3_add_scaled (&(ss->p1), thickness, &normal);
1718 v3_sum_scaled (&(ss->p2), &pos, width, &side);
1719 v3_add_scaled (&(ss->p2), -thickness, &normal);
1720 v3_sum_scaled (&(ss->p3), &pos, -width, &side);
1721 v3_add_scaled (&(ss->p3), -thickness, &normal);
1722 v3_sum_scaled (&(ss->p4), &pos, -width, &side);
1723 v3_add_scaled (&(ss->p4), thickness, &normal);
1724 ss->n1 = normal;
1725 ss->n2 = side;
1726 v3_scaled (&(ss->n3), -1.0, &normal);
1727 v3_scaled (&(ss->n4), -1.0, &side);
1728 if (col) ss->c = *col;
1729 ext3d_update (&pos, radius);
1730 }
1731 }
1732 /* strand body, last segment */
1733 dir = dir2;
1734 normal = *(normals + ch->length - 2);
1735 v3_cross_product (&side, &normal, &dir);
1736 v3_normalize (&side);
1737 pos = *(points + ch->length - 2);
1738 ss = strand_segment_next();
1739 v3_sum_scaled (&(ss->p1), &pos, width, &side);
1740 v3_add_scaled (&(ss->p1), thickness, &normal);
1741 v3_sum_scaled (&(ss->p2), &pos, width, &side);
1742 v3_add_scaled (&(ss->p2), -thickness, &normal);
1743 v3_sum_scaled (&(ss->p3), &pos, -width, &side);
1744 v3_add_scaled (&(ss->p3), -thickness, &normal);
1745 v3_sum_scaled (&(ss->p4), &pos, -width, &side);
1746 v3_add_scaled (&(ss->p4), thickness, &normal);
1747 ss->n1 = normal;
1748 ss->n2 = side;
1749 v3_scaled (&(ss->n3), -1.0, &normal);
1750 v3_scaled (&(ss->n4), -1.0, &side);
1751 if (col) ss->c = *col;
1752 ext3d_update (&pos, radius);
1753 /* strand arrow head: 3 segments */
1754 ss = strand_segment_next();
1755 v3_sum_scaled (&(ss->p1), &pos, 1.5 * width, &side);
1756 v3_add_scaled (&(ss->p1), thickness, &normal);
1757 v3_sum_scaled (&(ss->p2), &pos, 1.5 * width, &side);
1758 v3_add_scaled (&(ss->p2), -thickness, &normal);
1759 v3_sum_scaled (&(ss->p3), &pos, -1.5 * width, &side);
1760 v3_add_scaled (&(ss->p3), -thickness, &normal);
1761 v3_sum_scaled (&(ss->p4), &pos, -1.5 * width, &side);
1762 v3_add_scaled (&(ss->p4), thickness, &normal);
1763 ss->n1 = normal;
1764 ss->n2 = side;
1765 v3_scaled (&(ss->n3), -1.0, &normal);
1766 v3_scaled (&(ss->n4), -1.0, &side);
1767 if (col) ss->c = *col;
1768 ext3d_update (&pos, 1.5 * radius);
1769
1770 dir1 = dir2;
1771 v3_difference (&dir2, points + ch->length - 1, points + ch->length - 2);
1772 v3_normalize (&dir1);
1773 v3_middle (&dir, &dir1, &dir2);
1774 v3_middle (&pos, points + ch->length - 2, points + ch->length - 1);
1775 v3_middle (&normal, normals + ch->length - 2, normals + ch->length - 1);
1776 v3_normalize (&normal);
1777 v3_cross_product (&side, &normal, &dir);
1778 v3_normalize (&side);
1779 ss = strand_segment_next();
1780 v3_sum_scaled (&(ss->p1), &pos, 0.75 * width, &side);
1781 v3_add_scaled (&(ss->p1), thickness, &normal);
1782 v3_sum_scaled (&(ss->p2), &pos, 0.75 * width, &side);
1783 v3_add_scaled (&(ss->p2), -thickness, &normal);
1784 v3_sum_scaled (&(ss->p3), &pos, -0.75 * width, &side);
1785 v3_add_scaled (&(ss->p3), -thickness, &normal);
1786 v3_sum_scaled (&(ss->p4), &pos, -0.75 * width, &side);
1787 v3_add_scaled (&(ss->p4), thickness, &normal);
1788 ss->n1 = normal;
1789 ss->n2 = side;
1790 v3_scaled (&(ss->n3), -1.0, &normal);
1791 v3_scaled (&(ss->n4), -1.0, &side);
1792 if (current_state->colourparts) ss->c = ch->residues[ch->length-1]->colour;
1793 ext3d_update (&pos, 0.75 * radius);
1794
1795 pos = *(points + ch->length - 1);
1796 normal = *(normals + ch->length - 1);
1797 ss = strand_segment_next();
1798 v3_sum_scaled (&(ss->p1), &pos, thickness, &normal);
1799 v3_sum_scaled (&(ss->p2), &pos, -thickness, &normal);
1800 ss->n1 = normal;
1801 v3_scaled (&(ss->n2), -1.0, &normal);
1802 ext3d_update (points + ch->length - 1, thickness);
1803
1804 output_strand();
1805
1806 free (normals);
1807 free (points);
1808 }
1809
1810 if (first_ch) mol3d_chain_delete (first_ch);
1811
1812 assert (count_residue_selections() == 0);
1813 }
1814
1815
1816 /*------------------------------------------------------------*/
1817 void
trace(void)1818 trace (void)
1819 {
1820 int slot;
1821 mol3d_chain *first_ch, *ch;
1822 line_segment *ls;
1823
1824 assert (count_residue_selections() == 1);
1825
1826 first_ch = get_peptide_chains();
1827
1828 line_segment_init();
1829
1830 for (ch = first_ch; ch; ch = ch->next) {
1831 if (ch->length < 2) continue;
1832
1833 ls = line_segment_next();
1834 ls->new = TRUE;
1835 ls->p = ch->atoms[0]->xyz;
1836 ext3d_update (&(ls->p), 0.0);
1837
1838 if (current_state->colourparts) {
1839 vector3 middle;
1840
1841 ls->c = ch->residues[0]->colour;
1842 for (slot = 1; slot < ch->length; slot++) {
1843 if (colour_unequal (&(ch->residues[slot-1]->colour),
1844 &(ch->residues[slot]->colour))) {
1845 ls = line_segment_next();
1846 v3_middle (&middle, &(ch->atoms[slot-1]->xyz),
1847 &(ch->atoms[slot]->xyz));
1848 ls->p = middle;
1849 ls->c = ch->residues[slot-1]->colour;
1850 ls = line_segment_next();
1851 ls->new = TRUE;
1852 ls->p = middle;
1853 ls->c = ch->residues[slot]->colour;
1854 }
1855 ls = line_segment_next();
1856 ls->p = ch->atoms[slot]->xyz;
1857 ls->c = ch->residues[slot]->colour;
1858 ext3d_update (&(ls->p), 0.0);
1859 }
1860
1861 } else {
1862 ls->c = current_state->linecolour;
1863 for (slot = 1; slot < ch->length; slot++) {
1864 ls = line_segment_next();
1865 ls->p = ch->atoms[slot]->xyz;
1866 ext3d_update (&(ls->p), 0.0);
1867 }
1868 }
1869
1870 if (output_pickable) {
1871 for (slot = 0; slot < ch->length; slot++) {
1872 output_pickable (ch->atoms[slot]);
1873 }
1874 }
1875
1876 msg_chain ("trace", ch);
1877 }
1878
1879 output_line (TRUE);
1880
1881 if (first_ch) mol3d_chain_delete (first_ch);
1882
1883 assert (count_residue_selections() == 0);
1884 }
1885