| 1 | 
< | 
/* | 
| 2 | 
< | 
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project | 
| 3 | 
< | 
 *  | 
| 4 | 
< | 
 * Contact: oopse@oopse.org | 
| 5 | 
< | 
 *  | 
| 6 | 
< | 
 * This program is free software; you can redistribute it and/or | 
| 7 | 
< | 
 * modify it under the terms of the GNU Lesser General Public License | 
| 8 | 
< | 
 * as published by the Free Software Foundation; either version 2.1 | 
| 9 | 
< | 
 * of the License, or (at your option) any later version. | 
| 10 | 
< | 
 * All we ask is that proper credit is given for our work, which includes | 
| 11 | 
< | 
 * - but is not limited to - adding the above copyright notice to the beginning | 
| 12 | 
< | 
 * of your source code files, and to any copyright notice that you may distribute | 
| 13 | 
< | 
 * with programs based on this work. | 
| 14 | 
< | 
 *  | 
| 15 | 
< | 
 * This program is distributed in the hope that it will be useful, | 
| 16 | 
< | 
 * but WITHOUT ANY WARRANTY; without even the implied warranty of | 
| 17 | 
< | 
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the | 
| 18 | 
< | 
 * GNU Lesser General Public License for more details. | 
| 19 | 
< | 
 *  | 
| 20 | 
< | 
 * You should have received a copy of the GNU Lesser General Public License | 
| 21 | 
< | 
 * along with this program; if not, write to the Free Software | 
| 22 | 
< | 
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA. | 
| 1 | 
> | 
 /* | 
| 2 | 
> | 
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. | 
| 3 | 
  | 
 * | 
| 4 | 
+ | 
 * The University of Notre Dame grants you ("Licensee") a | 
| 5 | 
+ | 
 * non-exclusive, royalty free, license to use, modify and | 
| 6 | 
+ | 
 * redistribute this software in source and binary code form, provided | 
| 7 | 
+ | 
 * that the following conditions are met: | 
| 8 | 
+ | 
 * | 
| 9 | 
+ | 
 * 1. Acknowledgement of the program authors must be made in any | 
| 10 | 
+ | 
 *    publication of scientific results based in part on use of the | 
| 11 | 
+ | 
 *    program.  An acceptable form of acknowledgement is citation of | 
| 12 | 
+ | 
 *    the article in which the program was described (Matthew | 
| 13 | 
+ | 
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher | 
| 14 | 
+ | 
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented | 
| 15 | 
+ | 
 *    Parallel Simulation Engine for Molecular Dynamics," | 
| 16 | 
+ | 
 *    J. Comput. Chem. 26, pp. 252-271 (2005)) | 
| 17 | 
+ | 
 * | 
| 18 | 
+ | 
 * 2. Redistributions of source code must retain the above copyright | 
| 19 | 
+ | 
 *    notice, this list of conditions and the following disclaimer. | 
| 20 | 
+ | 
 * | 
| 21 | 
+ | 
 * 3. Redistributions in binary form must reproduce the above copyright | 
| 22 | 
+ | 
 *    notice, this list of conditions and the following disclaimer in the | 
| 23 | 
+ | 
 *    documentation and/or other materials provided with the | 
| 24 | 
+ | 
 *    distribution. | 
| 25 | 
+ | 
 * | 
| 26 | 
+ | 
 * This software is provided "AS IS," without a warranty of any | 
| 27 | 
+ | 
 * kind. All express or implied conditions, representations and | 
| 28 | 
+ | 
 * warranties, including any implied warranty of merchantability, | 
| 29 | 
+ | 
 * fitness for a particular purpose or non-infringement, are hereby | 
| 30 | 
+ | 
 * excluded.  The University of Notre Dame and its licensors shall not | 
| 31 | 
+ | 
 * be liable for any damages suffered by licensee as a result of | 
| 32 | 
+ | 
 * using, modifying or distributing the software or its | 
| 33 | 
+ | 
 * derivatives. In no event will the University of Notre Dame or its | 
| 34 | 
+ | 
 * licensors be liable for any lost revenue, profit or data, or for | 
| 35 | 
+ | 
 * direct, indirect, special, consequential, incidental or punitive | 
| 36 | 
+ | 
 * damages, however caused and regardless of the theory of liability, | 
| 37 | 
+ | 
 * arising out of the use of or inability to use software, even if the | 
| 38 | 
+ | 
 * University of Notre Dame has been advised of the possibility of | 
| 39 | 
+ | 
 * such damages. | 
| 40 | 
  | 
 */ | 
