001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021/**
022 *
023 */
024package org.biojava.nbio.structure.symmetry.jmolScript;
025
026import org.biojava.nbio.structure.symmetry.axis.RotationAxisAligner;
027import org.biojava.nbio.structure.symmetry.core.Rotation;
028import org.biojava.nbio.structure.symmetry.core.RotationGroup;
029import org.biojava.nbio.structure.symmetry.core.QuatSymmetrySubunits;
030import org.biojava.nbio.structure.symmetry.geometry.Polyhedron;
031import org.jcolorbrewer.ColorBrewer;
032
033import javax.vecmath.*;
034
035import java.awt.*;
036import java.util.*;
037import java.util.List;
038
039/**
040 * @author Peter
041 *
042 */
043public abstract class JmolSymmetryScriptGeneratorPointGroup extends JmolSymmetryScriptGenerator {
044        private static String N_FOLD_AXIS_COLOR = "red";
045        private static String TWO_FOLD_AXIS_COLOR = "deepskyblue";
046        private static String THREE_FOLD_AXIS_COLOR = "lime";
047        private static double AXIS_SCALE_FACTOR = 1.2;
048
049        private RotationAxisAligner rotationAxisAligner = null;
050        private RotationGroup rotationGroup = null;
051        private Polyhedron polyhedron = null;
052        private String name = "";
053        private String defaultColoring = "";
054        private boolean onTheFly = true;
055
056        public JmolSymmetryScriptGeneratorPointGroup(RotationAxisAligner rotationAxisAligner, String name) {
057                this.rotationAxisAligner = rotationAxisAligner;
058                this.rotationGroup = this.rotationAxisAligner.getRotationGroup();
059                this.name = name;
060        }
061
062        /* (non-Javadoc)
063         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getZoom()
064         */
065        @Override
066        abstract public int getZoom();
067
068        @Override
069        public void setOnTheFly(boolean onTheFly) {
070                this.onTheFly = onTheFly;
071        }
072
073        /* (non-Javadoc)
074         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getDefaultOrientation()
075         */
076        @Override
077        public String getDefaultOrientation() {
078                StringBuilder s = new StringBuilder();
079                s.append(setCentroid());
080
081                // calculate  orientation
082                Quat4d q = new Quat4d();
083                q.set(polyhedron.getViewMatrix(0));
084                q.normalize();
085                Quat4d r = new Quat4d();
086                r.set(rotationAxisAligner.getRotationMatrix());
087                r.normalize();
088                q.mul(r);
089                q.normalize();
090
091                // set orientation
092                s.append("moveto 0 quaternion{");
093                s.append(jMolFloat(q.x));
094                s.append(",");
095                s.append(jMolFloat(q.y));
096                s.append(",");
097                s.append(jMolFloat(q.z));
098                s.append(",");
099                s.append(jMolFloat(q.w));
100                s.append("};");
101                return s.toString();
102        }
103
104        /* (non-Javadoc)
105         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getOrientationCount()
106         */
107        @Override
108        public int getOrientationCount() {
109                return polyhedron.getViewCount();
110        }
111
112        /* (non-Javadoc)
113         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getOrientation(int)
114         */
115        @Override
116        public String getOrientation(int index) {
117                StringBuilder s = new StringBuilder();
118                s.append(setCentroid());
119
120                // calculate  orientation
121                Quat4d q = new Quat4d();
122                q.set(polyhedron.getViewMatrix(index));
123                q.normalize();
124                Quat4d r = new Quat4d();
125                r.set(rotationAxisAligner.getRotationMatrix());
126                r.normalize();
127                q.mul(r);
128                q.normalize();
129
130                // set orientation
131                s.append("moveto 4 quaternion{");
132                s.append(jMolFloat(q.x));
133                s.append(",");
134                s.append(jMolFloat(q.y));
135                s.append(",");
136                s.append(jMolFloat(q.z));
137                s.append(",");
138                s.append(jMolFloat(q.w));
139                s.append("}");
140                s.append(";");
141                return s.toString();
142        }
143
144        /* (non-Javadoc)
145         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getOrientationWithZoom(int)
146         */
147        @Override
148        public String getOrientationWithZoom(int index) {
149                StringBuilder s = new StringBuilder();
150                s.append(getOrientation(index));
151                s.insert(s.length()-1, getZoom());
152                return s.toString();
153
154        }
155        /* (non-Javadoc)
156         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getOrientationName(int)
157         */
158        @Override
159        public String getOrientationName(int index) {
160                return polyhedron.getViewName(index);
161        }
162
163        /* (non-Javadoc)
164         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getTransformation()
165         */
166        @Override
167        public Matrix4d getTransformation() {
168                return rotationAxisAligner.getTransformation();
169        }
170
171        @Override
172        public void setDefaultColoring(String colorScript) {
173                this.defaultColoring = colorScript;
174        }
175
176        /* (non-Javadoc)
177         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#drawPolyhedron()
178         */
179        @Override
180        public String drawPolyhedron() {
181                StringBuilder s = new StringBuilder();
182
183                Point3d[] vertices = getPolyhedronVertices();
184
185                int index = 0;
186                double width = getMaxExtension()*0.015;
187
188                for (int[] lineLoop: polyhedron.getLineLoops()) {
189                        s.append("draw polyhedron");
190                        s.append(name);
191                        s.append(index++);
192                        s.append(" line");
193                        for (int i: lineLoop) {
194                                s.append(getJmolPoint(vertices[i]));
195                        }
196                        s.append("width ");
197                        s.append(fDot2(width));
198                        s.append(" color");
199                        Color4f c = getPolyhedronColor();
200                        s.append(getJmolColor(c));
201                        s.append(" off;");
202                }
203
204                return s.toString();
205        }
206
207        /* (non-Javadoc)
208         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#hidePolyhedron()
209         */
210        @Override
211        public String hidePolyhedron() {
212                return "draw polyhedron* off;";
213        }
214
215        /* (non-Javadoc)
216         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#showPolyhedron()
217         */
218        @Override
219        public String showPolyhedron() {
220                return "draw polyhedron* on;";
221        }
222
223        /* (non-Javadoc)
224         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#drawAxes()
225         */
226        @Override
227        public String drawAxes() {
228                if (rotationGroup.getPointGroup().equals("C1")) {
229                        return drawInertiaAxes();
230                } else {
231                        return drawSymmetryAxes();
232                }
233        }
234
235        /* (non-Javadoc)
236         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#hideAxes()
237         */
238        @Override
239        public String hideAxes() {
240                return "draw axes* off;";
241        }
242
243        /* (non-Javadoc)
244         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#showAxes()
245         */
246        @Override
247        public String showAxes() {
248                return "draw axes* on;";
249        }
250
251        /* (non-Javadoc)
252         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#playOrientations()
253         */
254        @Override
255        public String playOrientations() {
256                StringBuilder s = new StringBuilder();
257
258                // draw point group
259
260                if ( rotationGroup.getPointGroup().equals("C1")) {
261                        s.append(drawFooter("Asymmetric", "white"));
262                } else {
263                        s.append(drawFooter("Point group " + rotationGroup.getPointGroup(), "white"));
264                }
265
266                // draw polygon
267                s.append(drawPolyhedron()); // draw invisibly
268                s.append(showPolyhedron());
269
270                // draw axes
271                s.append(drawAxes());
272                s.append(showAxes());
273
274                // loop over all orientations with 4 sec. delay
275                for (int i = 0; i < getOrientationCount(); i++) {
276                        s.append(deleteHeader());
277                        s.append(getOrientationWithZoom(i));
278                        s.append(drawHeader(polyhedron.getViewName(i), "white"));
279                        s.append("delay 4;");
280                }
281
282                // go back to first orientation
283                s.append(deleteHeader());
284                s.append(getOrientationWithZoom(0));
285                s.append(drawHeader(polyhedron.getViewName(0), "white"));
286
287                return s.toString();
288        }
289
290        /* (non-Javadoc)
291         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#colorBySubunit()
292         */
293        @Override
294        public String colorBySubunit() {
295                QuatSymmetrySubunits subunits = rotationAxisAligner.getSubunits();
296                List<Integer> modelNumbers = subunits.getModelNumbers();
297                List<String> chainIds = subunits.getChainIds();
298                List<List<Integer>> orbits = rotationAxisAligner.getOrbits();
299                int fold = rotationGroup.getRotation(0).getFold();
300
301                Color[] col = null;
302                Color4f[] colors = null;
303                if (fold > 1) {
304                        col = ColorBrewer.Spectral.getColorPalette(2*fold);
305                        colors = ColorConverter.convertColor4f(col);
306                } else {
307                        col = ColorBrewer.Spectral.getColorPalette(orbits.size());
308                        colors = ColorConverter.convertColor4f(col);
309                }
310                int half = colors.length/2;
311                for (int i = 0; i < half; i++) {
312                        if (i % 2 != 0) {
313                           Color4f temp = colors[i];
314                           colors[i] = colors[half+i];
315                           colors[half+i] = temp;
316                        }
317                }
318                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
319
320                for (int i = 0; i < orbits.size(); i++) {
321                        for (int j = 0; j < fold; j++) {
322                                // assign alternating color sets to adjacent orbits
323                                int colorIndex = i;
324                                if (fold > 1) {
325                                        if (i % 2 == 0) {
326                                                colorIndex = j;
327                                        } else {
328                                                colorIndex = fold + j;
329                                        }
330                                }
331                                int subunit = orbits.get(i).get(j);
332                                Color4f c = colors[colorIndex];
333                                List<String> ids = colorMap.get(c);
334                                if (ids == null) {
335                                        ids = new ArrayList<String>();
336                                        colorMap.put(c,  ids);
337                                }
338                                String id = getChainSpecification(modelNumbers, chainIds, subunit);
339                                ids.add(id);
340                        }
341                }
342                return defaultColoring + getJmolColorScript(colorMap) + getJmolLigandScript();
343        }
344
345
346        /* (non-Javadoc)
347         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#colorBySequenceCluster()
348         */
349        @Override
350        public String colorBySequenceCluster() {
351                QuatSymmetrySubunits subunits = rotationAxisAligner.getSubunits();
352                int n = subunits.getSubunitCount();
353                List<Integer> modelNumbers = subunits.getModelNumbers();
354                List<String> chainIds = subunits.getChainIds();
355                List<Integer> seqClusterIds = subunits.getClusterIds();
356                int clusters = Collections.max(seqClusterIds) + 1;
357                Color[] col = ColorBrewer.BrBG.getColorPalette(clusters);
358                Color4f[] colors = ColorConverter.convertColor4f(col);
359                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
360
361                for (int i = 0; i < n; i++) {
362                        Color4f c = colors[seqClusterIds.get(i)];
363                        List<String> ids = colorMap.get(c);
364                        if (ids == null) {
365                                ids = new ArrayList<String>();
366                                colorMap.put(c,  ids);
367                        }
368                        String id = getChainSpecification(modelNumbers, chainIds, i);
369                        ids.add(id);
370
371                }
372
373                return defaultColoring + getJmolColorScript(colorMap) + getJmolLigandScript();
374        }
375
376        /* (non-Javadoc)
377         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#colorBySymmetry()
378         */
379        @Override
380        public String colorBySymmetry() {
381                // TODO needs some refactoring
382                String pointGroup = rotationGroup.getPointGroup();
383                QuatSymmetrySubunits subunits = rotationAxisAligner.getSubunits();
384                List<Integer> modelNumbers = subunits.getModelNumbers();
385                List<String> chainIds = subunits.getChainIds();
386                List<List<Integer>> orbits = rotationAxisAligner.getOrbits();
387
388                int n = subunits.getSubunitCount();
389                int fold = rotationGroup.getRotation(0).getFold();
390
391                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
392
393                // Simple Cn symmetry
394                if (pointGroup.startsWith("C") && n == fold) {
395                        colorMap = getCnColorMap();
396                        // complex cases
397                } else if ((pointGroup.startsWith("D") && orbits.size() > 2) ||
398                                pointGroup.equals("T")|| pointGroup.equals("O") || pointGroup.equals("I")) {
399                        int nColor = 0;
400                        if (orbits.size() % 2 == 0) {
401                                nColor = orbits.size()/2;
402                        } else {
403                                nColor = (orbits.size() + 1)/2;
404                        }
405                        Color4f[] colors = getSymmetryColors(nColor);
406
407                        for (int i = 0; i < orbits.size(); i++) {
408                                int colorIndex = i;
409
410                                // reverse colors once the center of the structure has been reached
411                                if (i >= nColor) {
412                                        colorIndex = orbits.size() - 1 - i;
413                                }
414                                Color4f c = colors[colorIndex];
415                                List<String> ids = colorMap.get(c);
416                                if (ids == null) {
417                                        ids = new ArrayList<String>();
418                                        colorMap.put(c,  ids);
419                                }
420                                for (int subunit: orbits.get(i)) {
421                                        String id = getChainSpecification(modelNumbers, chainIds, subunit);
422                                        ids.add(id);
423                                }
424                        }
425
426                        // Simple Dn symmetry
427                } else {
428                        Color4f[] colors = getSymmetryColors(orbits.size());
429
430                        for (int i = 0; i < orbits.size(); i++) {
431                                Color4f c = new Color4f(colors[i]);
432                                List<String> ids = colorMap.get(c);
433                                if (ids == null) {
434                                        ids = new ArrayList<String>();
435                                        colorMap.put(c,  ids);
436                                }
437                                List<Integer> orbit = orbits.get(i);
438                                for (int j = 0; j < orbit.size(); j++) {
439                                        String id = getChainSpecification(modelNumbers, chainIds, orbit.get(j));
440                                        ids.add(id);
441                                }
442                        }
443                }
444
445                return defaultColoring + getJmolColorScript(colorMap) + getJmolLigandScript();
446        }
447
448        // --- protected methods ---
449        /**
450         * Returns the maximum extension (length) of structure
451         * @return
452         */
453        protected double getMaxExtension() {
454                Vector3d dimension = rotationAxisAligner.getDimension();
455                double maxExtension = Math.max(dimension.x, dimension.y);
456                maxExtension = Math.max(maxExtension, dimension.z);
457                return maxExtension;
458        }
459
460        /**
461         * Returns the mean extension (length) of structure
462         * @return
463         */
464        protected double getMeanExtension() {
465                Vector3d dimension = rotationAxisAligner.getDimension();
466                return (dimension.x+dimension.y+dimension.z)/3;
467        }
468
469        /**
470         * @return the axisTransformation
471         */
472        protected RotationAxisAligner getAxisTransformation() {
473                return rotationAxisAligner;
474        }
475
476        /**
477         * @param axisTransformation the axisTransformation to set
478         */
479        protected void setAxisTransformation(RotationAxisAligner axisTransformation) {
480                this.rotationAxisAligner = axisTransformation;
481        }
482
483        /**
484         * @return the rotationGroup
485         */
486        protected RotationGroup getRotationGroup() {
487                return rotationGroup;
488        }
489
490        /**
491         * @param rotationGroup the rotationGroup to set
492         */
493        protected void setRotationGroup(RotationGroup rotationGroup) {
494                this.rotationGroup = rotationGroup;
495        }
496
497        /**
498         * @return the polyhedron
499         */
500        protected Polyhedron getPolyhedron() {
501                return polyhedron;
502        }
503
504        /**
505         * @param polyhedron the polyhedron to set
506         */
507        protected void setPolyhedron(Polyhedron polyhedron) {
508                this.polyhedron = polyhedron;
509        }
510
511//  --- private methods ---
512
513        private String getChainSpecification(List<Integer> modelNumbers, List<String> chainIds, int subunit) {
514                if (onTheFly) {
515                        if (Collections.max(modelNumbers) > 1) {
516                                return chainIds.get(subunit) + "&symop=" + (modelNumbers.get(subunit)+1);
517                        } else {
518                                // if there is only a single symop, Jmol does not accept the symop syntax
519                                return chainIds.get(subunit);
520                        }
521                } else {
522                        return chainIds.get(subunit) + "/" + (modelNumbers.get(subunit)+1);
523                }
524        }
525
526        private Map<Color4f, List<String>> getCnColorMap() {
527                QuatSymmetrySubunits subunits = rotationAxisAligner.getSubunits();
528                List<Integer> modelNumbers = subunits.getModelNumbers();
529                List<String> chainIds = subunits.getChainIds();
530                List<List<Integer>> orbits = rotationAxisAligner.getOrbits();
531
532                int fold = rotationGroup.getRotation(0).getFold();
533
534                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
535                Color4f[] colors = getSymmetryColors(fold);
536
537                for (List<Integer> orbit: orbits) {
538                        for (int i = 0; i < fold; i++) {
539                                int subunit = orbit.get(i);
540                                String id = null;
541                                id = getChainSpecification(modelNumbers, chainIds, subunit);
542                                Color4f c = colors[i];
543                                List<String> ids = colorMap.get(c);
544                                if (ids == null) {
545                                        ids = new ArrayList<String>();
546                                        colorMap.put(c, ids);
547                                }
548                                ids.add(id);
549                        }
550                }
551
552                return colorMap;
553        }
554
555        private Point3d[] getPolyhedronVertices() {
556                Point3d[] vertices = polyhedron.getVertices();
557                Matrix4d reverseTransformation = rotationAxisAligner.getGeometicCenterTransformation();
558                for (int i = 0; i < vertices.length; i++) {
559                        reverseTransformation.transform(vertices[i]);
560                }
561                return vertices;
562        }
563
564        /**
565         * Return a color that is complementary to the symmetry color
566         * @return
567         */
568        private Color4f getPolyhedronColor() {
569                Color4f[] colors = getSymmetryColors(5);
570                Color4f strongestColor = colors[4];
571                Color4f complement = new Color4f(Color.WHITE);
572                complement.sub(strongestColor);
573                return complement;
574        }
575
576        /**
577         * Returns a unique color palette based on point group
578         * @param nColors
579         * @return
580         */
581        private Color4f[] getSymmetryColors(int nColors) {
582                String pointGroup = rotationGroup.getPointGroup();
583                Color[] col = null;
584                Color4f[] colors = null;
585                if (pointGroup.equals("C1")) {
586                        col = ColorBrewer.Greys.getColorPalette(nColors);
587                        colors = ColorConverter.convertColor4f(col);
588                } else if (pointGroup.startsWith("C")) {
589                        col = ColorBrewer.YlGnBu.getColorPalette(nColors);
590                        colors = ColorConverter.convertColor4f(col);
591                } else if (pointGroup.startsWith("D")) {
592                        col = ColorBrewer.YlOrRd.getColorPalette(nColors);
593                        colors = ColorConverter.convertColor4f(col);
594                } else if (pointGroup.equals("T")) {
595                        col = ColorBrewer.Greens.getColorPalette(nColors);
596                        colors = ColorConverter.convertColor4f(col);
597                } else if (pointGroup.equals("O")) {
598                        col = ColorBrewer.Blues.getColorPalette(nColors);
599                        colors = ColorConverter.convertColor4f(col);
600                } else if (pointGroup.equals("I")) {
601                        col = ColorBrewer.BuPu.getColorPalette(nColors);
602                        colors = ColorConverter.convertColor4f(col);
603                } else {
604                        col = ColorBrewer.Greys.getColorPalette(nColors);
605                        colors = ColorConverter.convertColor4f(col);
606                }
607                System.arraycopy(colors, 0, colors, 0, nColors);
608                return colors;
609        }
610
611        private String drawInertiaAxes() {
612                StringBuilder s = new StringBuilder();
613                Point3d centroid = rotationAxisAligner.getGeometricCenter();
614                Vector3d[] axes = rotationAxisAligner.getPrincipalAxesOfInertia();
615
616                for (int i = 0; i < axes.length; i++) {
617                        s.append("draw axesInertia");
618                        s.append(name);
619                        s.append(i);
620                        s.append(" ");
621                        s.append("line");
622                        Point3d v1 = new Point3d(axes[i]);
623                        if (i == 0) {
624                                v1.scale(AXIS_SCALE_FACTOR*rotationAxisAligner.getDimension().y);
625                        } else if (i == 1) {
626                                v1.scale(AXIS_SCALE_FACTOR*rotationAxisAligner.getDimension().x);
627                        } else if (i == 2) {
628                                v1.scale(AXIS_SCALE_FACTOR*rotationAxisAligner.getDimension().z);
629                        }
630                        Point3d v2 = new Point3d(v1);
631                        v2.negate();
632                        v1.add(centroid);
633                        v2.add(centroid);
634                        s.append(getJmolPoint(v1));
635                        s.append(getJmolPoint(v2));
636                        s.append("width 0.5 ");
637                        s.append(" color white");
638                        s.append(" off;");
639                }
640                return s.toString();
641        };
642
643        private String drawSymmetryAxes() {
644                StringBuilder s = new StringBuilder();
645
646                int n = rotationGroup.getOrder();
647                if (n == 0) {
648                        return s.toString();
649                }
650
651                float diameter = 0.5f;
652                double radius = 0;
653                String color = "";
654
655                List<Rotation> axes = getUniqueAxes();
656
657                int i = 0;
658                for (Rotation r: axes) {
659                        if (rotationGroup.getPointGroup().startsWith("C") || (rotationGroup.getPointGroup().startsWith("D") && r.getDirection() == 0)) {
660                                radius =  rotationAxisAligner.getDimension().z; // principal axis uses z-dimension
661                                color = N_FOLD_AXIS_COLOR;
662                        } else {
663                                radius = polyhedron.getCirumscribedRadius();
664
665                                if (r.getFold() == 2) {
666                                        color = TWO_FOLD_AXIS_COLOR;
667                                } else if (r.getFold() == 3) {
668                                        color = THREE_FOLD_AXIS_COLOR;
669                                } else {
670                                        color = N_FOLD_AXIS_COLOR;
671                                }
672                        }
673
674
675                        Point3d center = rotationAxisAligner.getGeometricCenter();
676                        AxisAngle4d axisAngle = r.getAxisAngle();
677                        Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
678                        Vector3d refAxis = rotationAxisAligner.getRotationReferenceAxis();
679
680                        s.append(getSymmetryAxis(i, i+axes.size(), rotationGroup.getPointGroup(), r.getFold(), refAxis, radius, diameter, color, center, axis));
681                        i++;
682                }
683
684                return s.toString();
685        }
686
687        private Vector3d getAligmentVector(Point3d point, Vector3d axis) {
688                // for system with a single Cn axis
689                if (rotationGroup.getPointGroup().startsWith("C") || rotationGroup.getPointGroup().startsWith("D")) {
690                        // if axis is orthogonal to principal axis, use principal axis as reference axis
691                        if (axis.dot(rotationAxisAligner.getPrincipalRotationAxis()) < 0.1) {
692                                return rotationAxisAligner.getPrincipalRotationAxis();
693                        } else {
694                                return rotationAxisAligner.getRotationReferenceAxis();
695                        }
696                }
697
698                // for T, O, and I point groups find reference axis with respect to
699                // nearest polyhedron vertex
700                Vector3d ref = new Vector3d();
701                double dSqThreshold = 25;
702                double dSqMin = Double.MAX_VALUE;
703                Point3d refPoint = null;
704                // find nearest point on polyhedron as reference point,
705                // but do not choose a point on the same axis (therefore, we
706                // apply a distance threshold squared 5A*5A = 25A^2
707                for (Point3d v: getPolyhedronVertices()) {
708                        double dSq = point.distanceSquared(v);
709                        if (dSq > dSqThreshold) {
710                                if (dSq < dSqMin) {
711                                        dSqMin = dSq;
712                                        refPoint = v;
713                                }
714                        }
715                }
716
717
718                ref.sub(point, refPoint);
719
720                // this ref vector is usually not orthogonal to the
721                // axis. The following double-cross product makes it
722                // orthogonal.
723                ref.cross(axis, ref);
724                ref.cross(axis, ref); // note, done twice on purpose
725                ref.normalize();
726                return ref;
727        }
728
729        private String getSymmetryAxis(int i, int j, String pointGroup, int n, Vector3d referenceAxis, double radius, float diameter, String color, Point3d center, Vector3d axis) {
730                boolean drawPolygon = true;
731
732                Point3d p1 = new Point3d(axis);
733                p1.scaleAdd(-AXIS_SCALE_FACTOR * radius, center);
734
735                Point3d p2 = new Point3d(axis);
736                p2.scaleAdd(AXIS_SCALE_FACTOR * radius, center);
737
738                StringBuilder s = new StringBuilder();
739                s.append("draw");
740                s.append(" axesSymmetry");
741                s.append(name);
742                s.append(i);
743                s.append(" cylinder");
744                s.append(getJmolPoint(p1));
745                s.append(getJmolPoint(p2));
746                s.append("diameter ");
747                s.append(diameter);
748                s.append(" color ");
749                s.append(color);
750                s.append(" off;");
751
752                // calc. point to center symmetry symbols. They are offset by 0.01
753                // to avoid overlap with the polyhedron
754                p1 = new Point3d(axis);
755                p1.scaleAdd(-1.01*radius, center);
756
757                p2 = new Point3d(axis);
758                p2.scaleAdd(1.01*radius, center);
759
760                if (drawPolygon) {
761                        double polygonRadius = getMeanExtension() * 0.06;
762                        if (n == 2) {
763                                referenceAxis = getAligmentVector(p1, axis);
764                                s.append(getC2PolygonJmol(i, p1, referenceAxis, axis, color, polygonRadius, name));
765                                referenceAxis = getAligmentVector(p2, axis);
766                                s.append(getC2PolygonJmol(j, p2,  referenceAxis, axis, color, polygonRadius, name));
767                        } else if (n > 2) {
768                                referenceAxis = getAligmentVector(p1, axis);
769                                s.append(getPolygonJmol(i, p1, referenceAxis, axis, n, color, polygonRadius, name));
770                                referenceAxis = getAligmentVector(p2, axis);
771                                s.append(getPolygonJmol(j, p2, referenceAxis, axis, n, color, polygonRadius, name));
772                        }
773                }
774
775                return s.toString();
776        }
777
778        private static String getPolygonJmol(int index, Point3d center, Vector3d referenceAxis, Vector3d axis, int n, String color, double radius, String name) {
779                StringBuilder s = new StringBuilder();
780                s.append("draw axesSymbol");
781                s.append(name);
782                s.append(index);
783                s.append(" ");
784                s.append("polygon");
785                s.append(" ");
786                s.append(n+1);
787                s.append(getJmolPoint(center));
788
789                Vector3d[] vertexes = getPolygonVertices(axis, referenceAxis, center, n, radius);
790                // create vertex list
791                for (Vector3d v: vertexes) {
792                        s.append(getJmolPoint(v));
793                }
794
795                // create face list
796                s.append(n);
797                for (int i = 1; i <= n; i++) {
798                        s.append("[");
799                        s.append(0);
800                        s.append(" ");
801                        s.append(i);
802                        s.append(" ");
803                        if (i < n) {
804                                s.append(i+1);
805                        } else {
806                                s.append(1);
807                        }
808                        s.append(" ");
809                        s.append(7);
810                        s.append("]");
811                }
812
813                if (n == 2) {
814                        s.append("mesh off");
815                }
816                s.append(" color ");
817                s.append(color);
818                s.append(" off;");
819
820                return s.toString();
821        }
822
823        private static Vector3d[] getPolygonVertices(Vector3d axis, Vector3d referenceAxis, Point3d center, int n, double radius) {
824                Vector3d ref = new Vector3d(referenceAxis);
825                ref.scale(radius);
826
827                AxisAngle4d axisAngle = new AxisAngle4d(axis, 0);
828                Vector3d[] vectors = new Vector3d[n];
829                Matrix4d m = new Matrix4d();
830
831                for (int i = 0; i < n; i++) {
832                        axisAngle.angle = i * 2 * Math.PI/n;
833                        vectors[i] = new Vector3d(ref);
834                        m.set(axisAngle);
835                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
836                        m.setElement(3, 3, 1.0);
837                        m.transform(vectors[i]);
838                        vectors[i].add(center);
839                }
840                return vectors;
841        }
842
843        private static String getC2PolygonJmol(int index, Point3d center, Vector3d referenceAxis, Vector3d axis, String color, double radius, String name) {
844                StringBuilder s = new StringBuilder();
845                int n = 10;
846                s.append("draw axesSymbol");
847                s.append(name);
848                s.append(index);
849                s.append(" ");
850                s.append("polygon");
851                s.append(" ");
852                s.append(n-1);
853                s.append(getJmolPoint(center));
854
855                Vector3d[] vertexes = getC2PolygonVertices(axis, referenceAxis, center, n, radius);
856                // create vertex list
857                for (Vector3d v: vertexes) {
858                        s.append(getJmolPoint(v));
859                }
860
861                // create face list
862                s.append(n-2);
863
864                for (int i = 1; i < n-1; i++) {
865                        s.append("[");
866                        s.append(0);
867                        s.append(" ");
868                        s.append(i);
869                        s.append(" ");
870                        if (i < n-2) {
871                                s.append(i+1);
872                        } else {
873                                s.append(1);
874                        }
875                        s.append(" ");
876                        s.append(7);
877                        s.append("]");
878                }
879
880                s.append("color ");
881                s.append(color);
882                s.append(" off;");
883
884                return s.toString();
885        }
886        private static Vector3d[] getC2PolygonVertices(Vector3d axis, Vector3d referenceAxis, Point3d center, int n, double radius) {
887                Vector3d ref = new Vector3d(referenceAxis);
888                ref.scale(4*radius);
889
890                AxisAngle4d axisAngle = new AxisAngle4d(axis, 0);
891                int k = n / 2;
892                // draw arc 1/6 of a full circle
893                int f = 6;
894                Vector3d[] vectors = new Vector3d[n-2];
895                Matrix4d m = new Matrix4d();
896
897                // first point of arc
898                axisAngle.angle = (k+0.5) * 2 * Math.PI/(f*k);
899                Vector3d begin = new Vector3d(ref);
900                m.set(axisAngle);
901                // make sure matrix element m33 is 1.0. It's 0 on Linux.
902                m.setElement(3, 3, 1.0);
903                m.transform(begin);
904
905                // last point of arc
906                axisAngle.angle = (2*k-1+0.5) * 2 * Math.PI/(f*k);
907                Vector3d end = new Vector3d(ref);
908                m.set(axisAngle);
909                // make sure matrix element m33 is 1.0. It's 0 on Linux.
910                m.setElement(3, 3, 1.0);
911                m.transform(end);
912
913                // center of arc
914                Vector3d arcCenter = new Vector3d();
915                arcCenter.interpolate(begin, end, 0.5);
916                arcCenter.negate();
917
918                // add translation component
919                Vector3d trans =  new Vector3d();
920                trans.sub(center, arcCenter);
921
922                // draw arc (1/6 of a full circle)
923                for (int i = 0; i < k; i++) {
924                        axisAngle.angle = (k + i + 0.5) * 2 * Math.PI/(f*k);
925                        vectors[i] = new Vector3d(ref);
926                        m.set(axisAngle);
927                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
928                        m.setElement(3, 3, 1.0);
929                        m.transform(vectors[i]);
930                        vectors[i].add(arcCenter);
931                        vectors[i].add(center);
932                }
933                // in reverse order, draw reflected half of arc (1/6 of full circle)
934                // don't draw first and last element, since the are already part of the previous arc
935                for (int i = k; i < 2*k-2; i++) {
936                        axisAngle.angle = (f/2*k + i + 1.5) * 2 * Math.PI/(f*k);
937                        vectors[i] = new Vector3d(ref);
938                        m.set(axisAngle);
939                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
940                        m.setElement(3, 3, 1.0);
941                        m.transform(vectors[i]);
942                        vectors[i].sub(arcCenter);
943                        vectors[i].add(center);
944                }
945                return vectors;
946        }
947
948
949        private List<Rotation> getUniqueAxes() {
950                List<Rotation> uniqueRotations = new ArrayList<Rotation>();
951
952                for (int i = 0, n = rotationGroup.getOrder(); i < n; i++) {
953                        Rotation rotationI = rotationGroup.getRotation(i);
954                        AxisAngle4d axisAngleI = rotationI.getAxisAngle();
955                        Vector3d axisI = new Vector3d(axisAngleI.x, axisAngleI.y, axisAngleI.z);
956
957                        boolean redundant = false;
958                        for (Rotation r: uniqueRotations) {
959                                AxisAngle4d axisAngleJ = r.getAxisAngle();
960                                Vector3d axisJ = new Vector3d(axisAngleJ.x, axisAngleJ.y, axisAngleJ.z);
961                                if (Math.abs(axisI.dot(axisJ)) > 0.99) {
962                                        redundant = true;
963                                        break;
964                                }
965                        }
966
967                        if (! redundant) {
968                                uniqueRotations.add(rotationI);
969                        }
970                }
971                return uniqueRotations;
972        }
973
974        private String drawHeader(String text, String color) {
975                StringBuilder s = new StringBuilder();
976                s.append("set echo top center;");
977                s.append("color echo ");
978                s.append(color);
979                s.append(";");
980                s.append("font echo 24 sanserif;");
981                s.append("echo ");
982                s.append(text);
983                s.append(";");
984                return s.toString();
985        }
986
987        private String deleteHeader() {
988                return "set echo top center;echo ;";
989        }
990
991        private String drawFooter(String text, String color) {
992                StringBuilder s = new StringBuilder();
993                s.append("set echo bottom center;");
994                s.append("color echo ");
995                s.append(color);
996                s.append(";");
997                s.append("font echo 24 sanserif;");
998                s.append("echo ").append(text);
999                //s.append("echo Point group ");
1000                //s.append(rotationGroup.getPointGroup());
1001                s.append(";");
1002                return s.toString();
1003        }
1004
1005        private String setCentroid() {
1006                // calculate center of rotation
1007        //      Point3d centroid = axisTransformation.getGeometricCenter();
1008                Point3d centroid = rotationAxisAligner.getCentroid();
1009
1010                // set centroid
1011                StringBuilder s = new StringBuilder();
1012                s.append("center");
1013                s.append(getJmolPoint(centroid));
1014                s.append(";");
1015                return s.toString();
1016        }
1017}