Back to index

wims  3.65+svn20090927
Molecule.java
Go to the documentation of this file.
00001 /*
00002     WIMSchem Elements: Chemistry molecular diagram drawing tool.
00003     
00004     (c) 2005 Dr. Alex M. Clark
00005     
00006     Released as GNUware, under the Gnu Public License (GPL)
00007     
00008     See www.gnu.org for details.
00009 */
00010 
00011 package WIMSchem;
00012 
00013 import java.util.*;
00014 
00015 /*
00016     Internal molecule datastructure, which captures the basic information about a molecule as envisaged by the traditional diagram
00017     notation. The core data is stored as a list of labelled nodes (atoms) and labelled edges (bonds). Besides very basic editing and
00018     convenient object manipulation, a few computed properties are available, some of which are cached.
00019 */
00020 
00021 public class Molecule
00022 {
00023     class Atom
00024     {
00025        String Element;
00026        double X,Y;
00027        int Charge,Unpaired;
00028        int HExplicit; 
00029        int MapNum;
00030     }
00031     class Bond
00032     {
00033        int From,To;
00034        int Order,Type;
00035     }
00036     
00037     public static final int HEXPLICIT_UNKNOWN=-1;
00038     
00039     public static final int BONDTYPE_NORMAL=0;
00040     public static final int BONDTYPE_INCLINED=1;
00041     public static final int BONDTYPE_DECLINED=2;
00042     public static final int BONDTYPE_UNKNOWN=3;
00043     
00044     public static final String[] ELEMENTS={null,
00045        "H","He","Li","Be","B","C","N","O","F","Ne","Na","Mg","Al","Si","P","S","Cl","Ar","K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co",
00046        "Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb","Sr","Y","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I",
00047        "Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W","Re","Os","Ir","Pt","Au",
00048        "Hg","Tl","Pb","Bi","Po","At","Rn","Fr","Ra","Ac","Th","Pa","U","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr"};
00049 
00050     public static final double[] WEIGHTS={0.0,
00051        1.00794,4.002602,6.941,9.01218,10.811,12.011,14.00674,15.9994,18.998403,20.1797,
00052        22.989768,24.305,26.981539,28.0855,30.973762,32.066,35.4527,39.948,39.0983,40.078,
00053        44.95591,47.88,50.9415,51.9961,54.93805,55.847,58.9332,58.6934,63.546,65.39,
00054        69.723,72.61,74.92159,78.96,79.904,83.8,85.4678,87.62,88.90585,91.224,
00055        92.90638,95.94,97.9072,101.07,102.9055,106.42,107.8682,112.411,114.818,118.71,
00056        121.760,127.6,126.90447,131.29,132.90543,137.327,138.9055,140.115,140.90765,144.24,
00057        144.9127,150.36,151.965,157.25,158.92534,162.50,164.93032,167.26,168.93421,173.04,
00058        174.967,178.49,180.9479,183.84,186.207,190.23,192.22,195.08,196.96654,200.59,
00059        204.3833,207.2,208.98037,208.9824,209.9871,222.0176,223.0197,226.0254,227.0278,232.0381,
00060        231.03588,238.0289,237.048,244.0642,243.0614,247.0703,247.0703,251.0796,252.083,257.0951,
00061        258.1,259.1009,262.11};                                                                                                              
00062                                 
00063        
00064 
00065     boolean invalMinMax=false;
00066     double minX=0,minY=0,maxX=0,maxY=0;
00067     ArrayList<Atom> atoms=new ArrayList<Atom>();
00068     ArrayList<Bond> bonds=new ArrayList<Bond>();
00069     
00070     // computed graph topology properties
00071     int graph[][]=null;
00072     int ringID[]=null,compID[]=null,priority[]=null;
00073     ArrayList<Object> ring5=null,ring6=null,ring7=null; // common ring sizes
00074     
00075     // computed stereochemistry properties
00076     public static final int STEREO_NONE=0; // topology does not allow any stereoisomers
00077     public static final int STEREO_POS=1; // R or Z, depending on type
00078     public static final int STEREO_NEG=2; // S or E, depending on type
00079     public static final int STEREO_UNKNOWN=3; // topology allows stereoisomers, but wedges/geometry do not resolve which
00080     int chiral[]=null; // (per atom)
00081     int cistrans[]=null; // (per bond)
00082 
00083     // data for calculation of implicit hydrogens
00084     static final String[] HYVALENCE_EL={"C","N","O","S","P"};
00085     static final int[]   HYVALENCE_VAL={ 4,  3,  2,  2,  3 };
00086     
00087     // ------------------ public functions --------------------
00088 
00089     public Molecule()
00090     {
00091     }
00092 
00093     // returns an equivalent replica of the molecule
00094     public Molecule Clone()
00095     {
00096        Molecule mol=new Molecule();
00097        for( int n=1;n<=NumAtoms();n++)
00098        {
00099            Atom a=atoms.get(n-1);
00100            int num=mol.AddAtom(a.Element,a.X,a.Y,a.Charge,a.Unpaired);
00101            mol.SetAtomHExplicit(num,AtomHExplicit(n));
00102            mol.SetAtomMapNum(num,AtomMapNum(n));
00103        }
00104        for( int n=1;n<=NumBonds();n++)
00105        {
00106            Bond b=bonds.get(n-1);
00107            mol.AddBond(b.From,b.To,b.Order,b.Type);
00108        }
00109        return mol;
00110     }
00111 
00112     
00113     // ----- fetching information directly from the datastructures ----
00114         
00115     public int NumAtoms() {return atoms.size();}
00116     public int NumBonds() {return bonds.size();}
00117     
00118     public String AtomElement(int N) {return atoms.get(N-1).Element;}
00119     public double AtomX(int N) {return atoms.get(N-1).X;}
00120     public double AtomY(int N) {return atoms.get(N-1).Y;}
00121     public int AtomCharge(int N) {return atoms.get(N-1).Charge;}
00122     public int AtomUnpaired(int N) {return atoms.get(N-1).Unpaired;}
00123     public int AtomHExplicit(int N) {return atoms.get(N-1).HExplicit;}
00124     public int AtomMapNum(int N) {return atoms.get(N-1).MapNum;}
00125     
00126     public int BondFrom(int N) {return bonds.get(N-1).From;}
00127     public int BondTo(int N) {return bonds.get(N-1).To;}
00128     public int BondOrder(int N) {return bonds.get(N-1).Order;}
00129     public int BondType(int N) {return bonds.get(N-1).Type;}
00130     
00131     // ----- modifying information in the datastructures ----
00132 
00133     // adds an atom to the end of the molecule; it is fairly reasonable to assume that the new one will be at position==NumAtoms(),
00134     // and that all other atoms will be in the same places, but the position is returned for convenience
00135     public int AddAtom(String Element,double X,double Y) {return AddAtom(Element,X,Y,0,0);}
00136     public int AddAtom(String Element,double X,double Y,int Charge,int Unpaired)
00137     {
00138        if( X<minX || NumAtoms()==0) minX=X;
00139        if( X>maxX || NumAtoms()==0) maxX=X;
00140        if( Y<minY || NumAtoms()==0) minY=Y;
00141        if( Y>maxY || NumAtoms()==0) maxY=Y;
00142        
00143        Atom a=new Atom();
00144        a.Element=Element;
00145        a.X=X;
00146        a.Y=Y;
00147        a.Charge=Charge;
00148        a.Unpaired=Unpaired;
00149        a.HExplicit=HEXPLICIT_UNKNOWN;
00150        atoms.add(a);
00151        TrashGraph();
00152        return atoms.size();
00153     }
00154     public void SetAtomElement(int N,String Element) {atoms.get(N-1).Element=Element; TrashPriority();}
00155     public void SetAtomPos(int N,double X,double Y) 
00156     {
00157        atoms.get(N-1).X=X; 
00158        atoms.get(N-1).Y=Y;
00159        invalMinMax=true;
00160        TrashStereo();
00161     }
00162     public void SetAtomCharge(int N,int Charge) {atoms.get(N-1).Charge=Charge; TrashPriority();}
00163     public void SetAtomUnpaired(int N,int Unpaired) {atoms.get(N-1).Unpaired=Unpaired; TrashPriority();}
00164     public void SetAtomHExplicit(int N,int HExplicit) {atoms.get(N-1).HExplicit=HExplicit; TrashPriority();}
00165     public void SetAtomMapNum(int N,int MapNum) {atoms.get(N-1).MapNum=MapNum;}
00166     
00167     // adds an atom the end of the molecule; the position of the new bond is returned for convenience
00168     public int AddBond(int From,int To,int Order) {return AddBond(From,To,Order,BONDTYPE_NORMAL);}
00169     public int AddBond(int From,int To,int Order,int Type)
00170     {
00171        Bond b=new Bond();
00172        b.From=From;
00173        b.To=To;
00174        b.Order=Order;
00175        b.Type=Type;
00176        bonds.add(b);
00177        TrashGraph();
00178        return bonds.size();
00179     }
00180     
00181     public void SetBondFromTo(int N,int From,int To) 
00182     {
00183        bonds.get(N-1).From=From; 
00184        bonds.get(N-1).To=To;
00185        TrashGraph();
00186     }
00187     public void SetBondOrder(int N,int Order) {bonds.get(N-1).Order=Order; TrashPriority();}
00188     public void SetBondType(int N,int Type) {bonds.get(N-1).Type=Type; TrashStereo();}
00189     
00190     // deletes the indicated atom, and also any bonds which were connected to it
00191     public void DeleteAtomAndBonds(int N)
00192     {
00193        for( int n=0;n<NumBonds();)
00194        {
00195            Bond b=bonds.get(n);
00196            if( b.From==N || b.To==N) bonds.remove(n);
00197            else
00198            {
00199               if( b.From>N) b.From--;
00200               if( b.To>N) b.To--;
00201               n++;
00202            }
00203        }
00204        atoms.remove(N-1);
00205        invalMinMax=true;
00206        TrashGraph();
00207     }
00208     public void DeleteBond(int N)
00209     {
00210        bonds.remove(N-1);
00211        TrashGraph();
00212     }
00213      
00214     // ----- fetching information which is computed from the underlying data, usually cached -----
00215      
00216     public double MinX() {if( invalMinMax) DetermineMinMax(); return minX;}
00217     public double MaxX() {if( invalMinMax) DetermineMinMax(); return maxX;}
00218     public double MinY() {if( invalMinMax) DetermineMinMax(); return minY;}
00219     public double MaxY() {if( invalMinMax) DetermineMinMax(); return maxY;}
00220     public double RangeX() {return maxX-minX;}
00221     public double RangeY() {return maxY-minY;}
00222     
00223     // return the index of the bond, if any, which connects the two indicated atoms
00224     public int FindBond(int A1,int A2)
00225     {
00226        for( int n=1;n<=NumBonds();n++) 
00227        {
00228            if( BondFrom(n)==A1 && BondTo(n)==A2) return n;
00229            if( BondTo(n)==A1 && BondFrom(n)==A2) return n;
00230        }
00231        return 0;
00232     }
00233 
00234     // for a bond, returns the end which is not==Ref; return value will be From,To or 0    
00235     public int BondOther(int N,int Ref)
00236     {
00237        if( BondFrom(N)==Ref) return BondTo(N);
00238        if( BondTo(N)==Ref) return BondFrom(N);
00239        return 0;
00240     }
00241 
00242     // returns whether an atom logically ought to be drawn "explicitly";if false, it is a carbon which should be just a dot
00243     public boolean AtomExplicit(int N) 
00244     {
00245        if( atoms.get(N-1).Element.compareTo("C")!=0 || atoms.get(N-1).Charge!=0 || atoms.get(N-1).Unpaired!=0) return true;
00246        for( int n=0;n<bonds.size();n++) if( bonds.get(n).From==N || bonds.get(n).To==N) return false;
00247        return true;
00248     }
00249     
00250     // uses either explicit or computed number to determine how many hydrogens the atom has; the field for explicit hydrogens takes
00251     // absolute preference, if it has its default value of 'unknown', the number is computed by looking up the hydrogen capacity for
00252     // the element (most of which are zero), subtracting the total of bond orders, then returning the difference, or zero; the calculation
00253     // tends to err on the side of too few, since the concept is just an aid to drawing organic structures, not a cheminformatic attempt
00254     // to compensate for 2 1/2 decades of bad file formats
00255     // (note: returns "implicit"+"explicit", but does NOT count "actual" hydrogens, i.e. those which have their own atom nodes)
00256     public int AtomHydrogens(int N) 
00257     {
00258        int hy=AtomHExplicit(N);
00259        if( hy!=HEXPLICIT_UNKNOWN) return hy;
00260        for( int n=0;n<HYVALENCE_EL.length;n++) if( HYVALENCE_EL[n].compareTo(AtomElement(N))==0) {hy=HYVALENCE_VAL[n]; break;}
00261        if( hy==HEXPLICIT_UNKNOWN) return 0;
00262        int ch=AtomCharge(N);
00263        if( AtomElement(N).compareTo("C")==0) ch=-Math.abs(ch);
00264        hy+=ch-AtomUnpaired(N);
00265        for( int n=1;n<=NumBonds();n++) if( BondFrom(n)==N || BondTo(n)==N) hy-=BondOrder(n);
00266        return hy<0 ? 0 : hy;
00267     }
00268      
00269     // returns the numerical ID of the ring block in which the atom resides, or 0 if it is not in a ring   
00270     public int AtomRingBlock(int N)
00271     {
00272        if( graph==null) BuildGraph();
00273        if( ringID==null) BuildRingID();
00274        return ringID[N-1];
00275     }
00276     // returns whether or not a bond resides in a ring (i.e. ring block of each end is the same and nonzero)
00277     public boolean BondInRing(int N)
00278     {
00279        int r1=AtomRingBlock(BondFrom(N)),r2=AtomRingBlock(BondTo(N));
00280        return r1>0 && r1==r2;
00281     }
00282     
00283     // returns the connected component ID of the atom, always 1 or more
00284     public int AtomConnComp(int N)
00285     {
00286        if( graph==null) BuildGraph();
00287        if( compID==null) BuildConnComp();
00288        return compID[N-1];
00289     }
00290     
00291     // returns the number of neighbours of an atom
00292     public int AtomAdjCount(int N)
00293     {
00294        if( graph==null) BuildGraph();
00295        return graph[N-1].length;
00296     }
00297     
00298     // returns the actual adjacency list of an atom; unfortunately this has to make a clone of the array, since the numbers are converted
00299     // to 1-based indices, and we would not want callers modifying it anyway
00300     public int[] AtomAdjList(int N)
00301     {
00302        if( graph==null) BuildGraph();
00303        int[] adj=graph[N-1].clone();
00304        for( int n=0;n<adj.length;n++) adj[n]++;
00305        return adj;
00306     }
00307     
00308     // returns atom-based priority according to the Cahn-Ingold-Prelog rules
00309     public int AtomPriority(int N)
00310     {
00311        if( graph==null) BuildGraph();
00312        if( compID==null) BuildConnComp();
00313        if( priority==null) BuildPriority();
00314        return priority[N-1];
00315     }
00316     
00317     // return the tetrahedral chirality of the given atom
00318     public int AtomChirality(int N)
00319     {
00320        if( graph==null) BuildGraph();
00321        if( compID==null) BuildConnComp();
00322        if( priority==null) BuildPriority();
00323        if( chiral==null) BuildChirality();
00324        return chiral[N-1];
00325     }
00326     
00327     // return the cis/trans style stereochemistry of the given bond
00328     public int BondStereo(int N)
00329     {
00330        if( graph==null) BuildGraph();
00331        if( compID==null) BuildConnComp();
00332        if( priority==null) BuildPriority();
00333        if( cistrans==null) BuildCisTrans();
00334        return cistrans[N-1];
00335     }
00336     
00337     // returns _all_ rings of indicated size; each item in the array list is an array of int[Size], a consecutively ordered array of atom 
00338     // numbers; uses a recursive depth first search, which must be bounded above by Size being small in order to avoid exponential blowup
00339     public int[][] FindRingSize(int Size) 
00340     {
00341        ArrayList<Object> rings=null;
00342        if( Size==5 && ring5!=null) rings=ring5;
00343        if( Size==6 && ring6!=null) rings=ring6;
00344        if( Size==7 && ring7!=null) rings=ring7;
00345     
00346        if( rings==null)
00347        {
00348            if( graph==null) BuildGraph();
00349            if( ringID==null) BuildRingID();
00350 
00351            rings=new ArrayList<Object>();
00352            for( int n=1;n<=NumAtoms();n++) if( ringID[n-1]>0) 
00353            {
00354               int path[]=new int[Size];
00355               path[0]=n;
00356               RecursiveRingFind(path,1,Size,ringID[n-1],rings);
00357            }
00358 
00359            if( Size==5) ring5=rings;
00360            if( Size==6) ring6=rings;
00361            if( Size==7) ring7=rings;
00362        }
00363        
00364        int ret[][]=new int[rings.size()][];
00365        for( int n=0;n<rings.size();n++) ret[n]=((int[])rings.get(n)).clone();
00366        return ret;
00367     }
00368     
00369     // compare to another molecule instance; null is equivalent to empty; return values:
00370     //     -1 if Mol < Other
00371     //     0 if Mol == Other
00372     //     1 if Mol > Other 
00373     //  can be used to sort molecules, albeit crudely; does not do any kind of canonical preprocessing
00374     public int CompareTo(Molecule Other)
00375     {
00376        if( Other==null) return NumAtoms()==0 ? 0 : 1; // null is equivalent to empty
00377        if( NumAtoms()<Other.NumAtoms()) return -1;
00378        if( NumAtoms()>Other.NumAtoms()) return 1;
00379        if( NumBonds()<Other.NumBonds()) return -1;
00380        if( NumBonds()>Other.NumBonds()) return 1;
00381        
00382        for( int n=1;n<=NumAtoms();n++)
00383        {
00384            int c=AtomElement(n).compareTo(Other.AtomElement(n));
00385            if( c!=0) return c;
00386            if( AtomX(n)<Other.AtomX(n)) return -1; if( AtomX(n)>Other.AtomX(n)) return 1;
00387            if( AtomY(n)<Other.AtomY(n)) return -1; if( AtomY(n)>Other.AtomY(n)) return 1;
00388            if( AtomCharge(n)<Other.AtomCharge(n)) return -1; if( AtomCharge(n)>Other.AtomCharge(n)) return 1;
00389            if( AtomUnpaired(n)<Other.AtomUnpaired(n)) return -1; if( AtomUnpaired(n)>Other.AtomUnpaired(n)) return 1;
00390            if( AtomHExplicit(n)<Other.AtomHExplicit(n)) return -1; if( AtomHExplicit(n)>Other.AtomHExplicit(n)) return 1;
00391        }
00392        for( int n=1;n<=NumBonds();n++)
00393        {
00394            if( BondFrom(n)<Other.BondFrom(n)) return -1; if( BondFrom(n)>Other.BondFrom(n)) return 1;
00395            if( BondTo(n)<Other.BondTo(n)) return -1; if( BondTo(n)>Other.BondTo(n)) return 1;
00396            if( BondOrder(n)<Other.BondOrder(n)) return -1; if( BondOrder(n)>Other.BondOrder(n)) return 1;
00397            if( BondType(n)<Other.BondType(n)) return -1; if( BondType(n)>Other.BondType(n)) return 1;
00398        }
00399        
00400        return 0;
00401     }
00402 /*
00403     public static String GetScore(Molecule Q,Molecule A){// Q = correct A= student answer
00404        if( A == null ){return "empty_answer";}
00405        if( Q == null ){return "empty_question";}
00406        int max_a=A.NumAtoms();
00407        int max_q=Q.NumAtoms();
00408        if( max_a<max_q){return "too_few_atoms" + (max_q - max_a);}
00409        if( max_a>max_q){return "too_many_atoms" + (max_q - max_a) ;}
00410        if( A.NumBonds()<Q.NumBonds()){return "too_few_bonds";}
00411        if( A.NumBonds()>Q.NumBonds()){return "too_many_bonds";}
00412        int hit=0;
00413        int hydro_a;int hydro_q;
00414        int hydrosum_a=0;int hydrosum_q=0;
00415        String element_a="";String element_q="";
00416        String next_a="";String next_q="";
00417        String prev_a="";String prev_q="";
00418        int type_a;int type_q;
00419        int order_a;int order_q;
00420        int charge_a;int charge_q;
00421        int hexplicit_a;int hexplicit_q;
00422        int radical_a;int radical_q;
00423 
00424        for( int n=1 ; n<=max_q ; n++ ){
00425        // loop the Question
00426            charge_q=Q.AtomCharge(n);
00427            radical_q=Q.AtomUnpaired(n);
00428            hexplicit_q=Q.AtomHExplicit(n);
00429            element_q=Q.AtomElement(n);
00430            hydro_q=Q.AtomHydrogens(n);
00431            if(n<max_q){ // these are Bonds
00432               type_q=Q.BondType(n);
00433               order_q=Q.BondOrder(n);
00434               next_q=Q.AtomElement(Q.BondTo(n));
00435               prev_q=Q.AtomElement(Q.BondFrom(n));
00436            }
00437            else
00438            {
00439               type_q=100;
00440               order_q=100;
00441               next_q="not";
00442               prev_q="not";
00443            }
00444            hydrosum_q=hydrosum_q+hydro_q;
00445            hit=0;
00446            for( int m=1 ; m<=max_a ; m++ ){
00447               // loop the answer
00448               charge_a=A.AtomCharge(m);
00449               radical_a=A.AtomUnpaired(m);
00450               hexplicit_a=A.AtomHExplicit(m);
00451               element_a=A.AtomElement(m);
00452               hydro_a=A.AtomHydrogens(m);
00453               if(m<max_a){// these are Bonds
00454                   type_a=A.BondType(m);
00455                   order_a=A.BondOrder(m);
00456                   next_a=A.AtomElement(A.BondTo(m));
00457                   prev_a=A.AtomElement(A.BondFrom(m));
00458               }
00459               else
00460               {
00461                   type_a=100;
00462                   order_a=100;
00463                   next_a="not";
00464                   next_q="not";
00465               }
00466               if(n==1){hydrosum_a=hydrosum_a+hydro_a;}
00467               System.out.println("Question="+hydro_q+"=?="+hydro_a);
00468               if( element_a.equals(element_q) && 
00469               //    next_a.equals(next_q) && 
00470               //    prev_a.equals(prev_q) &&
00471                   hydro_a == hydro_q && 
00472                   order_a == order_q && 
00473                   type_a == type_q && 
00474                   radical_a == radical_q &&
00475                   charge_a == charge_q &&
00476                   hexplicit_a == hexplicit_q){   
00477                   hit=1;
00478                   //if(n!=1){break;}
00479               }
00480            }
00481            if(hit==0){return "found_error_for_element "+element_q;}
00482        }
00483        if(hydrosum_a != hydrosum_q){return "Hydrogen_mismatch h_a="+hydrosum_a+" h_q="+hydrosum_q;}
00484        return "OK";
00485     }
00486 */    
00487 
00488     // lookup the atomic number for the element, or return 0 if not in the periodic table
00489     int AtomicNumber(int N) {return AtomicNumber(AtomElement(N));}
00490     int AtomicNumber(String El)
00491     {
00492        for( int n=1;n<ELEMENTS.length;n++) if( ELEMENTS[n].compareTo(El)==0) return n;
00493        return 0;
00494     }
00495 
00496     // takes the indicated molecular fragment and appends it to the end of the current molecule, with order preserved
00497     public void Append(Molecule Frag)
00498     {
00499        int base=NumAtoms();
00500        for( int n=1;n<=Frag.NumAtoms();n++)
00501        {
00502            int num=AddAtom(Frag.AtomElement(n),Frag.AtomX(n),Frag.AtomY(n),Frag.AtomCharge(n),Frag.AtomUnpaired(n));
00503            SetAtomHExplicit(num,Frag.AtomHExplicit(n));
00504            SetAtomMapNum(num,Frag.AtomMapNum(n));
00505        }
00506        for( int n=1;n<=Frag.NumBonds();n++)
00507        {
00508            AddBond(Frag.BondFrom(n)+base,Frag.BondTo(n)+base,Frag.BondOrder(n),Frag.BondType(n));
00509        }
00510     }
00511 
00512     // returns a molecule which is smaller than the current one, according to the mask (which must be of size #atoms); boundary cases
00513     // are a null molecule or cloned copy
00514     public Molecule Subgraph(boolean[] Mask)
00515     {
00516        int[] invidx=new int[NumAtoms()];
00517        int sum=0;
00518        for( int n=0;n<NumAtoms();n++) if( Mask[n]) invidx[n]=++sum; else invidx[n]=0;
00519        if( sum==0) return new Molecule();
00520        if( sum==NumAtoms()) return Clone();
00521        
00522        Molecule frag=new Molecule();
00523        for( int n=1;n<=NumAtoms();n++) if( Mask[n-1])
00524        {
00525            int num=frag.AddAtom(AtomElement(n),AtomX(n),AtomY(n),AtomCharge(n),AtomUnpaired(n));
00526            frag.SetAtomHExplicit(num,AtomHExplicit(n));
00527            frag.SetAtomMapNum(num,AtomMapNum(n));
00528        }
00529        for( int n=1;n<=NumBonds();n++) 
00530        {
00531            int from=invidx[BondFrom(n)-1],to=invidx[BondTo(n)-1];
00532            if( from>0 && to>0) frag.AddBond(from,to,BondOrder(n),BondType(n));
00533        }
00534        
00535        return frag;
00536     }
00537 
00538     
00539     // convenient debugging mechanism
00540     public void Dump()
00541     {
00542        System.out.println("#Atoms="+NumAtoms()+" #Bonds="+NumBonds());
00543        for( int n=0;n<NumAtoms();n++)
00544        {
00545            Atom a=atoms.get(n);
00546            System.out.println(" A"+(n+1)+": "+a.Element+"="+a.X+","+a.Y+";"+a.Charge+","+a.Unpaired);
00547        }
00548        for( int n=0;n<NumBonds();n++)
00549        {
00550            Bond b=bonds.get(n);
00551            System.out.println(" B"+(n+1)+": "+b.From+"-"+b.To+"="+b.Order+","+b.Type);
00552        }
00553     }
00554 
00555     // ------------------ private functions --------------------
00556     
00557     // for when the connectivity changes somehow; abandon graph properties
00558     void TrashGraph()
00559     {
00560        graph=null;
00561        ringID=null;
00562        ring5=null;
00563        ring6=null;
00564        ring7=null;
00565        compID=null;
00566        TrashPriority();
00567     }
00568     void TrashPriority() 
00569     {
00570        priority=null;
00571        TrashStereo();
00572     }
00573     void TrashStereo()
00574     {
00575        chiral=null;
00576        cistrans=null;
00577     }
00578     
00579     // update the cache of atom range
00580     void DetermineMinMax()
00581     {
00582        invalMinMax=false;
00583        if( NumAtoms()==0) 
00584        {
00585            minX=maxX=minY=maxY=0; 
00586            return;
00587        }
00588        minX=maxX=AtomX(1);
00589        minY=maxY=AtomY(1);
00590        for( int n=2;n<=NumAtoms();n++)
00591        {
00592            double x=AtomX(n),y=AtomY(n);
00593            minX=Math.min(minX,x);
00594            maxX=Math.max(maxX,x);
00595            minY=Math.min(minY,y);
00596            maxY=Math.max(maxY,y);
00597        }
00598     }
00599 
00600     // update the cache of the molecular graph, in adjacency format, rather than the edge list format; the former is much better for
00601     // graph algorithms, the latter much more suited to capturing sketcher-ish information content; note that the graph indices are
00602     // zero-based
00603     void BuildGraph()
00604     {
00605        graph=new int[NumAtoms()][];
00606        for( int n=1;n<=NumBonds();n++)
00607        {
00608            int bf=BondFrom(n)-1,bt=BondTo(n)-1;
00609            if( bf==bt) continue; // ouch!
00610            int lf=graph[bf]==null ? 0 : graph[bf].length,lt=graph[bt]==null ? 0 : graph[bt].length;
00611            int bl[]=new int[lf+1];
00612            for( int i=0;i<lf;i++) bl[i]=graph[bf][i];
00613            bl[lf]=bt;
00614            graph[bf]=bl;
00615            bl=new int[lt+1];
00616            for( int i=0;i<lt;i++) bl[i]=graph[bt][i];
00617            bl[lt]=bf;
00618            graph[bt]=bl;
00619        }
00620        for( int n=0;n<NumAtoms();n++) if( graph[n]==null) graph[n]=new int[0];
00621        
00622        /*for( int n=0;n<NumAtoms();n++) 
00623        {
00624            System.out.print("#"+n+":");
00625            for( int i=0;i<graph[n].length;i++) System.out.print(" "+graph[n][i]);
00626            System.out.println();
00627        }*/
00628     }
00629     
00630     // passes over the graph establishing which component each atom belongs to
00631     void BuildConnComp()
00632     {
00633        compID=new int[NumAtoms()];
00634        for( int n=0;n<NumAtoms();n++) compID[n]=0;
00635        int comp=1;
00636        compID[0]=comp;
00637 
00638        // (not very efficient, should use a stack-based depth first search)
00639        while (true)
00640        {
00641            boolean anything=false;
00642            for( int n=0;n<NumAtoms();n++) if( compID[n]>0)
00643            {
00644               for( int i=0;i<graph[n].length;i++) if( compID[graph[n][i]]==0)
00645               {
00646                   compID[graph[n][i]]=comp;
00647                   anything=true;
00648               }
00649            }
00650            
00651            if( !anything)
00652            {
00653               for( int n=0;n<NumAtoms();n++) if( compID[n]==0) {compID[n]=++comp; anything=true; break;}
00654               if( !anything) break;
00655            }
00656        }
00657     }
00658 
00659     // generates Cahn-Ingold-Prelog priority for each atom, where degeneracies are indicated by having the same number
00660     void BuildPriority()
00661     {
00662        // build a graph representation which has entries replicated according to bond order
00663        int cipgr[][]=new int[NumAtoms()][];
00664        for( int n=0;n<NumAtoms();n++) if( cipgr[n]==null) 
00665        {
00666            cipgr[n]=new int[AtomHydrogens(n+1)];
00667            for( int i=0;i<cipgr[n].length;i++) cipgr[n][i]=-1;
00668        }
00669        for( int n=1;n<=NumBonds();n++)
00670        {
00671            int bf=BondFrom(n)-1,bt=BondTo(n)-1,bo=BondOrder(n);
00672            if( bf==bt || bo==0) continue;
00673            int lf=cipgr[bf].length,lt=cipgr[bt].length;
00674            int bl[]=new int[lf+bo];
00675            for( int i=0;i<lf;i++) bl[i]=cipgr[bf][i];
00676            for( int i=0;i<bo;i++) bl[lf++]=bt;
00677            cipgr[bf]=bl;
00678            bl=new int[lt+bo];
00679            for( int i=0;i<lt;i++) bl[i]=cipgr[bt][i];
00680            for( int i=0;i<bo;i++) bl[lt++]=bf;
00681            cipgr[bt]=bl;
00682        }
00683 
00684        // seed the priorities with atomic number
00685        priority=new int[NumAtoms()];
00686        boolean anyActualH=false;
00687        for( int n=0;n<NumAtoms();n++) {priority[n]=AtomicNumber(n+1); if( priority[n]==1) anyActualH=true;}
00688        
00689        // pass through and reassign priorities as many times as necessary, until no change
00690        int prigr[][]=new int[NumAtoms()][];
00691        while (true)
00692        {
00693            // make an equivalent to cipgr which has priorities instead of indices
00694            for( int n=0;n<NumAtoms();n++) 
00695            {
00696               prigr[n]=new int[cipgr[n].length];
00697               for( int i=0;i<prigr[n].length;i++) if( cipgr[n][i]<0) prigr[n][i]=1; else prigr[n][i]=priority[cipgr[n][i]];
00698               int p=0;
00699               while (p<prigr[n].length-1)
00700               {
00701                   if( prigr[n][p]<prigr[n][p+1]) {
00702                      int i=prigr[n][p]; prigr[n][p]=prigr[n][p+1]; prigr[n][p+1]=i;
00703                      if( p>0) p--;
00704                   } else p++;
00705               }
00706            }
00707            
00708            // divide each priority category into groups, then for each of these groups, split the contents out and reassign
00709            int groups[][]=SortAndGroup(priority);
00710            int nextpri=anyActualH ? 0 : 1;
00711            boolean repartitioned=false;
00712            
00713            for( int n=0;n<groups.length;n++)
00714            {
00715               // sort the groups according to their cipgr contents
00716               int p=0;
00717               while (p<groups[n].length-1)
00718               {
00719                   int i1=groups[n][p],i2=groups[n][p+1];
00720                   int cmp=0,sz=Math.max(prigr[i1].length,prigr[i2].length);
00721                   for( int i=0;i<sz;i++)
00722                   {
00723                      int v1=i<prigr[i1].length ? prigr[i1][i] : 0,v2=i<prigr[i2].length ? prigr[i2][i] : 0;
00724                      if( v1<v2) {cmp=-1; break;}
00725                      if( v1>v2) {cmp=1; break;}
00726                   }
00727                   if( cmp>0) {groups[n][p]=i2; groups[n][p+1]=i1; if( p>0) p--;}
00728                   else p++;
00729               }
00730               
00731               for( int i=0;i<groups[n].length;i++) 
00732               {
00733                   if( i==0) nextpri++;
00734                   else if( prigr[groups[n][i]].length!=prigr[groups[n][i-1]].length) {nextpri++; repartitioned=true;}
00735                   else
00736                   {
00737                      for( int j=0;j<prigr[groups[n][i]].length;j++) 
00738                          if( prigr[groups[n][i]][j]!=prigr[groups[n][i-1]][j]) {nextpri++; repartitioned=true; break;}
00739                   }
00740                   
00741                   priority[groups[n][i]]=nextpri;
00742               }
00743            }
00744 
00745            if( !repartitioned) break;
00746        }
00747     }
00748 
00749     // compute the chirality values for each atom centre
00750     void BuildChirality()
00751     {
00752        chiral=new int[NumAtoms()];
00753        
00754        boolean haswedge[]=new boolean[NumAtoms()];
00755        for( int n=0;n<NumAtoms();n++) haswedge[n]=false;
00756        for( int n=1;n<=NumBonds();n++) if( BondType(n)==BONDTYPE_INCLINED || BondType(n)==BONDTYPE_DECLINED) haswedge[BondFrom(n)-1]=true;
00757        
00758        int pri[]=new int[4];
00759        double x[]=new double[4],y[]=new double[4],z[]=new double[4];
00760        
00761        for( int n=0;n<NumAtoms();n++)
00762        {
00763            chiral[n]=STEREO_NONE;
00764            if( !(graph[n].length==4 || (graph[n].length==3 && AtomHydrogens(n+1)==1))) continue;
00765 
00766            for( int i=0;i<graph[n].length;i++)
00767            {
00768               pri[i]=priority[graph[n][i]];
00769               x[i]=AtomX(graph[n][i]+1)-AtomX(n+1);
00770               y[i]=AtomY(graph[n][i]+1)-AtomY(n+1);
00771               z[i]=0;
00772               if( haswedge[n]) 
00773                   for( int j=1;j<=NumBonds();j++) if( BondFrom(j)==n+1 && BondTo(j)==graph[n][i]+1)
00774               {
00775                   if( BondType(j)==BONDTYPE_INCLINED) z[i]=1;
00776                   if( BondType(j)==BONDTYPE_DECLINED) z[i]=-1;
00777                   break;
00778               }
00779            }
00780            if( graph[n].length==3) {pri[3]=0; x[3]=0; y[3]=0; z[3]=0;}
00781            
00782            int p=0;
00783            while (p<3) 
00784            {
00785               if( pri[p]>pri[p+1])
00786               {
00787                   int i; double d;
00788                   i=pri[p]; pri[p]=pri[p+1]; pri[p+1]=i;
00789                   d=x[p]; x[p]=x[p+1]; x[p+1]=d;
00790                   d=y[p]; y[p]=y[p+1]; y[p+1]=d;
00791                   d=z[p]; z[p]=z[p+1]; z[p+1]=d;
00792                   if( p>0) p--;
00793               }
00794               else p++;
00795            }
00796            if( (pri[0]==0 && pri[1]==1) || pri[0]==pri[1] || pri[1]==pri[2] || pri[2]==pri[3]) continue; // no topological chirality
00797            
00798            chiral[n]=STEREO_UNKNOWN;
00799            if( z[0]==0 && z[1]==0 && z[2]==0 && z[3]==0) continue; // all atoms are in-plane
00800            
00801            boolean sane=true;
00802            for( int i=0;i<4;i++) if( pri[0]!=0)
00803            {
00804               double r=x[i]*x[i]+y[i]*y[i]+z[i]*z[i];
00805               if( r<0.01*0.01) {sane=false; break;}
00806               r=1/Math.sqrt(r);
00807               x[i]=x[i]*r; y[i]=y[i]*r; z[i]=z[i]*r;
00808            }
00809            if( !sane) continue;
00810            
00811            if( pri[0]==0) // build a position for the implicit H
00812            {
00813               x[0]=-(x[1]+x[2]+x[3]);
00814               y[0]=-(y[1]+y[2]+y[3]);
00815               z[0]=-(z[1]+z[2]+z[3]);
00816               double r=x[0]*x[0]+y[0]*y[0]+z[0]*z[0];
00817               if( r<0.01*0.01) {sane=false; break;}
00818               r=1/Math.sqrt(r);
00819               x[0]=x[0]*r; y[0]=y[0]*r; z[0]=z[0]*r;
00820            }
00821            if( !sane) continue;
00822            
00823            double R=0,S=0;
00824            for( int i=1;i<=6;i++)
00825            {
00826               int a=0,b=0;
00827               if      (i==1) {a=1; b=2;} else if( i==2) {a=2; b=3;} else if( i==3) {a=3; b=1;}
00828               else if( i==4) {a=2; b=1;} else if( i==5) {a=3; b=2;} else if( i==6) {a=1; b=3;}
00829               double xx=y[a]*z[b]-y[b]*z[a] - x[0];
00830               double yy=z[a]*x[b]-z[b]*x[a] - y[0];
00831               double zz=x[a]*y[b]-x[b]*y[a] - z[0];
00832               if( i<=3) R+=xx*xx+yy*yy+zz*zz; else S+=xx*xx+yy*yy+zz*zz;
00833            }
00834            chiral[n]=R>S ? STEREO_POS : STEREO_NEG;
00835        }
00836     }
00837     
00838     // computer the cis/trans stereochemistry for each bond
00839     void BuildCisTrans()
00840     {
00841        cistrans=new int[NumBonds()];
00842        int sf[]=new int[2],st[]=new int[2];
00843        
00844        for( int n=0;n<NumBonds();n++) 
00845        {
00846            cistrans[n]=STEREO_NONE;
00847            int bf=BondFrom(n+1)-1,bt=BondTo(n+1)-1;
00848            if( BondOrder(n+1)!=2 || graph[bf].length<=1 || graph[bt].length<= 1 || graph[bf].length>3 || graph[bt].length>3) continue;
00849 
00850            int nf=0,nt=0;
00851            for( int i=0;i<graph[bf].length;i++) if( graph[bf][i]!=bt) sf[nf++]=graph[bf][i];
00852            for( int i=0;i<graph[bt].length;i++) if( graph[bt][i]!=bf) st[nt++]=graph[bt][i];
00853 
00854            if( nf==1) {if( AtomHydrogens(bf+1)!=1 || priority[sf[0]]==1) continue;}
00855            else 
00856            {
00857               if( priority[sf[0]]==priority[sf[1]]) continue;
00858               if( priority[sf[0]]<priority[sf[1]]) {int i=sf[0]; sf[0]=sf[1]; sf[1]=i;}
00859            }
00860            
00861            if( nt==1) {if( AtomHydrogens(bt+1)!=1 || priority[st[0]]==1) continue;}
00862            else 
00863            {
00864               if( priority[st[0]]==priority[st[1]]) continue;
00865               if( priority[st[0]]<priority[st[1]]) {int i=st[0]; st[0]=st[1]; st[1]=i;}
00866            }
00867        
00868            cistrans[n]=STEREO_UNKNOWN;
00869 
00870            double xa=AtomX(bf+1),ya=AtomY(bf+1),xb=AtomX(bt+1),yb=AtomY(bt+1);
00871            double tha0=Math.atan2(yb-ya,xb-xa),thb0=Math.PI+tha0;
00872            double tha1=Math.atan2(AtomY(sf[0]+1)-ya,AtomX(sf[0]+1)-xa);
00873            double tha2=nf==2 ? Math.atan2(AtomY(sf[1]+1)-ya,AtomX(sf[1]+1)-xa) : ThetaObtuse(tha0,tha1);
00874            double thb1=Math.atan2(AtomY(st[0]+1)-yb,AtomX(st[0]+1)-xb);
00875            double thb2=nt==2 ? Math.atan2(AtomY(st[1]+1)-yb,AtomX(st[1]+1)-xb) : ThetaObtuse(thb0,thb1);
00876            tha1-=tha0; tha2-=tha0; thb1-=thb0; thb2-=thb0;
00877            tha1+=(tha1<-Math.PI ? 2*Math.PI : 0)+(tha1>Math.PI ? -2*Math.PI : 0);
00878            tha2+=(tha2<-Math.PI ? 2*Math.PI : 0)+(tha2>Math.PI ? -2*Math.PI : 0);
00879            thb1+=(thb1<-Math.PI ? 2*Math.PI : 0)+(thb1>Math.PI ? -2*Math.PI : 0);
00880            thb2+=(thb2<-Math.PI ? 2*Math.PI : 0)+(thb2>Math.PI ? -2*Math.PI : 0);
00881 
00882            final double SMALL=5*Math.PI/180; // basically colinear
00883            if( Math.abs(tha1)<SMALL || Math.abs(tha2)<SMALL || Math.abs(thb1)<SMALL || Math.abs(thb2)<SMALL) continue;
00884            if( Math.abs(tha1)>Math.PI-SMALL || Math.abs(tha2)>Math.PI-SMALL) continue;
00885            if( Math.abs(thb1)>Math.PI-SMALL || Math.abs(thb2)>Math.PI-SMALL) continue;
00886            tha1=Math.signum(tha1); tha2=Math.signum(tha2); thb1=Math.signum(thb1); thb2=Math.signum(thb2);
00887            if( tha1==tha2 || thb1==thb2) continue;
00888            if( tha1*thb1<0) cistrans[n]=STEREO_POS;
00889            if( tha1*thb1>0) cistrans[n]=STEREO_NEG;
00890        }
00891     }
00892 
00893     // generally useful function which takes a list of numbers and sorts them, then bins the unique values into sub-arrays
00894     // (note: inefficient implementation, but could be improved easily enough)
00895     public static int[][] SortAndGroup(int[] Val)
00896     {
00897        ArrayList<Integer> unique=new ArrayList<Integer>();
00898        for( int n=0;n<Val.length;n++) unique.add(Val[n]);
00899        Collections.sort(unique);
00900        int p=0;
00901        while (p<unique.size()-1)
00902        {
00903            if( unique.get(p)==unique.get(p+1)) unique.remove(p); else p++;
00904        }
00905 
00906        int ret[][]=new int[unique.size()][];
00907        
00908        for( int n=0;n<Val.length;n++)
00909        {
00910            int grp=unique.indexOf(Val[n]);
00911            int[] cnt=new int[ret[grp]==null ? 1 : ret[grp].length+1];
00912            for( int i=0;i<cnt.length-1;i++) cnt[i]=ret[grp][i];
00913            cnt[cnt.length-1]=n;
00914            ret[grp]=cnt;
00915        }
00916        
00917        return ret;
00918     }
00919 
00920     // update the ring-block-identifier for each atom
00921     void BuildRingID()
00922     {
00923        ringID=new int[NumAtoms()];
00924        if( NumAtoms()==0) return;
00925        boolean visited[]=new boolean[NumAtoms()];
00926        for( int n=0;n<NumAtoms();n++) 
00927        {
00928            ringID[n]=0;
00929            visited[n]=false;
00930        }
00931 
00932        int path[]=new int[NumAtoms()+1],plen=0,numVisited=0;
00933        boolean rewound=false;
00934        while (true)
00935        {
00936            int last,current;
00937 
00938            if( plen==0) // find an element of a new component to visit
00939            {
00940               last=-1;
00941               for( current=0;visited[current];current++) {}
00942            }
00943            else
00944            {
00945               last=path[plen-1];
00946               current=-1;
00947               for( int n=0;n<graph[last].length;n++) if( !visited[graph[last][n]]) {current=graph[last][n]; break;}
00948            }
00949            
00950            /*System.out.print("numVisited="+numVisited+" last="+last+" current="+current+" path=");
00951            for( int n=0;n<plen;n++) System.out.print(path[n]+" ");
00952            System.out.println();*/
00953            
00954            if( current>=0 && plen>=2) // path is at least 2 items long, so look for any not-previous visited neighbours
00955            {
00956               int back=path[plen-1];
00957               //System.out.println(" back="+back);
00958               for( int n=0;n<graph[current].length;n++) 
00959               {
00960                   int join=graph[current][n];
00961                   //System.out.println(" join="+join+" visited[join]="+visited[join]);
00962                   if( join!=back && visited[join])
00963                   {
00964                      //System.out.print(" path:");
00965 
00966                      path[plen]=current;
00967                      for( int i=plen;i==plen || path[i+1]!=join;i--)
00968                      {
00969                          //System.out.print(" "+path[i]);
00970 
00971                          int id=ringID[path[i]];
00972                          if( id==0) ringID[path[i]]=last;
00973                          else if( id!=last)
00974                          {
00975                             for( int j=0;j<NumAtoms();j++) if( ringID[j]==id) ringID[j]=last;
00976                          }
00977                      }
00978                      //System.out.println();
00979                   }
00980               }
00981            }
00982            if( current>=0)  // can mark the new one as visited
00983            {
00984               visited[current]=true;
00985               path[plen++]=current;
00986               numVisited++;
00987               rewound=false;
00988            }
00989            else // otherwise, found nothing and must rewind the path
00990            {
00991               plen--;
00992               rewound=true;
00993            }
00994            
00995            if( numVisited==NumAtoms()) break;
00996        }
00997        
00998        // the ring ID's are not necessarily consecutive, so reassign them to 0=none, 1..NBlocks
00999        int nextID=0;
01000        for( int i=0;i<NumAtoms();i++) if( ringID[i]>0)
01001        {
01002            nextID--;
01003            for( int j=NumAtoms()-1;j>=i;j--) if( ringID[j]==ringID[i]) ringID[j]=nextID;
01004        }
01005        for( int i=0;i<NumAtoms();i++) ringID[i]=-ringID[i];
01006     }
01007     
01008     // ring hunter: recursive step; finds, compares and collects
01009     void RecursiveRingFind(int[] Path,int PSize,int Capacity,int RBlk,ArrayList<Object> Rings)
01010     {
01011        // not enough atoms yet, so look for new possibilities
01012        if( PSize<Capacity)
01013        {
01014            int last=Path[PSize-1];
01015            for( int n=0;n<graph[last-1].length;n++)
01016            {
01017               int adj=graph[last-1][n]+1;
01018               if( ringID[adj-1]!=RBlk) continue;
01019               boolean fnd=false;
01020               for( int i=0;i<PSize;i++) if( Path[i]==adj) {fnd=true; break;}
01021               if( !fnd)
01022               {
01023                   int newPath[]=Path.clone();
01024                   newPath[PSize]=adj;
01025                   RecursiveRingFind(newPath,PSize+1,Capacity,RBlk,Rings);
01026               }
01027            }
01028            return;
01029        }
01030        
01031        // path is full, so make sure it eats its tail
01032        int last=Path[PSize-1];
01033        boolean fnd=false;
01034        for( int n=0;n<graph[last-1].length;n++) if( graph[last-1][n]+1==Path[0]) {fnd=true; break;}
01035        if( !fnd) return;
01036        
01037        // reorder the array then look for duplicates
01038        int first=0;
01039        for( int n=1;n<PSize;n++) if( Path[n]<Path[first]) first=n;
01040        int fm=(first-1+PSize)%PSize,fp=(first+1)%PSize;
01041        boolean flip=Path[fm]<Path[fp];
01042        if( first!=0 || flip)
01043        {
01044            int newPath[]=new int[PSize];
01045            for( int n=0;n<PSize;n++) newPath[n]=Path[(first+(flip ? PSize-n : n))%PSize];
01046            Path=newPath;
01047        }
01048               
01049        for( int n=0;n<Rings.size();n++)
01050        {
01051            int[] look=(int[])Rings.get(n);
01052            boolean same=true;
01053            for( int i=0;i<PSize;i++) if( look[i]!=Path[i]) {same=false; break;}
01054            if( same) return;
01055        }
01056        
01057        Rings.add(Path);
01058     }
01059     
01060     // returns the angle maximally equidistant from Th1 and Th2
01061     double ThetaObtuse(double Th1,double Th2)
01062     {
01063        double dth=Th2-Th1;
01064        while (dth<-Math.PI) dth+=2*Math.PI;
01065        while (dth>Math.PI) dth-=2*Math.PI;
01066        return dth>0 ? Th1-0.5*(2*Math.PI-dth) : Th1+0.5*(2*Math.PI+dth);
01067     }
01068 }