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 (&center, &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