| 41 | 
< | 
 | 
| 41 | 
> | 
  | 
| 42 | 
  | 
/** | 
| 43 | 
  | 
 * @file SquareMatrix3.hpp | 
| 44 | 
  | 
 * @author Teng Lin | 
| 57 | 
  | 
    template<typename Real> | 
| 58 | 
  | 
    class SquareMatrix3 : public SquareMatrix<Real, 3> { | 
| 59 | 
  | 
        public: | 
| 60 | 
+ | 
 | 
| 61 | 
+ | 
            typedef Real ElemType; | 
| 62 | 
+ | 
            typedef Real* ElemPoinerType; | 
| 63 | 
  | 
             | 
| 64 | 
  | 
            /** default constructor */ | 
| 65 | 
  | 
            SquareMatrix3() : SquareMatrix<Real, 3>() { | 
| 66 | 
  | 
            } | 
| 67 | 
  | 
 | 
| 68 | 
+ | 
            /** Constructs and initializes every element of this matrix to a scalar */  | 
| 69 | 
+ | 
            SquareMatrix3(Real s) : SquareMatrix<Real,3>(s){ | 
| 70 | 
+ | 
            } | 
| 71 | 
+ | 
 | 
| 72 | 
+ | 
            /** Constructs and initializes from an array */  | 
| 73 | 
+ | 
            SquareMatrix3(Real* array) : SquareMatrix<Real,3>(array){ | 
| 74 | 
+ | 
            } | 
| 75 | 
+ | 
 | 
| 76 | 
+ | 
 | 
| 77 | 
  | 
            /** copy  constructor */ | 
| 78 | 
  | 
            SquareMatrix3(const SquareMatrix<Real, 3>& m)  : SquareMatrix<Real, 3>(m) { | 
| 79 | 
  | 
            } | 
| 80 | 
< | 
 | 
| 80 | 
> | 
             | 
| 81 | 
  | 
            SquareMatrix3( const Vector3<Real>& eulerAngles) { | 
| 82 | 
  | 
                setupRotMat(eulerAngles); | 
| 83 | 
  | 
            } | 
| 101 | 
  | 
                    return *this; | 
| 102 | 
  | 
                 SquareMatrix<Real, 3>::operator=(m); | 
| 103 | 
  | 
                 return *this; | 
| 104 | 
+ | 
            } | 
| 105 | 
+ | 
 | 
| 106 | 
+ | 
 | 
| 107 | 
+ | 
            SquareMatrix3<Real>& operator =(const Quaternion<Real>& q) { | 
| 108 | 
+ | 
                this->setupRotMat(q); | 
| 109 | 
+ | 
                return *this; | 
| 110 | 
  | 
            } | 
| 111 | 
  | 
 | 
| 112 | 
  | 
            /** | 
| 231 | 
  | 
            */             | 
| 232 | 
  | 
            Vector3<Real> toEulerAngles() { | 
| 233 | 
  | 
                Vector3<Real> myEuler; | 
| 234 | 
< | 
                Real phi,theta,psi,eps; | 
| 235 | 
< | 
                Real ctheta,stheta; | 
| 234 | 
> | 
                Real phi; | 
| 235 | 
> | 
                Real theta; | 
| 236 | 
> | 
                Real psi; | 
| 237 | 
> | 
                Real ctheta; | 
| 238 | 
> | 
                Real stheta; | 
| 239 | 
  | 
                 | 
| 240 | 
  | 
                // set the tolerance for Euler angles and rotation elements | 
| 241 | 
  | 
 | 
| 285 | 
  | 
 | 
| 286 | 
  | 
                return(x + y + z); | 
| 287 | 
  | 
            }             | 
| 288 | 
+ | 
 | 
| 289 | 
+ | 
            /** Returns the trace of this matrix. */ | 
| 290 | 
+ | 
            Real trace() const { | 
| 291 | 
+ | 
                return data_[0][0] + data_[1][1] + data_[2][2]; | 
| 292 | 
+ | 
            } | 
| 293 | 
  | 
             | 
| 294 | 
  | 
            /** | 
| 295 | 
  | 
             * Sets the value of this matrix to  the inversion of itself.  | 
| 296 | 
  | 
             * @note since simple algorithm can be applied to inverse the 3 by 3 matrix, we hide the  | 
| 297 | 
  | 
             * implementation of inverse in SquareMatrix class | 
| 298 | 
  | 
             */ | 
| 299 | 
< | 
            SquareMatrix3<Real>  inverse() { | 
| 299 | 
> | 
            SquareMatrix3<Real>  inverse() const { | 
| 300 | 
  | 
                SquareMatrix3<Real> m; | 
| 301 | 
  | 
                double det = determinant(); | 
| 302 | 
  | 
                if (fabs(det) <= oopse::epsilon) { | 
| 465 | 
  | 
        // transpose the eigenvectors back again | 
| 466 | 
  | 
        v = v.transpose(); | 
| 467 | 
  | 
        return ; | 
| 468 | 
+ | 
    } | 
| 469 | 
+ | 
 | 
| 470 | 
+ | 
    /** | 
| 471 | 
+ | 
    * Return the multiplication of two matrixes  (m1 * m2).  | 
| 472 | 
+ | 
    * @return the multiplication of two matrixes | 
| 473 | 
+ | 
    * @param m1 the first matrix | 
| 474 | 
+ | 
    * @param m2 the second matrix | 
| 475 | 
+ | 
    */ | 
| 476 | 
+ | 
    template<typename Real>  | 
| 477 | 
+ | 
    inline SquareMatrix3<Real> operator *(const SquareMatrix3<Real>& m1, const SquareMatrix3<Real>& m2) { | 
| 478 | 
+ | 
        SquareMatrix3<Real> result; | 
| 479 | 
+ | 
 | 
| 480 | 
+ | 
            for (unsigned int i = 0; i < 3; i++) | 
| 481 | 
+ | 
                for (unsigned int j = 0; j < 3; j++) | 
| 482 | 
+ | 
                    for (unsigned int k = 0; k < 3; k++) | 
| 483 | 
+ | 
                        result(i, j)  += m1(i, k) * m2(k, j);                 | 
| 484 | 
+ | 
 | 
| 485 | 
+ | 
        return result; | 
| 486 | 
  | 
    } | 
| 487 | 
+ | 
 | 
| 488 | 
+ | 
    template<typename Real>  | 
| 489 | 
+ | 
    inline SquareMatrix3<Real> outProduct(const Vector3<Real>& v1, const Vector3<Real>& v2) { | 
| 490 | 
+ | 
        SquareMatrix3<Real> result; | 
| 491 | 
+ | 
 | 
| 492 | 
+ | 
            for (unsigned int i = 0; i < 3; i++) { | 
| 493 | 
+ | 
                for (unsigned int j = 0; j < 3; j++) { | 
| 494 | 
+ | 
                        result(i, j)  = v1[i] * v2[j];                 | 
| 495 | 
+ | 
                } | 
| 496 | 
+ | 
            } | 
| 497 | 
+ | 
             | 
| 498 | 
+ | 
        return result;         | 
| 499 | 
+ | 
    } | 
| 500 | 
+ | 
 | 
| 501 | 
+ | 
     | 
| 502 | 
  | 
    typedef SquareMatrix3<double> Mat3x3d; | 
| 503 | 
  | 
    typedef SquareMatrix3<double> RotMat3x3d; | 
| 504 | 
  | 
 |