javax.vecmath/src/javax/vecmath/Matrix4d.java

   1 /*
   2  * $RCSfile: Matrix4d.java,v $
   3  *
   4  * Copyright 1996-2008 Sun Microsystems, Inc.  All Rights Reserved.
   5  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
   6  *
   7  * This code is free software; you can redistribute it and/or modify it
   8  * under the terms of the GNU General Public License version 2 only, as
   9  * published by the Free Software Foundation.  Sun designates this
  10  * particular file as subject to the "Classpath" exception as provided
  11  * by Sun in the LICENSE file that accompanied this code.
  12  *
  13  * This code is distributed in the hope that it will be useful, but WITHOUT
  14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16  * version 2 for more details (a copy is included in the LICENSE file that
  17  * accompanied this code).
  18  *
  19  * You should have received a copy of the GNU General Public License version
  20  * 2 along with this work; if not, write to the Free Software Foundation,
  21  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  22  *
  23  * Please contact Sun Microsystems, Inc., 4150 Network Circle, Santa Clara,
  24  * CA 95054 USA or visit www.sun.com if you need additional information or
  25  * have any questions.
  26  *
  27  * $Revision: 1.8 $
  28  * $Date: 2008/02/28 20:18:50 $
  29  * $State: Exp $
  30  */
  31
  32 package javax.vecmath;
  33
  34 import java.lang.Math;
  35
  36 /**
  37  * A double precision floating point 4 by 4 matrix.
  38  * Primarily to support 3D rotations.
  39  *
  40  */
  41 public class Matrix4d implements java.io.Serializable, Cloneable {
  42
  43     // Compatible with 1.1
  44     static final long serialVersionUID = 8223903484171633710L;
  45
  46     /**
  47      *  The first element of the first row.
  48      */
  49     public      double  m00;
  50
  51     /**
  52      *  The second element of the first row.
  53      */
  54     public      double  m01;
  55
  56     /**
  57      *  The third element of the first row.
  58      */
  59     public      double  m02;
  60
  61     /**
  62      *  The fourth element of the first row.
  63      */
  64     public      double  m03;
  65
  66     /**
  67      *  The first element of the second row.
  68      */
  69     public      double  m10;
  70
  71     /**
  72      *  The second element of the second row.
  73      */
  74     public      double  m11;
  75
  76     /**
  77      *  The third element of the second row.
  78      */
  79     public      double  m12;
  80
  81     /**
  82      *  The fourth element of the second row.
  83      */
  84     public      double  m13;
  85
  86     /**
  87      *  The first element of the third row.
  88      */
  89     public      double  m20;
  90
  91     /**
  92      *  The second element of the third row.
  93      */
  94     public      double  m21;
  95
  96     /**
  97      *  The third element of the third row.
  98      */
  99     public      double  m22;
 100
 101     /**
 102      *  The fourth element of the third row.
 103      */
 104     public      double  m23;
 105
 106     /**
 107      *  The first element of the fourth row.
 108      */
 109     public      double  m30;
 110
 111     /**
 112      *  The second element of the fourth row.
 113      */
 114     public      double  m31;
 115
 116     /**
 117      *  The third element of the fourth row.
 118      */
 119     public      double  m32;
 120
 121     /**
 122      *  The fourth element of the fourth row.
 123      */
 124     public      double  m33;
 125     /*
 126     double[] tmp = new double[16];
 127     double[] tmp_rot = new double[9];  // scratch matrix
 128     double[] tmp_scale = new double[3];  // scratch matrix
 129     */
 130     private static final double EPS = 1.0E-10;
 131
 132
 133     /**
 134      * Constructs and initializes a Matrix4d from the specified 16 values.
 135      * @param m00 the [0][0] element
 136      * @param m01 the [0][1] element
 137      * @param m02 the [0][2] element
 138      * @param m03 the [0][3] element
 139      * @param m10 the [1][0] element
 140      * @param m11 the [1][1] element
 141      * @param m12 the [1][2] element
 142      * @param m13 the [1][3] element
 143      * @param m20 the [2][0] element
 144      * @param m21 the [2][1] element
 145      * @param m22 the [2][2] element
 146      * @param m23 the [2][3] element
 147      * @param m30 the [3][0] element
 148      * @param m31 the [3][1] element
 149      * @param m32 the [3][2] element
 150      * @param m33 the [3][3] element
 151      */
 152     public Matrix4d(double m00, double m01, double m02, double m03,
 153                     double m10, double m11, double m12, double m13,
 154                     double m20, double m21, double m22, double m23,
 155                     double m30, double m31, double m32, double m33)
 156     {
 157         this.m00 = m00;
 158         this.m01 = m01;
 159         this.m02 = m02;
 160         this.m03 = m03;
 161
 162         this.m10 = m10;
 163         this.m11 = m11;
 164         this.m12 = m12;
 165         this.m13 = m13;
 166
 167         this.m20 = m20;
 168         this.m21 = m21;
 169         this.m22 = m22;
 170         this.m23 = m23;
 171
 172         this.m30 = m30;
 173         this.m31 = m31;
 174         this.m32 = m32;
 175         this.m33 = m33;
 176
 177     }
 178
 179     /**
 180      * Constructs and initializes a Matrix4d from the specified 16
 181      * element array.  this.m00 =v[0], this.m01=v[1], etc.
 182      * @param v the array of length 16 containing in order
 183      */
 184     public Matrix4d(double[] v)
 185     {
 186         this.m00 = v[ 0];
 187         this.m01 = v[ 1];
 188         this.m02 = v[ 2];
 189         this.m03 = v[ 3];
 190
 191         this.m10 = v[ 4];
 192         this.m11 = v[ 5];
 193         this.m12 = v[ 6];
 194         this.m13 = v[ 7];
 195
 196         this.m20 = v[ 8];
 197         this.m21 = v[ 9];
 198         this.m22 = v[10];
 199         this.m23 = v[11];
 200
 201         this.m30 = v[12];
 202         this.m31 = v[13];
 203         this.m32 = v[14];
 204         this.m33 = v[15];
 205
 206     }
 207
 208    /**
 209      * Constructs and initializes a Matrix4d from the quaternion,
 210      * translation, and scale values; the scale is applied only to the
 211      * rotational components of the matrix (upper 3x3) and not to the
 212      * translational components.
 213      * @param q1  the quaternion value representing the rotational component
 214      * @param t1  the translational component of the matrix
 215      * @param s   the scale value applied to the rotational components
 216      */
 217     public Matrix4d(Quat4d q1, Vector3d t1, double s)
 218     {
 219         m00 = s*(1.0 - 2.0*q1.y*q1.y - 2.0*q1.z*q1.z);
 220         m10 = s*(2.0*(q1.x*q1.y + q1.w*q1.z));
 221         m20 = s*(2.0*(q1.x*q1.z - q1.w*q1.y));
 222
 223         m01 = s*(2.0*(q1.x*q1.y - q1.w*q1.z));
 224         m11 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.z*q1.z);
 225         m21 = s*(2.0*(q1.y*q1.z + q1.w*q1.x));
 226
 227         m02 = s*(2.0*(q1.x*q1.z + q1.w*q1.y));
 228         m12 = s*(2.0*(q1.y*q1.z - q1.w*q1.x));
 229         m22 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.y*q1.y);
 230
 231         m03 = t1.x;
 232         m13 = t1.y;
 233         m23 = t1.z;
 234
 235         m30 = 0.0;
 236         m31 = 0.0;
 237         m32 = 0.0;
 238         m33 = 1.0;
 239
 240     }
 241
 242    /**
 243      * Constructs and initializes a Matrix4d from the quaternion,
 244      * translation, and scale values; the scale is applied only to the
 245      * rotational components of the matrix (upper 3x3) and not to the
 246      * translational components.
 247      * @param q1  the quaternion value representing the rotational component
 248      * @param t1  the translational component of the matrix
 249      * @param s   the scale value applied to the rotational components
 250      */
 251     public Matrix4d(Quat4f q1, Vector3d t1, double s)
 252     {
 253         m00 = s*(1.0 - 2.0*q1.y*q1.y - 2.0*q1.z*q1.z);
 254         m10 = s*(2.0*(q1.x*q1.y + q1.w*q1.z));
 255         m20 = s*(2.0*(q1.x*q1.z - q1.w*q1.y));
 256
 257         m01 = s*(2.0*(q1.x*q1.y - q1.w*q1.z));
 258         m11 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.z*q1.z);
 259         m21 = s*(2.0*(q1.y*q1.z + q1.w*q1.x));
 260
 261         m02 = s*(2.0*(q1.x*q1.z + q1.w*q1.y));
 262         m12 = s*(2.0*(q1.y*q1.z - q1.w*q1.x));
 263         m22 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.y*q1.y);
 264
 265         m03 = t1.x;
 266         m13 = t1.y;
 267         m23 = t1.z;
 268
 269         m30 = 0.0;
 270         m31 = 0.0;
 271         m32 = 0.0;
 272         m33 = 1.0;
 273
 274     }
 275
 276    /**
 277      *  Constructs a new matrix with the same values as the
 278      *  Matrix4d parameter.
 279      *  @param m1  the source matrix
 280      */
 281    public Matrix4d(Matrix4d m1)
 282    {
 283         this.m00 = m1.m00;
 284         this.m01 = m1.m01;
 285         this.m02 = m1.m02;
 286         this.m03 = m1.m03;
 287
 288         this.m10 = m1.m10;
 289         this.m11 = m1.m11;
 290         this.m12 = m1.m12;
 291         this.m13 = m1.m13;
 292
 293         this.m20 = m1.m20;
 294         this.m21 = m1.m21;
 295         this.m22 = m1.m22;
 296         this.m23 = m1.m23;
 297
 298         this.m30 = m1.m30;
 299         this.m31 = m1.m31;
 300         this.m32 = m1.m32;
 301         this.m33 = m1.m33;
 302
 303    }
 304
 305    /**
 306      *  Constructs a new matrix with the same values as the
 307      *  Matrix4f parameter.
 308      *  @param m1  the source matrix
 309      */
 310    public Matrix4d(Matrix4f m1)
 311    {
 312         this.m00 = m1.m00;
 313         this.m01 = m1.m01;
 314         this.m02 = m1.m02;
 315         this.m03 = m1.m03;
 316
 317         this.m10 = m1.m10;
 318         this.m11 = m1.m11;
 319         this.m12 = m1.m12;
 320         this.m13 = m1.m13;
 321
 322         this.m20 = m1.m20;
 323         this.m21 = m1.m21;
 324         this.m22 = m1.m22;
 325         this.m23 = m1.m23;
 326
 327         this.m30 = m1.m30;
 328         this.m31 = m1.m31;
 329         this.m32 = m1.m32;
 330         this.m33 = m1.m33;
 331
 332    }
 333
 334    /**
 335      * Constructs and initializes a Matrix4d from the rotation matrix,
 336      * translation, and scale values; the scale is applied only to the
 337      * rotational components of the matrix (upper 3x3) and not to the
 338      * translational components of the matrix.
 339      * @param m1  the rotation matrix representing the rotational components
 340      * @param t1  the translational components of the matrix
 341      * @param s   the scale value applied to the rotational components
 342      */
 343     public Matrix4d(Matrix3f m1, Vector3d t1, double s)
 344     {
 345         this.m00 = m1.m00*s;
 346         this.m01 = m1.m01*s;
 347         this.m02 = m1.m02*s;
 348         this.m03 = t1.x;
 349
 350         this.m10 = m1.m10*s;
 351         this.m11 = m1.m11*s;
 352         this.m12 = m1.m12*s;
 353         this.m13 = t1.y;
 354
 355         this.m20 = m1.m20*s;
 356         this.m21 = m1.m21*s;
 357         this.m22 = m1.m22*s;
 358         this.m23 = t1.z;
 359
 360         this.m30 = 0.0;
 361         this.m31 = 0.0;
 362         this.m32 = 0.0;
 363         this.m33 = 1.0;
 364
 365     }
 366
 367    /**
 368      * Constructs and initializes a Matrix4f from the rotation matrix,
 369      * translation, and scale values; the scale is applied only to the
 370      * rotational components of the matrix (upper 3x3) and not to the
 371      * translational components of the matrix.
 372      * @param m1  the rotation matrix representing the rotational components
 373      * @param t1  the translational components of the matrix
 374      * @param s   the scale value applied to the rotational components
 375      */
 376     public Matrix4d(Matrix3d m1, Vector3d t1, double s)
 377     {
 378         this.m00 = m1.m00*s;
 379         this.m01 = m1.m01*s;
 380         this.m02 = m1.m02*s;
 381         this.m03 = t1.x;
 382
 383         this.m10 = m1.m10*s;
 384         this.m11 = m1.m11*s;
 385         this.m12 = m1.m12*s;
 386         this.m13 = t1.y;
 387
 388         this.m20 = m1.m20*s;
 389         this.m21 = m1.m21*s;
 390         this.m22 = m1.m22*s;
 391         this.m23 = t1.z;
 392
 393         this.m30 = 0.0;
 394         this.m31 = 0.0;
 395         this.m32 = 0.0;
 396         this.m33 = 1.0;
 397
 398     }
 399
 400     /**
 401      * Constructs and initializes a Matrix4d to all zeros.
 402      */
 403     public Matrix4d()
 404     {
 405         this.m00 = 0.0;
 406         this.m01 = 0.0;
 407         this.m02 = 0.0;
 408         this.m03 = 0.0;
 409
 410         this.m10 = 0.0;
 411         this.m11 = 0.0;
 412         this.m12 = 0.0;
 413         this.m13 = 0.0;
 414
 415         this.m20 = 0.0;
 416         this.m21 = 0.0;
 417         this.m22 = 0.0;
 418         this.m23 = 0.0;
 419
 420         this.m30 = 0.0;
 421         this.m31 = 0.0;
 422         this.m32 = 0.0;
 423         this.m33 = 0.0;
 424
 425     }
 426
 427    /**
 428      * Returns a string that contains the values of this Matrix4d.
 429      * @return the String representation
 430      */
 431     public String toString() {
 432       return
 433         this.m00 + ", " + this.m01 + ", " + this.m02 + ", " + this.m03 + "\n" +
 434         this.m10 + ", " + this.m11 + ", " + this.m12 + ", " + this.m13 + "\n" +
 435         this.m20 + ", " + this.m21 + ", " + this.m22 + ", " + this.m23 + "\n" +
 436         this.m30 + ", " + this.m31 + ", " + this.m32 + ", " + this.m33 + "\n";
 437     }
 438
 439     /**
 440      * Sets this Matrix4d to identity.
 441      */
 442     public final void setIdentity()
 443     {
 444         this.m00 = 1.0;
 445         this.m01 = 0.0;
 446         this.m02 = 0.0;
 447         this.m03 = 0.0;
 448
 449         this.m10 = 0.0;
 450         this.m11 = 1.0;
 451         this.m12 = 0.0;
 452         this.m13 = 0.0;
 453
 454         this.m20 = 0.0;
 455         this.m21 = 0.0;
 456         this.m22 = 1.0;
 457         this.m23 = 0.0;
 458
 459         this.m30 = 0.0;
 460         this.m31 = 0.0;
 461         this.m32 = 0.0;
 462         this.m33 = 1.0;
 463     }
 464
 465     /**
 466      * Sets the specified element of this matrix4f to the value provided.
 467      * @param row the row number to be modified (zero indexed)
 468      * @param column the column number to be modified (zero indexed)
 469      * @param value the new value
 470      */
 471     public final void setElement(int row, int column, double value)
 472     {
 473         switch (row)
 474           {
 475           case 0:
 476             switch(column)
 477               {
 478               case 0:
 479                 this.m00 = value;
 480                 break;
 481               case 1:
 482                 this.m01 = value;
 483                 break;
 484               case 2:
 485                 this.m02 = value;
 486                 break;
 487               case 3:
 488                 this.m03 = value;
 489                 break;
 490               default:
 491                 throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d0"));
 492               }
 493             break;
 494
 495           case 1:
 496             switch(column)
 497               {
 498               case 0:
 499                 this.m10 = value;
 500                 break;
 501               case 1:
 502                 this.m11 = value;
 503                 break;
 504               case 2:
 505                 this.m12 = value;
 506                 break;
 507               case 3:
 508                 this.m13 = value;
 509                 break;
 510               default:
 511                 throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d0"));
 512               }
 513             break;
 514
 515           case 2:
 516             switch(column)
 517               {
 518               case 0:
 519                 this.m20 = value;
 520                 break;
 521               case 1:
 522                 this.m21 = value;
 523                 break;
 524               case 2:
 525                 this.m22 = value;
 526                 break;
 527               case 3:
 528                 this.m23 = value;
 529                 break;
 530               default:
 531                 throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d0"));
 532               }
 533             break;
 534
 535           case 3:
 536             switch(column)
 537               {
 538               case 0:
 539                 this.m30 = value;
 540                 break;
 541               case 1:
 542                 this.m31 = value;
 543                 break;
 544               case 2:
 545                 this.m32 = value;
 546                 break;
 547               case 3:
 548                 this.m33 = value;
 549                 break;
 550               default:
 551                 throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d0"));
 552               }
 553             break;
 554
 555           default:
 556                 throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d0"));
 557           }
 558     }
 559
 560     /**
 561      * Retrieves the value at the specified row and column of this matrix.
 562      * @param row the row number to be retrieved (zero indexed)
 563      * @param column the column number to be retrieved (zero indexed)
 564      * @return the value at the indexed element
 565      */
 566     public final double getElement(int row, int column)
 567     {
 568         switch (row)
 569           {
 570           case 0:
 571             switch(column)
 572               {
 573               case 0:
 574                 return(this.m00);
 575               case 1:
 576                 return(this.m01);
 577               case 2:
 578                 return(this.m02);
 579               case 3:
 580                 return(this.m03);
 581               default:
 582                 break;
 583               }
 584             break;
 585           case 1:
 586             switch(column)
 587               {
 588               case 0:
 589                 return(this.m10);
 590               case 1:
 591                 return(this.m11);
 592               case 2:
 593                 return(this.m12);
 594               case 3:
 595                 return(this.m13);
 596               default:
 597                 break;
 598               }
 599             break;
 600
 601           case 2:
 602             switch(column)
 603               {
 604               case 0:
 605                 return(this.m20);
 606               case 1:
 607                 return(this.m21);
 608               case 2:
 609                 return(this.m22);
 610               case 3:
 611                 return(this.m23);
 612               default:
 613                 break;
 614               }
 615             break;
 616
 617           case 3:
 618             switch(column)
 619               {
 620               case 0:
 621                 return(this.m30);
 622               case 1:
 623                 return(this.m31);
 624               case 2:
 625                 return(this.m32);
 626               case 3:
 627                 return(this.m33);
 628               default:
 629                 break;
 630               }
 631             break;
 632
 633           default:
 634             break;
 635           }
 636         throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d1"));
 637     }
 638
 639     /**
 640      * Copies the matrix values in the specified row into the vector parameter.
 641      * @param row  the matrix row
 642      * @param v    the vector into which the matrix row values will be copied
 643      */
 644     public final void getRow(int row, Vector4d v) {
 645         if( row == 0 ) {
 646            v.x = m00;
 647            v.y = m01;
 648            v.z = m02;
 649            v.w = m03;
 650         } else if(row == 1) {
 651            v.x = m10;
 652            v.y = m11;
 653            v.z = m12;
 654            v.w = m13;
 655         } else if(row == 2) {
 656            v.x = m20;
 657            v.y = m21;
 658            v.z = m22;
 659            v.w = m23;
 660         } else if(row == 3) {
 661            v.x = m30;
 662            v.y = m31;
 663            v.z = m32;
 664            v.w = m33;
 665         } else {
 666           throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d2"));
 667         }
 668     }
 669
 670
 671     /**
 672      * Copies the matrix values in the specified row into the array parameter.
 673      * @param row  the matrix row
 674      * @param v    the array into which the matrix row values will be copied
 675      */
 676     public final void getRow(int row, double v[]) {
 677         if( row == 0 ) {
 678                 v[0] = m00;
 679                 v[1] = m01;
 680                 v[2] = m02;
 681                 v[3] = m03;
 682         } else if(row == 1) {
 683                 v[0] = m10;
 684                 v[1] = m11;
 685                 v[2] = m12;
 686                 v[3] = m13;
 687         } else if(row == 2) {
 688                 v[0] = m20;
 689                 v[1] = m21;
 690                 v[2] = m22;
 691                 v[3] = m23;
 692         } else if(row == 3) {
 693                 v[0] = m30;
 694                 v[1] = m31;
 695                 v[2] = m32;
 696                 v[3] = m33;
 697         } else {
 698
 699           throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d2"));
 700         }
 701     }
 702
 703
 704
 705     /**
 706      * Copies the matrix values in the specified column into the vector
 707      * parameter.
 708      * @param column  the matrix column
 709      * @param v    the vector into which the matrix column values will be copied
 710      */
 711     public final void getColumn(int column, Vector4d v) {
 712         if( column == 0 ) {
 713            v.x = m00;
 714            v.y = m10;
 715            v.z = m20;
 716            v.w = m30;
 717         } else if(column == 1) {
 718            v.x = m01;
 719            v.y = m11;
 720            v.z = m21;
 721            v.w = m31;
 722         } else if(column == 2) {
 723            v.x = m02;
 724            v.y = m12;
 725            v.z = m22;
 726            v.w = m32;
 727         } else if(column == 3) {
 728            v.x = m03;
 729            v.y = m13;
 730            v.z = m23;
 731            v.w = m33;
 732         } else {
 733           throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d3"));
 734
 735         }
 736
 737     }
 738
 739
 740
 741     /**
 742      * Copies the matrix values in the specified column into the array
 743      * parameter.
 744      * @param column the matrix column
 745      * @param v    the array into which the matrix column values will be copied
 746      */
 747     public final void getColumn(int column, double v[]) {
 748         if( column == 0 ) {
 749            v[0] = m00;
 750            v[1] = m10;
 751            v[2] = m20;
 752            v[3] = m30;
 753         } else if(column == 1) {
 754            v[0] = m01;
 755            v[1] = m11;
 756            v[2] = m21;
 757            v[3] = m31;
 758         } else if(column == 2) {
 759            v[0] = m02;
 760            v[1] = m12;
 761            v[2] = m22;
 762            v[3] = m32;
 763         } else if(column == 3) {
 764            v[0] = m03;
 765            v[1] = m13;
 766            v[2] = m23;
 767            v[3] = m33;
 768         } else {
 769           throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d3"));
 770
 771         }
 772
 773     }
 774
 775
 776     /**
 777      * Performs an SVD normalization of this matrix in order to acquire
 778      * the normalized rotational component; the values are placed into
 779      * the Matrix3d parameter.
 780      * @param m1 the matrix into which the rotational component is placed
 781      */
 782     public final void get(Matrix3d m1)
 783     {
 784
 785         double[]    tmp_rot = new double[9];  // scratch matrix
 786         double[]    tmp_scale = new double[3];  // scratch matrix
 787         getScaleRotate( tmp_scale, tmp_rot );
 788
 789         m1.m00 = tmp_rot[0];
 790         m1.m01 = tmp_rot[1];
 791         m1.m02 = tmp_rot[2];
 792
 793         m1.m10 = tmp_rot[3];
 794         m1.m11 = tmp_rot[4];
 795         m1.m12 = tmp_rot[5];
 796
 797         m1.m20 = tmp_rot[6];
 798         m1.m21 = tmp_rot[7];
 799         m1.m22 = tmp_rot[8];
 800
 801     }
 802
 803
 804     /**
 805      * Performs an SVD normalization of this matrix in order to acquire
 806      * the normalized rotational component; the values are placed into
 807      * the Matrix3f parameter.
 808      * @param m1 the matrix into which the rotational component is placed
 809      */
 810     public final void get(Matrix3f m1)
 811     {
 812         double[]    tmp_rot = new double[9];  // scratch matrix
 813         double[]    tmp_scale = new double[3];  // scratch matrix
 814
 815         getScaleRotate( tmp_scale, tmp_rot );
 816
 817         m1.m00 = (float)tmp_rot[0];
 818         m1.m01 = (float)tmp_rot[1];
 819         m1.m02 = (float)tmp_rot[2];
 820
 821         m1.m10 = (float)tmp_rot[3];
 822         m1.m11 = (float)tmp_rot[4];
 823         m1.m12 = (float)tmp_rot[5];
 824
 825         m1.m20 = (float)tmp_rot[6];
 826         m1.m21 = (float)tmp_rot[7];
 827         m1.m22 = (float)tmp_rot[8];
 828     }
 829
 830    /**
 831      * Performs an SVD normalization of this matrix to calculate
 832      * the rotation as a 3x3 matrix, the translation, and the scale.
 833      * None of the matrix values are modified.
 834      * @param m1  the normalized matrix representing the rotation
 835      * @param t1  the translation component
 836      * @return  the scale component of this transform
 837      */
 838     public final double get(Matrix3d m1, Vector3d t1)
 839     {
 840
 841         double[]    tmp_rot = new double[9];  // scratch matrix
 842         double[]    tmp_scale = new double[3];  // scratch matrix
 843         getScaleRotate( tmp_scale, tmp_rot );
 844
 845         m1.m00 = tmp_rot[0];
 846         m1.m01 = tmp_rot[1];
 847         m1.m02 = tmp_rot[2];
 848
 849         m1.m10 = tmp_rot[3];
 850         m1.m11 = tmp_rot[4];
 851         m1.m12 = tmp_rot[5];
 852
 853         m1.m20 = tmp_rot[6];
 854         m1.m21 = tmp_rot[7];
 855         m1.m22 = tmp_rot[8];
 856
 857         t1.x = m03;
 858         t1.y = m13;
 859         t1.z = m23;
 860
 861         return( Matrix3d.max3( tmp_scale ));
 862
 863     }
 864
 865    /**
 866      * Performs an SVD normalization of this matrix to calculate
 867      * the rotation as a 3x3 matrix, the translation, and the scale.
 868      * None of the matrix values are modified.
 869      * @param m1  the normalized matrix representing the rotation
 870      * @param t1  the translation component
 871      * @return  the scale component of this transform
 872      */
 873     public final double get(Matrix3f m1, Vector3d t1){
 874
 875         double[]    tmp_rot = new double[9];  // scratch matrix
 876         double[]    tmp_scale = new double[3];  // scratch matrix
 877         getScaleRotate( tmp_scale, tmp_rot );
 878
 879         m1.m00 = (float)tmp_rot[0];
 880         m1.m01 = (float)tmp_rot[1];
 881         m1.m02 = (float)tmp_rot[2];
 882
 883         m1.m10 = (float)tmp_rot[3];
 884         m1.m11 = (float)tmp_rot[4];
 885         m1.m12 = (float)tmp_rot[5];
 886
 887         m1.m20 = (float)tmp_rot[6];
 888         m1.m21 = (float)tmp_rot[7];
 889         m1.m22 = (float)tmp_rot[8];
 890
 891         t1.x = m03;
 892         t1.y = m13;
 893         t1.z = m23;
 894
 895         return( Matrix3d.max3( tmp_scale ));
 896
 897     }
 898
 899     /**
 900      * Performs an SVD normalization of this matrix in order to acquire
 901      * the normalized rotational component; the values are placed into
 902      * the Quat4f parameter.
 903      * @param q1  quaternion into which the rotation component is placed
 904      */
 905     public final void get(Quat4f q1)
 906     {
 907
 908         double[]    tmp_rot = new double[9];  // scratch matrix
 909         double[]    tmp_scale = new double[3];  // scratch matrix
 910         getScaleRotate( tmp_scale, tmp_rot );
 911
 912         double ww;
 913
 914         ww = 0.25*(1.0 + tmp_rot[0] + tmp_rot[4] + tmp_rot[8]);
 915         if(!((ww<0?-ww:ww) < 1.0e-30)) {
 916           q1.w = (float)Math.sqrt(ww);
 917           ww = 0.25/q1.w;
 918           q1.x = (float)((tmp_rot[7] - tmp_rot[5])*ww);
 919           q1.y = (float)((tmp_rot[2] - tmp_rot[6])*ww);
 920           q1.z = (float)((tmp_rot[3] - tmp_rot[1])*ww);
 921           return;
 922         }
 923
 924         q1.w = 0.0f;
 925         ww = -0.5*(tmp_rot[4] + tmp_rot[8]);
 926         if(!((ww<0?-ww:ww) < 1.0e-30)) {
 927           q1.x =  (float)Math.sqrt(ww);
 928           ww = 0.5/q1.x;
 929           q1.y = (float)(tmp_rot[3]*ww);
 930           q1.z = (float)(tmp_rot[6]*ww);
 931           return;
 932         }
 933
 934         q1.x = 0.0f;
 935         ww = 0.5*(1.0 - tmp_rot[8]);
 936         if(!((ww<0?-ww:ww) < 1.0e-30)) {
 937           q1.y =  (float)(Math.sqrt(ww));
 938           q1.z = (float)(tmp_rot[7]/(2.0*q1.y));
 939           return;
 940         }
 941
 942         q1.y = 0.0f;
 943         q1.z = 1.0f;
 944
 945     }
 946
 947     /**
 948      * Performs an SVD normalization of q1 matrix in order to acquire
 949      * the normalized rotational component; the values are placed into
 950      * the Quat4d parameter.
 951      * @param q1  the quaternion into which the rotation component is placed
 952      */
 953     public final void get(Quat4d q1)
 954     {
 955         double[]    tmp_rot = new double[9];  // scratch matrix
 956         double[]    tmp_scale = new double[3];  // scratch matrix
 957
 958         getScaleRotate( tmp_scale, tmp_rot );
 959
 960         double ww;
 961
 962         ww = 0.25*(1.0 + tmp_rot[0] + tmp_rot[4] + tmp_rot[8]);
 963         if(!((ww<0?-ww:ww) < 1.0e-30)) {
 964           q1.w = Math.sqrt(ww);
 965           ww = 0.25/q1.w;
 966           q1.x = (tmp_rot[7] - tmp_rot[5])*ww;
 967           q1.y = (tmp_rot[2] - tmp_rot[6])*ww;
 968           q1.z = (tmp_rot[3] - tmp_rot[1])*ww;
 969           return;
 970         }
 971
 972         q1.w = 0.0f;
 973         ww = -0.5*(tmp_rot[4] + tmp_rot[8]);
 974         if(!((ww<0?-ww:ww) < 1.0e-30)) {
 975           q1.x =  Math.sqrt(ww);
 976           ww = 0.5/q1.x;
 977           q1.y = tmp_rot[3]*ww;
 978           q1.z = tmp_rot[6]*ww;
 979           return;
 980         }
 981
 982         q1.x = 0.0;
 983         ww = 0.5*(1.0 - tmp_rot[8]);
 984         if(!((ww<0?-ww:ww) < 1.0e-30)) {
 985           q1.y =  Math.sqrt(ww);
 986           q1.z = tmp_rot[7]/(2.0*q1.y);
 987           return;
 988         }
 989
 990         q1.y = 0.0;
 991         q1.z = 1.0;
 992     }
 993
 994    /**
 995      * Retrieves the translational components of this matrix.
 996      * @param trans  the vector that will receive the translational component
 997      */
 998     public final void get(Vector3d trans)
 999     {
1000         trans.x = m03;
1001         trans.y = m13;
1002         trans.z = m23;
1003     }
1004
1005   /**
1006     * Gets the upper 3x3 values of this matrix and places them into
1007     * the matrix m1.
1008     * @param m1  the matrix that will hold the values
1009     */
1010     public final void getRotationScale(Matrix3f m1)
1011     {
1012       m1.m00 = (float)m00; m1.m01 = (float)m01; m1.m02 = (float)m02;
1013       m1.m10 = (float)m10; m1.m11 = (float)m11; m1.m12 = (float)m12;
1014       m1.m20 = (float)m20; m1.m21 = (float)m21; m1.m22 = (float)m22;
1015     }
1016
1017   /**
1018     * Gets the upper 3x3 values of this matrix and places them into
1019     * the matrix m1.
1020     * @param m1  the matrix that will hold the values
1021     */
1022     public final void getRotationScale(Matrix3d m1)
1023     {
1024       m1.m00 = m00; m1.m01 = m01; m1.m02 = m02;
1025       m1.m10 = m10; m1.m11 = m11; m1.m12 = m12;
1026       m1.m20 = m20; m1.m21 = m21; m1.m22 = m22;
1027     }
1028
1029    /**
1030      * Performs an SVD normalization of this matrix to calculate
1031      * and return the uniform scale factor. If the matrix has non-uniform
1032      * scale factors, the largest of the x, y, and z scale factors will
1033      * be returned. This matrix is not modified.
1034      * @return  the scale factor of this matrix
1035      */
1036     public final double getScale()
1037     {
1038
1039         double[]    tmp_rot = new double[9];  // scratch matrix
1040         double[]    tmp_scale = new double[3];  // scratch matrix
1041         getScaleRotate( tmp_scale, tmp_rot );
1042
1043         return( Matrix3d.max3( tmp_scale ));
1044
1045     }
1046
1047   /**
1048    * Replaces the upper 3x3 matrix values of this matrix with the
1049    * values in the matrix m1.
1050    * @param m1  the matrix that will be the new upper 3x3
1051    */
1052     public final void setRotationScale(Matrix3d m1)
1053     {
1054       m00 = m1.m00; m01 = m1.m01; m02 = m1.m02;
1055       m10 = m1.m10; m11 = m1.m11; m12 = m1.m12;
1056       m20 = m1.m20; m21 = m1.m21; m22 = m1.m22;
1057     }
1058
1059   /**
1060    * Replaces the upper 3x3 matrix values of this matrix with the
1061    * values in the matrix m1.
1062    * @param m1  the matrix that will be the new upper 3x3
1063    */
1064     public final void setRotationScale(Matrix3f m1)
1065     {
1066       m00 = m1.m00; m01 = m1.m01; m02 = m1.m02;
1067       m10 = m1.m10; m11 = m1.m11; m12 = m1.m12;
1068       m20 = m1.m20; m21 = m1.m21; m22 = m1.m22;
1069     }
1070
1071    /**
1072      * Sets the scale component of the current matrix by factoring
1073      * out the current scale (by doing an SVD) from the rotational
1074      * component and multiplying by the new scale.
1075      * @param scale  the new scale amount
1076      */
1077     public final void setScale(double scale)
1078     {
1079         double[]    tmp_rot = new double[9];  // scratch matrix
1080         double[]    tmp_scale = new double[3];  // scratch matrix
1081
1082         getScaleRotate( tmp_scale, tmp_rot );
1083
1084         m00 = tmp_rot[0]*scale;
1085         m01 = tmp_rot[1]*scale;
1086         m02 = tmp_rot[2]*scale;
1087
1088         m10 = tmp_rot[3]*scale;
1089         m11 = tmp_rot[4]*scale;
1090         m12 = tmp_rot[5]*scale;
1091
1092         m20 = tmp_rot[6]*scale;
1093         m21 = tmp_rot[7]*scale;
1094         m22 = tmp_rot[8]*scale;
1095
1096     }
1097
1098     /**
1099      * Sets the specified row of this matrix4d to the four values provided.
1100      * @param row the row number to be modified (zero indexed)
1101      * @param x the first column element
1102      * @param y the second column element
1103      * @param z the third column element
1104      * @param w the fourth column element
1105      */
1106     public final void setRow(int row, double x, double y, double z, double w)
1107     {
1108         switch (row) {
1109         case 0:
1110             this.m00 = x;
1111             this.m01 = y;
1112             this.m02 = z;
1113             this.m03 = w;
1114             break;
1115
1116         case 1:
1117             this.m10 = x;
1118             this.m11 = y;
1119             this.m12 = z;
1120             this.m13 = w;
1121             break;
1122
1123         case 2:
1124             this.m20 = x;
1125             this.m21 = y;
1126             this.m22 = z;
1127             this.m23 = w;
1128             break;
1129
1130         case 3:
1131             this.m30 = x;
1132             this.m31 = y;
1133             this.m32 = z;
1134             this.m33 = w;
1135             break;
1136
1137         default:
1138             throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d4"));
1139
1140         }
1141     }
1142
1143     /**
1144      * Sets the specified row of this matrix4d to the Vector provided.
1145      * @param row the row number to be modified (zero indexed)
1146      * @param v the replacement row
1147      */
1148     public final void setRow(int row, Vector4d v)
1149     {
1150         switch (row) {
1151         case 0:
1152             this.m00 = v.x;
1153             this.m01 = v.y;
1154             this.m02 = v.z;
1155             this.m03 = v.w;
1156             break;
1157
1158         case 1:
1159             this.m10 = v.x;
1160             this.m11 = v.y;
1161             this.m12 = v.z;
1162             this.m13 = v.w;
1163             break;
1164
1165         case 2:
1166             this.m20 = v.x;
1167             this.m21 = v.y;
1168             this.m22 = v.z;
1169             this.m23 = v.w;
1170             break;
1171
1172         case 3:
1173             this.m30 = v.x;
1174             this.m31 = v.y;
1175             this.m32 = v.z;
1176             this.m33 = v.w;
1177             break;
1178
1179         default:
1180             throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d4"));
1181         }
1182     }
1183
1184     /**
1185      * Sets the specified row of this matrix4d to the four values provided.
1186      * @param row the row number to be modified (zero indexed)
1187      * @param v the replacement row
1188      */
1189     public final void setRow(int row, double v[])
1190     {
1191         switch (row) {
1192         case 0:
1193             this.m00 = v[0];
1194             this.m01 = v[1];
1195             this.m02 = v[2];
1196             this.m03 = v[3];
1197             break;
1198
1199         case 1:
1200             this.m10 = v[0];
1201             this.m11 = v[1];
1202             this.m12 = v[2];
1203             this.m13 = v[3];
1204             break;
1205
1206         case 2:
1207             this.m20 = v[0];
1208             this.m21 = v[1];
1209             this.m22 = v[2];
1210             this.m23 = v[3];
1211             break;
1212
1213         case 3:
1214             this.m30 = v[0];
1215             this.m31 = v[1];
1216             this.m32 = v[2];
1217             this.m33 = v[3];
1218             break;
1219
1220         default:
1221             throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d4"));
1222         }
1223     }
1224
1225     /**
1226      * Sets the specified column of this matrix4d to the four values provided.
1227      * @param column the column number to be modified (zero indexed)
1228      * @param x the first row element
1229      * @param y the second row element
1230      * @param z the third row element
1231      * @param w the fourth row element
1232      */
1233     public final void setColumn(int column, double x, double y, double z, double w)
1234     {
1235         switch (column) {
1236         case 0:
1237             this.m00 = x;
1238             this.m10 = y;
1239             this.m20 = z;
1240             this.m30 = w;
1241             break;
1242
1243         case 1:
1244             this.m01 = x;
1245             this.m11 = y;
1246             this.m21 = z;
1247             this.m31 = w;
1248             break;
1249
1250         case 2:
1251             this.m02 = x;
1252             this.m12 = y;
1253             this.m22 = z;
1254             this.m32 = w;
1255             break;
1256
1257         case 3:
1258             this.m03 = x;
1259             this.m13 = y;
1260             this.m23 = z;
1261             this.m33 = w;
1262             break;
1263
1264         default:
1265             throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d7"));
1266         }
1267     }
1268
1269     /**
1270      * Sets the specified column of this matrix4d to the vector provided.
1271      * @param column the column number to be modified (zero indexed)
1272      * @param v the replacement column
1273      */
1274     public final void setColumn(int column, Vector4d v)
1275     {
1276         switch (column) {
1277         case 0:
1278             this.m00 = v.x;
1279             this.m10 = v.y;
1280             this.m20 = v.z;
1281             this.m30 = v.w;
1282             break;
1283
1284         case 1:
1285             this.m01 = v.x;
1286             this.m11 = v.y;
1287             this.m21 = v.z;
1288             this.m31 = v.w;
1289             break;
1290
1291         case 2:
1292             this.m02 = v.x;
1293             this.m12 = v.y;
1294             this.m22 = v.z;
1295             this.m32 = v.w;
1296             break;
1297
1298         case 3:
1299             this.m03 = v.x;
1300             this.m13 = v.y;
1301             this.m23 = v.z;
1302             this.m33 = v.w;
1303             break;
1304
1305         default:
1306             throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d7"));
1307         }
1308     }
1309
1310     /**
1311      * Sets the specified column of this matrix4d to the four values provided.
1312      * @param column the column number to be modified (zero indexed)
1313      * @param v the replacement column
1314      */
1315     public final void setColumn(int column, double v[])
1316     {
1317         switch (column) {
1318         case 0:
1319             this.m00 = v[0];
1320             this.m10 = v[1];
1321             this.m20 = v[2];
1322             this.m30 = v[3];
1323             break;
1324
1325         case 1:
1326             this.m01 = v[0];
1327             this.m11 = v[1];
1328             this.m21 = v[2];
1329             this.m31 = v[3];
1330             break;
1331
1332         case 2:
1333             this.m02 = v[0];
1334             this.m12 = v[1];
1335             this.m22 = v[2];
1336             this.m32 = v[3];
1337             break;
1338
1339         case 3:
1340             this.m03 = v[0];
1341             this.m13 = v[1];
1342             this.m23 = v[2];
1343             this.m33 = v[3];
1344             break;
1345
1346         default:
1347             throw new ArrayIndexOutOfBoundsException(VecMathI18N.getString("Matrix4d7"));
1348         }
1349     }
1350
1351    /**
1352      *  Adds a scalar to each component of this matrix.
1353      *  @param scalar  the scalar adder
1354      */
1355     public final void add(double scalar)
1356     {
1357         m00 += scalar;
1358         m01 += scalar;
1359         m02 += scalar;
1360         m03 += scalar;
1361         m10 += scalar;
1362         m11 += scalar;
1363         m12 += scalar;
1364         m13 += scalar;
1365         m20 += scalar;
1366         m21 += scalar;
1367         m22 += scalar;
1368         m23 += scalar;
1369         m30 += scalar;
1370         m31 += scalar;
1371         m32 += scalar;
1372         m33 += scalar;
1373     }
1374
1375    /**
1376      *  Adds a scalar to each component of the matrix m1 and places
1377      *  the result into this.  Matrix m1 is not modified.
1378      *  @param scalar  the scalar adder
1379      *  @param m1  the original matrix values
1380      */
1381     public final void add(double scalar, Matrix4d m1)
1382     {
1383         this.m00 = m1.m00 +  scalar;
1384         this.m01 = m1.m01 +  scalar;
1385         this.m02 = m1.m02 +  scalar;
1386         this.m03 = m1.m03 +  scalar;
1387         this.m10 = m1.m10 +  scalar;
1388         this.m11 = m1.m11 +  scalar;
1389         this.m12 = m1.m12 +  scalar;
1390         this.m13 = m1.m13 +  scalar;
1391         this.m20 = m1.m20 +  scalar;
1392         this.m21 = m1.m21 +  scalar;
1393         this.m22 = m1.m22 +  scalar;
1394         this.m23 = m1.m23 +  scalar;
1395         this.m30 = m1.m30 +  scalar;
1396         this.m31 = m1.m31 +  scalar;
1397         this.m32 = m1.m32 +  scalar;
1398         this.m33 = m1.m33 +  scalar;
1399     }
1400
1401     /**
1402      * Sets the value of this matrix to the matrix sum of matrices m1 and m2.
1403      * @param m1 the first matrix
1404      * @param m2 the second matrix
1405      */
1406     public final void add(Matrix4d m1, Matrix4d m2)
1407     {
1408         this.m00 = m1.m00 + m2.m00;
1409         this.m01 = m1.m01 + m2.m01;
1410         this.m02 = m1.m02 + m2.m02;
1411         this.m03 = m1.m03 + m2.m03;
1412
1413         this.m10 = m1.m10 + m2.m10;
1414         this.m11 = m1.m11 + m2.m11;
1415         this.m12 = m1.m12 + m2.m12;
1416         this.m13 = m1.m13 + m2.m13;
1417
1418         this.m20 = m1.m20 + m2.m20;
1419         this.m21 = m1.m21 + m2.m21;
1420         this.m22 = m1.m22 + m2.m22;
1421         this.m23 = m1.m23 + m2.m23;
1422
1423         this.m30 = m1.m30 + m2.m30;
1424         this.m31 = m1.m31 + m2.m31;
1425         this.m32 = m1.m32 + m2.m32;
1426         this.m33 = m1.m33 + m2.m33;
1427     }
1428
1429     /**
1430      * Sets the value of this matrix to sum of itself and matrix m1.
1431      * @param m1 the other matrix
1432      */
1433     public final void add(Matrix4d m1)
1434     {
1435         this.m00 += m1.m00;
1436         this.m01 += m1.m01;
1437         this.m02 += m1.m02;
1438         this.m03 += m1.m03;
1439
1440         this.m10 += m1.m10;
1441         this.m11 += m1.m11;
1442         this.m12 += m1.m12;
1443         this.m13 += m1.m13;
1444
1445         this.m20 += m1.m20;
1446         this.m21 += m1.m21;
1447         this.m22 += m1.m22;
1448         this.m23 += m1.m23;
1449
1450         this.m30 += m1.m30;
1451         this.m31 += m1.m31;
1452         this.m32 += m1.m32;
1453         this.m33 += m1.m33;
1454     }
1455
1456     /**
1457      * Sets the value of this matrix to the matrix difference
1458      * of matrices m1 and m2.
1459      * @param m1 the first matrix
1460      * @param m2 the second matrix
1461      */
1462     public final void sub(Matrix4d m1, Matrix4d m2)
1463     {
1464         this.m00 = m1.m00 - m2.m00;
1465         this.m01 = m1.m01 - m2.m01;
1466         this.m02 = m1.m02 - m2.m02;
1467         this.m03 = m1.m03 - m2.m03;
1468
1469         this.m10 = m1.m10 - m2.m10;
1470         this.m11 = m1.m11 - m2.m11;
1471         this.m12 = m1.m12 - m2.m12;
1472         this.m13 = m1.m13 - m2.m13;
1473
1474         this.m20 = m1.m20 - m2.m20;
1475         this.m21 = m1.m21 - m2.m21;
1476         this.m22 = m1.m22 - m2.m22;
1477         this.m23 = m1.m23 - m2.m23;
1478
1479         this.m30 = m1.m30 - m2.m30;
1480         this.m31 = m1.m31 - m2.m31;
1481         this.m32 = m1.m32 - m2.m32;
1482         this.m33 = m1.m33 - m2.m33;
1483     }
1484
1485
1486     /**
1487      * Sets the value of this matrix to the matrix difference of itself
1488      * and matrix m1 (this = this - m1).
1489      * @param m1 the other matrix
1490      */
1491     public final void sub(Matrix4d m1)
1492     {
1493         this.m00 -= m1.m00;
1494         this.m01 -= m1.m01;
1495         this.m02 -= m1.m02;
1496         this.m03 -= m1.m03;
1497
1498         this.m10 -= m1.m10;
1499         this.m11 -= m1.m11;
1500         this.m12 -= m1.m12;
1501         this.m13 -= m1.m13;
1502
1503         this.m20 -= m1.m20;
1504         this.m21 -= m1.m21;
1505         this.m22 -= m1.m22;
1506         this.m23 -= m1.m23;
1507
1508         this.m30 -= m1.m30;
1509         this.m31 -= m1.m31;
1510         this.m32 -= m1.m32;
1511         this.m33 -= m1.m33;
1512     }
1513
1514     /**
1515      * Sets the value of this matrix to its transpose.
1516      */
1517     public final void transpose()
1518     {
1519         double temp;
1520
1521         temp = this.m10;
1522         this.m10 = this.m01;
1523         this.m01 = temp;
1524
1525         temp = this.m20;
1526         this.m20 = this.m02;
1527         this.m02 = temp;
1528
1529         temp = this.m30;
1530         this.m30 = this.m03;
1531         this.m03 = temp;
1532
1533         temp = this.m21;
1534         this.m21 = this.m12;
1535         this.m12 = temp;
1536
1537         temp = this.m31;
1538         this.m31 = this.m13;
1539         this.m13 = temp;
1540
1541         temp = this.m32;
1542         this.m32 = this.m23;
1543         this.m23 = temp;
1544     }
1545
1546     /**
1547      * Sets the value of this matrix to the transpose of the argument matrix
1548      * @param m1 the matrix to be transposed
1549      */
1550     public final void transpose(Matrix4d m1)
1551     {
1552         if (this != m1) {
1553             this.m00 = m1.m00;
1554             this.m01 = m1.m10;
1555             this.m02 = m1.m20;
1556             this.m03 = m1.m30;
1557
1558             this.m10 = m1.m01;
1559             this.m11 = m1.m11;
1560             this.m12 = m1.m21;
1561             this.m13 = m1.m31;
1562
1563             this.m20 = m1.m02;
1564             this.m21 = m1.m12;
1565             this.m22 = m1.m22;
1566             this.m23 = m1.m32;
1567
1568             this.m30 = m1.m03;
1569             this.m31 = m1.m13;
1570             this.m32 = m1.m23;
1571             this.m33 = m1.m33;
1572         } else
1573             this.transpose();
1574     }
1575
1576     /**
1577      *  Sets the values in this Matrix4d equal to the row-major
1578      *  array parameter (ie, the first four elements of the
1579      *  array will be copied into the first row of this matrix, etc.).
1580      *  @param m  the double precision array of length 16
1581      */
1582     public final void set(double[] m)
1583     {
1584           m00 = m[0];
1585           m01 = m[1];
1586           m02 = m[2];
1587           m03 = m[3];
1588           m10 = m[4];
1589           m11 = m[5];
1590           m12 = m[6];
1591           m13 = m[7];
1592           m20 = m[8];
1593           m21 = m[9];
1594           m22 = m[10];
1595           m23 = m[11];
1596           m30 = m[12];
1597           m31 = m[13];
1598           m32 = m[14];
1599           m33 = m[15];
1600     }
1601
1602     /**
1603      * Sets the rotational component (upper 3x3) of this matrix to the
1604      * matrix values in the single precision Matrix3f argument; the other
1605      * elements of this matrix are initialized as if this were an identity
1606      * matrix (i.e., affine matrix with no translational component).
1607      * @param m1   the double precision 3x3 matrix
1608      */
1609     public final void set(Matrix3f m1)
1610     {
1611        m00 = m1.m00; m01 = m1.m01; m02 = m1.m02; m03 = 0.0;
1612        m10 = m1.m10; m11 = m1.m11; m12 = m1.m12; m13 = 0.0;
1613        m20 = m1.m20; m21 = m1.m21; m22 = m1.m22; m23 = 0.0;
1614        m30 = 0.0;    m31 = 0.0   ; m32 = 0.0   ; m33 = 1.0;
1615     }
1616
1617     /**
1618      * Sets the rotational component (upper 3x3) of this matrix to the
1619      * matrix values in the double precision Matrix3d argument; the other
1620      * elements of this matrix are initialized as if this were an identity
1621      * matrix (i.e., affine matrix with no translational component).
1622      * @param m1   the double precision 3x3 matrix
1623      */
1624     public final void set(Matrix3d m1)
1625     {
1626        m00 = m1.m00; m01 = m1.m01; m02 = m1.m02; m03 = 0.0;
1627        m10 = m1.m10; m11 = m1.m11; m12 = m1.m12; m13 = 0.0;
1628        m20 = m1.m20; m21 = m1.m21; m22 = m1.m22; m23 = 0.0;
1629        m30 = 0.0;    m31 = 0.0   ; m32 = 0.0   ; m33 = 1.0;
1630     }
1631
1632     /**
1633      * Sets the value of this matrix to the matrix conversion of the
1634      * (double precision) quaternion argument.
1635      * @param q1 the quaternion to be converted
1636      */
1637     public final void set(Quat4d q1)
1638     {
1639         this.m00 = (1.0 - 2.0*q1.y*q1.y - 2.0*q1.z*q1.z);
1640         this.m10 = (2.0*(q1.x*q1.y + q1.w*q1.z));
1641         this.m20 = (2.0*(q1.x*q1.z - q1.w*q1.y));
1642
1643         this.m01 = (2.0*(q1.x*q1.y - q1.w*q1.z));
1644         this.m11 = (1.0 - 2.0*q1.x*q1.x - 2.0*q1.z*q1.z);
1645         this.m21 = (2.0*(q1.y*q1.z + q1.w*q1.x));
1646
1647         this.m02 = (2.0*(q1.x*q1.z + q1.w*q1.y));
1648         this.m12 = (2.0*(q1.y*q1.z - q1.w*q1.x));
1649         this.m22 = (1.0 - 2.0*q1.x*q1.x - 2.0*q1.y*q1.y);
1650
1651         this.m03 = 0.0;
1652         this.m13 = 0.0;
1653         this.m23 = 0.0;
1654
1655         this.m30 = 0.0;
1656         this.m31 = 0.0;
1657         this.m32 = 0.0;
1658         this.m33 = 1.0;
1659     }
1660
1661     /**
1662      * Sets the value of this matrix to the matrix conversion of the
1663      * double precision axis and angle argument.
1664      * @param a1 the axis and angle to be converted
1665      */
1666     public final void set(AxisAngle4d a1)
1667     {
1668       double mag = Math.sqrt( a1.x*a1.x + a1.y*a1.y + a1.z*a1.z);
1669
1670       if( mag < EPS ) {
1671          m00 = 1.0;
1672          m01 = 0.0;
1673          m02 = 0.0;
1674
1675          m10 = 0.0;
1676          m11 = 1.0;
1677          m12 = 0.0;
1678
1679          m20 = 0.0;
1680          m21 = 0.0;
1681          m22 = 1.0;
1682       } else {
1683          mag = 1.0/mag;
1684          double ax = a1.x*mag;
1685          double ay = a1.y*mag;
1686          double az = a1.z*mag;
1687
1688          double sinTheta = Math.sin(a1.angle);
1689          double cosTheta = Math.cos(a1.angle);
1690          double t = 1.0 - cosTheta;
1691
1692          double xz = ax * az;
1693          double xy = ax * ay;
1694          double yz = ay * az;
1695
1696          m00 = t * ax * ax + cosTheta;
1697          m01 = t * xy - sinTheta * az;
1698          m02 = t * xz + sinTheta * ay;
1699
1700          m10 = t * xy + sinTheta * az;
1701          m11 = t * ay * ay + cosTheta;
1702          m12 = t * yz - sinTheta * ax;
1703
1704          m20 = t * xz - sinTheta * ay;
1705          m21 = t * yz + sinTheta * ax;
1706          m22 = t * az * az + cosTheta;
1707       }
1708
1709       m03 = 0.0;
1710       m13 = 0.0;
1711       m23 = 0.0;
1712
1713       m30 = 0.0;
1714       m31 = 0.0;
1715       m32 = 0.0;
1716       m33 = 1.0;
1717     }
1718
1719     /**
1720      * Sets the value of this matrix to the matrix conversion of the
1721      * single precision quaternion argument.
1722      * @param q1 the quaternion to be converted
1723      */
1724     public final void set(Quat4f q1)
1725     {
1726         this.m00 = (1.0 - 2.0*q1.y*q1.y - 2.0*q1.z*q1.z);
1727         this.m10 = (2.0*(q1.x*q1.y + q1.w*q1.z));
1728         this.m20 = (2.0*(q1.x*q1.z - q1.w*q1.y));
1729
1730         this.m01 = (2.0*(q1.x*q1.y - q1.w*q1.z));
1731         this.m11 = (1.0 - 2.0*q1.x*q1.x - 2.0*q1.z*q1.z);
1732         this.m21 = (2.0*(q1.y*q1.z + q1.w*q1.x));
1733
1734         this.m02 = (2.0*(q1.x*q1.z + q1.w*q1.y));
1735         this.m12 = (2.0*(q1.y*q1.z - q1.w*q1.x));
1736         this.m22 = (1.0 - 2.0*q1.x*q1.x - 2.0*q1.y*q1.y);
1737
1738         this.m03 = 0.0;
1739         this.m13 = 0.0;
1740         this.m23 = 0.0;
1741
1742         this.m30 = 0.0;
1743         this.m31 = 0.0;
1744         this.m32 = 0.0;
1745         this.m33 = 1.0;
1746     }
1747
1748     /**
1749      * Sets the value of this matrix to the matrix conversion of the
1750      * single precision axis and angle argument.
1751      * @param a1 the axis and angle to be converted
1752      */
1753     public final void set(AxisAngle4f a1)
1754     {
1755       double mag = Math.sqrt( a1.x*a1.x + a1.y*a1.y + a1.z*a1.z);
1756
1757       if( mag < EPS ) {
1758          m00 = 1.0;
1759          m01 = 0.0;
1760          m02 = 0.0;
1761
1762          m10 = 0.0;
1763          m11 = 1.0;
1764          m12 = 0.0;
1765
1766          m20 = 0.0;
1767          m21 = 0.0;
1768          m22 = 1.0;
1769       } else {
1770          mag = 1.0/mag;
1771          double ax = a1.x*mag;
1772          double ay = a1.y*mag;
1773          double az = a1.z*mag;
1774
1775          double sinTheta = Math.sin((double)a1.angle);
1776          double cosTheta = Math.cos((double)a1.angle);
1777          double t = 1.0 - cosTheta;
1778
1779          double xz = ax * az;
1780          double xy = ax * ay;
1781          double yz = ay * az;
1782
1783          m00 = t * ax * ax + cosTheta;
1784          m01 = t * xy - sinTheta * az;
1785          m02 = t * xz + sinTheta * ay;
1786
1787          m10 = t * xy + sinTheta * az;
1788          m11 = t * ay * ay + cosTheta;
1789          m12 = t * yz - sinTheta * ax;
1790
1791          m20 = t * xz - sinTheta * ay;
1792          m21 = t * yz + sinTheta * ax;
1793          m22 = t * az * az + cosTheta;
1794       }
1795       m03 = 0.0;
1796       m13 = 0.0;
1797       m23 = 0.0;
1798
1799       m30 = 0.0;
1800       m31 = 0.0;
1801       m32 = 0.0;
1802       m33 = 1.0;
1803     }
1804
1805     /**
1806      * Sets the value of this matrix from the rotation expressed
1807      * by the quaternion q1, the translation t1, and the scale s.
1808      * @param q1 the rotation expressed as a quaternion
1809      * @param t1 the translation
1810      * @param s the scale value
1811      */
1812     public final void set(Quat4d q1, Vector3d t1, double s)
1813     {
1814         this.m00 = s*(1.0 - 2.0*q1.y*q1.y - 2.0*q1.z*q1.z);
1815         this.m10 = s*(2.0*(q1.x*q1.y + q1.w*q1.z));
1816         this.m20 = s*(2.0*(q1.x*q1.z - q1.w*q1.y));
1817
1818         this.m01 = s*(2.0*(q1.x*q1.y - q1.w*q1.z));
1819         this.m11 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.z*q1.z);
1820         this.m21 = s*(2.0*(q1.y*q1.z + q1.w*q1.x));
1821
1822         this.m02 = s*(2.0*(q1.x*q1.z + q1.w*q1.y));
1823         this.m12 = s*(2.0*(q1.y*q1.z - q1.w*q1.x));
1824         this.m22 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.y*q1.y);
1825
1826         this.m03 = t1.x;
1827         this.m13 = t1.y;
1828         this.m23 = t1.z;
1829
1830         this.m30 = 0.0;
1831         this.m31 = 0.0;
1832         this.m32 = 0.0;
1833         this.m33 = 1.0;
1834     }
1835
1836     /**
1837      * Sets the value of this matrix from the rotation expressed
1838      * by the quaternion q1, the translation t1, and the scale s.
1839      * @param q1 the rotation expressed as a quaternion
1840      * @param t1 the translation
1841      * @param s the scale value
1842      */
1843     public final void set(Quat4f q1, Vector3d t1, double s)
1844     {
1845         this.m00 = s*(1.0 - 2.0*q1.y*q1.y - 2.0*q1.z*q1.z);
1846         this.m10 = s*(2.0*(q1.x*q1.y + q1.w*q1.z));
1847         this.m20 = s*(2.0*(q1.x*q1.z - q1.w*q1.y));
1848
1849         this.m01 = s*(2.0*(q1.x*q1.y - q1.w*q1.z));
1850         this.m11 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.z*q1.z);
1851         this.m21 = s*(2.0*(q1.y*q1.z + q1.w*q1.x));
1852
1853         this.m02 = s*(2.0*(q1.x*q1.z + q1.w*q1.y));
1854         this.m12 = s*(2.0*(q1.y*q1.z - q1.w*q1.x));
1855         this.m22 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.y*q1.y);
1856
1857         this.m03 = t1.x;
1858         this.m13 = t1.y;
1859         this.m23 = t1.z;
1860
1861         this.m30 = 0.0;
1862         this.m31 = 0.0;
1863         this.m32 = 0.0;
1864         this.m33 = 1.0;
1865     }
1866
1867     /**
1868      * Sets the value of this matrix from the rotation expressed
1869      * by the quaternion q1, the translation t1, and the scale s.
1870      * @param q1 the rotation expressed as a quaternion
1871      * @param t1 the translation
1872      * @param s the scale value
1873      */
1874     public final void set(Quat4f q1, Vector3f t1, float s)
1875    {
1876         this.m00 = s*(1.0 - 2.0*q1.y*q1.y - 2.0*q1.z*q1.z);
1877         this.m10 = s*(2.0*(q1.x*q1.y + q1.w*q1.z));
1878         this.m20 = s*(2.0*(q1.x*q1.z - q1.w*q1.y));
1879
1880         this.m01 = s*(2.0*(q1.x*q1.y - q1.w*q1.z));
1881         this.m11 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.z*q1.z);
1882         this.m21 = s*(2.0*(q1.y*q1.z + q1.w*q1.x));
1883
1884         this.m02 = s*(2.0*(q1.x*q1.z + q1.w*q1.y));
1885         this.m12 = s*(2.0*(q1.y*q1.z - q1.w*q1.x));
1886         this.m22 = s*(1.0 - 2.0*q1.x*q1.x - 2.0*q1.y*q1.y);
1887
1888         this.m03 = t1.x;
1889         this.m13 = t1.y;
1890         this.m23 = t1.z;
1891
1892         this.m30 = 0.0;
1893         this.m31 = 0.0;
1894         this.m32 = 0.0;
1895         this.m33 = 1.0;
1896    }
1897
1898     /**
1899      * Sets the value of this matrix to a copy of the
1900      * passed matrix m1.
1901      * @param m1 the matrix4f
1902      */
1903     public final void set(Matrix4f m1)
1904     {
1905         this.m00 = m1.m00;
1906         this.m01 = m1.m01;
1907         this.m02 = m1.m02;
1908         this.m03 = m1.m03;
1909
1910         this.m10 = m1.m10;
1911         this.m11 = m1.m11;
1912         this.m12 = m1.m12;
1913         this.m13 = m1.m13;
1914
1915         this.m20 = m1.m20;
1916         this.m21 = m1.m21;
1917         this.m22 = m1.m22;
1918         this.m23 = m1.m23;
1919
1920         this.m30 = m1.m30;
1921         this.m31 = m1.m31;
1922         this.m32 = m1.m32;
1923         this.m33 = m1.m33;
1924     }
1925
1926     /**
1927      * Sets the value of this matrix to a copy of the
1928      * passed matrix m1.
1929      * @param m1 the matrix to be copied
1930      */
1931     public final void set(Matrix4d m1)
1932     {
1933         this.m00 = m1.m00;
1934         this.m01 = m1.m01;
1935         this.m02 = m1.m02;
1936         this.m03 = m1.m03;
1937
1938         this.m10 = m1.m10;
1939         this.m11 = m1.m11;
1940         this.m12 = m1.m12;
1941         this.m13 = m1.m13;
1942
1943         this.m20 = m1.m20;
1944         this.m21 = m1.m21;
1945         this.m22 = m1.m22;
1946         this.m23 = m1.m23;
1947
1948         this.m30 = m1.m30;
1949         this.m31 = m1.m31;
1950         this.m32 = m1.m32;
1951         this.m33 = m1.m33;
1952     }
1953
1954   /**
1955    * Sets the value of this matrix to the matrix inverse
1956    * of the passed (user declared) matrix m1.
1957    * @param m1 the matrix to be inverted
1958    */
1959   public final void invert(Matrix4d m1)
1960   {
1961
1962      invertGeneral( m1);
1963   }
1964
1965   /**
1966    * Inverts this matrix in place.
1967    */
1968   public final void invert()
1969   {
1970      invertGeneral( this );
1971   }
1972
1973     /**
1974      * General invert routine.  Inverts m1 and places the result in "this".
1975      * Note that this routine handles both the "this" version and the
1976      * non-"this" version.
1977      *
1978      * Also note that since this routine is slow anyway, we won't worry
1979      * about allocating a little bit of garbage.
1980      */
1981     final void invertGeneral(Matrix4d  m1) {
1982         double result[] = new double[16];
1983         int row_perm[] = new int[4];
1984         int i, r, c;
1985
1986         // Use LU decomposition and backsubstitution code specifically
1987         // for floating-point 4x4 matrices.
1988         double[]    tmp = new double[16];  // scratch matrix
1989         // Copy source matrix to t1tmp
1990         tmp[0] = m1.m00;
1991         tmp[1] = m1.m01;
1992         tmp[2] = m1.m02;
1993         tmp[3] = m1.m03;
1994
1995         tmp[4] = m1.m10;
1996         tmp[5] = m1.m11;
1997         tmp[6] = m1.m12;
1998         tmp[7] = m1.m13;
1999
2000         tmp[8] = m1.m20;
2001         tmp[9] = m1.m21;
2002         tmp[10] = m1.m22;
2003         tmp[11] = m1.m23;
2004
2005         tmp[12] = m1.m30;
2006         tmp[13] = m1.m31;
2007         tmp[14] = m1.m32;
2008         tmp[15] = m1.m33;
2009
2010         // Calculate LU decomposition: Is the matrix singular?
2011         if (!luDecomposition(tmp, row_perm)) {
2012             // Matrix has no inverse
2013             throw new SingularMatrixException(VecMathI18N.getString("Matrix4d10"));
2014         }
2015
2016         // Perform back substitution on the identity matrix
2017         for(i=0;i<16;i++) result[i] = 0.0;
2018         result[0] = 1.0; result[5] = 1.0; result[10] = 1.0; result[15] = 1.0;
2019         luBacksubstitution(tmp, row_perm, result);
2020
2021         this.m00 = result[0];
2022         this.m01 = result[1];
2023         this.m02 = result[2];
2024         this.m03 = result[3];
2025
2026         this.m10 = result[4];
2027         this.m11 = result[5];
2028         this.m12 = result[6];
2029         this.m13 = result[7];
2030
2031         this.m20 = result[8];
2032         this.m21 = result[9];
2033         this.m22 = result[10];
2034         this.m23 = result[11];
2035
2036         this.m30 = result[12];
2037         this.m31 = result[13];
2038         this.m32 = result[14];
2039         this.m33 = result[15];
2040
2041     }
2042
2043     /**
2044      * Given a 4x4 array "matrix0", this function replaces it with the
2045      * LU decomposition of a row-wise permutation of itself.  The input
2046      * parameters are "matrix0" and "dimen".  The array "matrix0" is also
2047      * an output parameter.  The vector "row_perm[4]" is an output
2048      * parameter that contains the row permutations resulting from partial
2049      * pivoting.  The output parameter "even_row_xchg" is 1 when the
2050      * number of row exchanges is even, or -1 otherwise.  Assumes data
2051      * type is always double.
2052      *
2053      * This function is similar to luDecomposition, except that it
2054      * is tuned specifically for 4x4 matrices.
2055      *
2056      * @return true if the matrix is nonsingular, or false otherwise.
2057      */
2058     //
2059     // Reference: Press, Flannery, Teukolsky, Vetterling,
2060     //        _Numerical_Recipes_in_C_, Cambridge University Press,
2061     //        1988, pp 40-45.
2062     //
2063     static boolean luDecomposition(double[] matrix0,
2064                                    int[] row_perm) {
2065
2066         double row_scale[] = new double[4];
2067
2068         // Determine implicit scaling information by looping over rows
2069         {
2070             int i, j;
2071             int ptr, rs;
2072             double big, temp;
2073
2074             ptr = 0;
2075             rs = 0;
2076
2077             // For each row ...
2078             i = 4;
2079             while (i-- != 0) {
2080                 big = 0.0;
2081
2082                 // For each column, find the largest element in the row
2083                 j = 4;
2084                 while (j-- != 0) {
2085                     temp = matrix0[ptr++];
2086                     temp = Math.abs(temp);
2087                     if (temp > big) {
2088                         big = temp;
2089                     }
2090                 }
2091
2092                 // Is the matrix singular?
2093                 if (big == 0.0) {
2094                     return false;
2095                 }
2096                 row_scale[rs++] = 1.0 / big;
2097             }
2098         }
2099
2100         {
2101             int j;
2102             int mtx;
2103
2104             mtx = 0;
2105
2106             // For all columns, execute Crout's method
2107             for (j = 0; j < 4; j++) {
2108                 int i, imax, k;
2109                 int target, p1, p2;
2110                 double sum, big, temp;
2111
2112                 // Determine elements of upper diagonal matrix U
2113                 for (i = 0; i < j; i++) {
2114                     target = mtx + (4*i) + j;
2115                     sum = matrix0[target];
2116                     k = i;
2117                     p1 = mtx + (4*i);
2118                     p2 = mtx + j;
2119                     while (k-- != 0) {
2120                         sum -= matrix0[p1] * matrix0[p2];
2121                         p1++;
2122                         p2 += 4;
2123                     }
2124                     matrix0[target] = sum;
2125                 }
2126
2127                 // Search for largest pivot element and calculate
2128                 // intermediate elements of lower diagonal matrix L.
2129                 big = 0.0;
2130                 imax = -1;
2131                 for (i = j; i < 4; i++) {
2132                     target = mtx + (4*i) + j;
2133                     sum = matrix0[target];
2134                     k = j;
2135                     p1 = mtx + (4*i);
2136                     p2 = mtx + j;
2137                     while (k-- != 0) {
2138                         sum -= matrix0[p1] * matrix0[p2];
2139                         p1++;
2140                         p2 += 4;
2141                     }
2142                     matrix0[target] = sum;
2143
2144                     // Is this the best pivot so far?
2145                     if ((temp = row_scale[i] * Math.abs(sum)) >= big) {
2146                         big = temp;
2147                         imax = i;
2148                     }
2149                 }
2150
2151                 if (imax < 0) {
2152                     throw new RuntimeException(VecMathI18N.getString("Matrix4d11"));
2153                 }
2154
2155                 // Is a row exchange necessary?
2156                 if (j != imax) {
2157                     // Yes: exchange rows
2158                     k = 4;
2159                     p1 = mtx + (4*imax);
2160                     p2 = mtx + (4*j);
2161                     while (k-- != 0) {
2162                         temp = matrix0[p1];
2163                         matrix0[p1++] = matrix0[p2];
2164                         matrix0[p2++] = temp;
2165                     }
2166
2167                     // Record change in scale factor
2168                     row_scale[imax] = row_scale[j];
2169                 }
2170
2171                 // Record row permutation
2172                 row_perm[j] = imax;
2173
2174                 // Is the matrix singular
2175                 if (matrix0[(mtx + (4*j) + j)] == 0.0) {
2176                     return false;
2177                 }
2178
2179                 // Divide elements of lower diagonal matrix L by pivot
2180                 if (j != (4-1)) {
2181                     temp = 1.0 / (matrix0[(mtx + (4*j) + j)]);
2182                     target = mtx + (4*(j+1)) + j;
2183                     i = 3 - j;
2184                     while (i-- != 0) {
2185                         matrix0[target] *= temp;
2186                         target += 4;
2187                     }
2188                 }
2189             }
2190         }
2191
2192         return true;
2193     }
2194
2195     /**
2196      * Solves a set of linear equations.  The input parameters "matrix1",
2197      * and "row_perm" come from luDecompostionD4x4 and do not change
2198      * here.  The parameter "matrix2" is a set of column vectors assembled
2199      * into a 4x4 matrix of floating-point values.  The procedure takes each
2200      * column of "matrix2" in turn and treats it as the right-hand side of the
2201      * matrix equation Ax = LUx = b.  The solution vector replaces the
2202      * original column of the matrix.
2203      *
2204      * If "matrix2" is the identity matrix, the procedure replaces its contents
2205      * with the inverse of the matrix from which "matrix1" was originally
2206      * derived.
2207      */
2208     //
2209     // Reference: Press, Flannery, Teukolsky, Vetterling,
2210     //        _Numerical_Recipes_in_C_, Cambridge University Press,
2211     //        1988, pp 44-45.
2212     //
2213     static void luBacksubstitution(double[] matrix1,
2214                                    int[] row_perm,
2215                                    double[] matrix2) {
2216
2217         int i, ii, ip, j, k;
2218         int rp;
2219         int cv, rv;
2220
2221         //      rp = row_perm;
2222         rp = 0;
2223
2224         // For each column vector of matrix2 ...
2225         for (k = 0; k < 4; k++) {
2226             //      cv = &(matrix2[0][k]);
2227             cv = k;
2228             ii = -1;
2229
2230             // Forward substitution
2231             for (i = 0; i < 4; i++) {
2232                 double sum;
2233
2234                 ip = row_perm[rp+i];
2235                 sum = matrix2[cv+4*ip];
2236                 matrix2[cv+4*ip] = matrix2[cv+4*i];
2237                 if (ii >= 0) {
2238                     //              rv = &(matrix1[i][0]);
2239                     rv = i*4;
2240                     for (j = ii; j <= i-1; j++) {
2241                         sum -= matrix1[rv+j] * matrix2[cv+4*j];
2242                     }
2243                 }
2244                 else if (sum != 0.0) {
2245                     ii = i;
2246                 }
2247                 matrix2[cv+4*i] = sum;
2248             }
2249
2250             // Backsubstitution
2251             //      rv = &(matrix1[3][0]);
2252             rv = 3*4;
2253             matrix2[cv+4*3] /= matrix1[rv+3];
2254
2255             rv -= 4;
2256             matrix2[cv+4*2] = (matrix2[cv+4*2] -
2257                             matrix1[rv+3] * matrix2[cv+4*3]) / matrix1[rv+2];
2258
2259             rv -= 4;
2260             matrix2[cv+4*1] = (matrix2[cv+4*1] -
2261                             matrix1[rv+2] * matrix2[cv+4*2] -
2262                             matrix1[rv+3] * matrix2[cv+4*3]) / matrix1[rv+1];
2263
2264             rv -= 4;
2265             matrix2[cv+4*0] = (matrix2[cv+4*0] -
2266                             matrix1[rv+1] * matrix2[cv+4*1] -
2267                             matrix1[rv+2] * matrix2[cv+4*2] -
2268                             matrix1[rv+3] * matrix2[cv+4*3]) / matrix1[rv+0];
2269         }
2270     }
2271
2272     /**
2273      * Computes the determinant of this matrix.
2274      * @return the determinant of the matrix
2275      */
2276     public final double determinant()
2277     {
2278        double det;
2279
2280        // cofactor exapainsion along first row
2281
2282         det = m00*(m11*m22*m33+ m12*m23*m31 + m13*m21*m32
2283                  - m13*m22*m31 -m11*m23*m32 - m12*m21*m33);
2284         det -= m01*(m10*m22*m33+ m12*m23*m30 + m13*m20*m32
2285                   - m13*m22*m30 -m10*m23*m32 - m12*m20*m33);
2286         det += m02*(m10*m21*m33+ m11*m23*m30 + m13*m20*m31
2287                   - m13*m21*m30 -m10*m23*m31 - m11*m20*m33);
2288         det -= m03*(m10*m21*m32+ m11*m22*m30 + m12*m20*m31
2289                   - m12*m21*m30 -m10*m22*m31 - m11*m20*m32);
2290
2291         return( det );
2292     }
2293
2294     /**
2295      * Sets the value of this matrix to a scale matrix with the
2296      * passed scale amount.
2297      * @param scale the scale factor for the matrix
2298      */
2299     public final void set(double scale)
2300     {
2301         this.m00 = scale;
2302         this.m01 = 0.0;
2303         this.m02 = 0.0;
2304         this.m03 = 0.0;
2305
2306         this.m10 = 0.0;
2307         this.m11 = scale;
2308         this.m12 = 0.0;
2309         this.m13 = 0.0;
2310
2311         this.m20 = 0.0;
2312         this.m21 = 0.0;
2313         this.m22 = scale;
2314         this.m23 = 0.0;
2315
2316         this.m30 = 0.0;
2317         this.m31 = 0.0;
2318         this.m32 = 0.0;
2319         this.m33 = 1.0;
2320     }
2321
2322     /**
2323      * Sets the value of this matrix to a translate matrix by the
2324      * passed translation value.
2325      * @param v1 the translation amount
2326      */
2327     public final void set(Vector3d v1)
2328     {
2329         this.m00 = 1.0;
2330         this.m01 = 0.0;
2331         this.m02 = 0.0;
2332         this.m03 = v1.x;
2333
2334         this.m10 = 0.0;
2335         this.m11 = 1.0;
2336         this.m12 = 0.0;
2337         this.m13 = v1.y;
2338
2339         this.m20 = 0.0;
2340         this.m21 = 0.0;
2341         this.m22 = 1.0;
2342         this.m23 = v1.z;
2343
2344         this.m30 = 0.0;
2345         this.m31 = 0.0;
2346         this.m32 = 0.0;
2347         this.m33 = 1.0;
2348     }
2349
2350     /**
2351      * Sets the value of this transform to a scale and translation matrix;
2352      * the scale is not applied to the translation and all of the matrix
2353      * values are modified.
2354      * @param scale the scale factor for the matrix
2355      * @param v1 the translation amount
2356      */
2357     public final void set(double scale, Vector3d v1)
2358     {
2359         this.m00 = scale;
2360         this.m01 = 0.0;
2361         this.m02 = 0.0;
2362         this.m03 = v1.x;
2363
2364         this.m10 = 0.0;
2365         this.m11 = scale;
2366         this.m12 = 0.0;
2367         this.m13 = v1.y;
2368
2369         this.m20 = 0.0;
2370         this.m21 = 0.0;
2371         this.m22 = scale;
2372         this.m23 = v1.z;
2373
2374         this.m30 = 0.0;
2375         this.m31 = 0.0;
2376         this.m32 = 0.0;
2377         this.m33 = 1.0;
2378     }
2379
2380    /**
2381      * Sets the value of this transform to a scale and translation matrix;
2382      * the translation is scaled by the scale factor and all of the matrix
2383      * values are modified.
2384      * @param v1 the translation amount
2385      * @param scale the scale factor for the matrix
2386      */
2387     public final void set(Vector3d v1, double scale)
2388     {
2389         this.m00 = scale;
2390         this.m01 = 0.0;
2391         this.m02 = 0.0;
2392         this.m03 = scale*v1.x;
2393
2394         this.m10 = 0.0;
2395         this.m11 = scale;
2396         this.m12 = 0.0;
2397         this.m13 = scale*v1.y;
2398
2399         this.m20 = 0.0;
2400         this.m21 = 0.0;
2401         this.m22 = scale;
2402         this.m23 = scale*v1.z;
2403
2404         this.m30 = 0.0;
2405         this.m31 = 0.0;
2406         this.m32 = 0.0;
2407         this.m33 = 1.0;
2408     }
2409
2410    /**
2411      * Sets the value of this matrix from the rotation expressed by
2412      * the rotation matrix m1, the translation t1, and the scale factor.
2413      * The translation is not modified by the scale.
2414      * @param m1  the rotation component
2415      * @param t1  the translation component
2416      * @param scale  the scale component
2417      */
2418     public final void set(Matrix3f m1, Vector3f t1, float scale)
2419     {
2420         this.m00 = m1.m00*scale;
2421         this.m01 = m1.m01*scale;
2422         this.m02 = m1.m02*scale;
2423         this.m03 = t1.x;
2424
2425         this.m10 = m1.m10*scale;
2426         this.m11 = m1.m11*scale;
2427         this.m12 = m1.m12*scale;
2428         this.m13 = t1.y;
2429
2430         this.m20 = m1.m20*scale;
2431         this.m21 = m1.m21*scale;
2432         this.m22 = m1.m22*scale;
2433         this.m23 = t1.z;
2434
2435         this.m30 = 0.0;
2436         this.m31 = 0.0;
2437         this.m32 = 0.0;
2438         this.m33 = 1.0;
2439     }
2440
2441
2442    /**
2443      * Sets the value of this matrix from the rotation expressed by
2444      * the rotation matrix m1, the translation t1, and the scale factor.
2445      * The translation is not modified by the scale.
2446      * @param m1  the rotation component
2447      * @param t1  the translation component
2448      * @param scale  the scale component
2449      */
2450     public final void set(Matrix3d m1, Vector3d t1, double scale)
2451     {
2452         this.m00 = m1.m00*scale;
2453         this.m01 = m1.m01*scale;
2454         this.m02 = m1.m02*scale;
2455         this.m03 = t1.x;
2456
2457         this.m10 = m1.m10*scale;
2458         this.m11 = m1.m11*scale;
2459         this.m12 = m1.m12*scale;
2460         this.m13 = t1.y;
2461
2462         this.m20 = m1.m20*scale;
2463         this.m21 = m1.m21*scale;
2464         this.m22 = m1.m22*scale;
2465         this.m23 = t1.z;
2466
2467         this.m30 = 0.0;
2468         this.m31 = 0.0;
2469         this.m32 = 0.0;
2470         this.m33 = 1.0;
2471     }
2472
2473    /**
2474      * Modifies the translational components of this matrix to the values
2475      * of the Vector3d argument; the other values of this matrix are not
2476      * modified.
2477      * @param trans  the translational component
2478      */
2479     public final void setTranslation(Vector3d trans)
2480     {
2481        m03 = trans.x;
2482        m13 = trans.y;
2483        m23 = trans.z;
2484     }
2485
2486     /**
2487      * Sets the value of this matrix to a counter-clockwise rotation
2488      * about the x axis.
2489      * @param angle the angle to rotate about the X axis in radians
2490      */
2491     public final void rotX(double angle)
2492     {
2493         double  sinAngle, cosAngle;
2494
2495         sinAngle = Math.sin(angle);
2496         cosAngle = Math.cos(angle);
2497
2498         this.m00 = 1.0;
2499         this.m01 = 0.0;
2500         this.m02 = 0.0;
2501         this.m03 = 0.0;
2502
2503         this.m10 = 0.0;
2504         this.m11 = cosAngle;
2505         this.m12 = -sinAngle;
2506         this.m13 = 0.0;
2507
2508         this.m20 = 0.0;
2509         this.m21 = sinAngle;
2510         this.m22 = cosAngle;
2511         this.m23 = 0.0;
2512
2513         this.m30 = 0.0;
2514         this.m31 = 0.0;
2515         this.m32 = 0.0;
2516         this.m33 = 1.0;
2517     }
2518
2519     /**
2520      * Sets the value of this matrix to a counter-clockwise rotation
2521      * about the y axis.
2522      * @param angle the angle to rotate about the Y axis in radians
2523      */
2524     public final void rotY(double angle)
2525     {
2526         double  sinAngle, cosAngle;
2527
2528         sinAngle = Math.sin(angle);
2529         cosAngle = Math.cos(angle);
2530
2531         this.m00 = cosAngle;
2532         this.m01 = 0.0;
2533         this.m02 = sinAngle;
2534         this.m03 = 0.0;
2535
2536         this.m10 = 0.0;
2537         this.m11 = 1.0;
2538         this.m12 = 0.0;
2539         this.m13 = 0.0;
2540
2541         this.m20 = -sinAngle;
2542         this.m21 = 0.0;
2543         this.m22 = cosAngle;
2544         this.m23 = 0.0;
2545
2546         this.m30 = 0.0;
2547         this.m31 = 0.0;
2548         this.m32 = 0.0;
2549         this.m33 = 1.0;
2550     }
2551
2552     /**
2553      * Sets the value of this matrix to a counter-clockwise rotation
2554      * about the z axis.
2555      * @param angle the angle to rotate about the Z axis in radians
2556      */
2557     public final void rotZ(double angle)
2558     {
2559         double  sinAngle, cosAngle;
2560
2561         sinAngle = Math.sin(angle);
2562         cosAngle = Math.cos(angle);
2563
2564         this.m00 = cosAngle;
2565         this.m01 = -sinAngle;
2566         this.m02 = 0.0;
2567         this.m03 = 0.0;
2568
2569         this.m10 = sinAngle;
2570         this.m11 = cosAngle;
2571         this.m12 = 0.0;
2572         this.m13 = 0.0;
2573
2574         this.m20 = 0.0;
2575         this.m21 = 0.0;
2576         this.m22 = 1.0;
2577         this.m23 = 0.0;
2578
2579         this.m30 = 0.0;
2580         this.m31 = 0.0;
2581         this.m32 = 0.0;
2582         this.m33 = 1.0;
2583     }
2584
2585    /**
2586      * Multiplies each element of this matrix by a scalar.
2587      * @param scalar  the scalar multiplier.
2588      */
2589     public final void mul(double scalar)
2590     {
2591       m00 *= scalar;
2592       m01 *= scalar;
2593       m02 *= scalar;
2594       m03 *= scalar;
2595       m10 *= scalar;
2596       m11 *= scalar;
2597       m12 *= scalar;
2598       m13 *= scalar;
2599       m20 *= scalar;
2600       m21 *= scalar;
2601       m22 *= scalar;
2602       m23 *= scalar;
2603       m30 *= scalar;
2604       m31 *= scalar;
2605       m32 *= scalar;
2606       m33 *= scalar;
2607     }
2608
2609    /**
2610      * Multiplies each element of matrix m1 by a scalar and places
2611      * the result into this.  Matrix m1 is not modified.
2612      * @param scalar  the scalar multiplier
2613      * @param m1  the original matrix
2614      */
2615     public final void mul(double scalar, Matrix4d m1)
2616     {
2617       this.m00 = m1.m00 * scalar;
2618       this.m01 = m1.m01 * scalar;
2619       this.m02 = m1.m02 * scalar;
2620       this.m03 = m1.m03 * scalar;
2621       this.m10 = m1.m10 * scalar;
2622       this.m11 = m1.m11 * scalar;
2623       this.m12 = m1.m12 * scalar;
2624       this.m13 = m1.m13 * scalar;
2625       this.m20 = m1.m20 * scalar;
2626       this.m21 = m1.m21 * scalar;
2627       this.m22 = m1.m22 * scalar;
2628       this.m23 = m1.m23 * scalar;
2629       this.m30 = m1.m30 * scalar;
2630       this.m31 = m1.m31 * scalar;
2631       this.m32 = m1.m32 * scalar;
2632       this.m33 = m1.m33 * scalar;
2633     }
2634
2635     /**
2636      * Sets the value of this matrix to the result of multiplying itself
2637      * with matrix m1.
2638      * @param m1 the other matrix
2639      */
2640     public final void mul(Matrix4d m1)
2641     {
2642             double      m00, m01, m02, m03,
2643                         m10, m11, m12, m13,
2644                         m20, m21, m22, m23,
2645                         m30, m31, m32, m33;  // vars for temp result matrix
2646
2647         m00 = this.m00*m1.m00 + this.m01*m1.m10 +
2648               this.m02*m1.m20 + this.m03*m1.m30;
2649         m01 = this.m00*m1.m01 + this.m01*m1.m11 +
2650               this.m02*m1.m21 + this.m03*m1.m31;
2651         m02 = this.m00*m1.m02 + this.m01*m1.m12 +
2652               this.m02*m1.m22 + this.m03*m1.m32;
2653         m03 = this.m00*m1.m03 + this.m01*m1.m13 +
2654               this.m02*m1.m23 + this.m03*m1.m33;
2655
2656         m10 = this.m10*m1.m00 + this.m11*m1.m10 +
2657               this.m12*m1.m20 + this.m13*m1.m30;
2658         m11 = this.m10*m1.m01 + this.m11*m1.m11 +
2659               this.m12*m1.m21 + this.m13*m1.m31;
2660         m12 = this.m10*m1.m02 + this.m11*m1.m12 +
2661               this.m12*m1.m22 + this.m13*m1.m32;
2662         m13 = this.m10*m1.m03 + this.m11*m1.m13 +
2663               this.m12*m1.m23 + this.m13*m1.m33;
2664
2665         m20 = this.m20*m1.m00 + this.m21*m1.m10 +
2666               this.m22*m1.m20 + this.m23*m1.m30;
2667         m21 = this.m20*m1.m01 + this.m21*m1.m11 +
2668               this.m22*m1.m21 + this.m23*m1.m31;
2669         m22 = this.m20*m1.m02 + this.m21*m1.m12 +
2670               this.m22*m1.m22 + this.m23*m1.m32;
2671         m23 = this.m20*m1.m03 + this.m21*m1.m13 +
2672               this.m22*m1.m23 + this.m23*m1.m33;
2673
2674         m30 = this.m30*m1.m00 + this.m31*m1.m10 +
2675               this.m32*m1.m20 + this.m33*m1.m30;
2676         m31 = this.m30*m1.m01 + this.m31*m1.m11 +
2677               this.m32*m1.m21 + this.m33*m1.m31;
2678         m32 = this.m30*m1.m02 + this.m31*m1.m12 +
2679               this.m32*m1.m22 + this.m33*m1.m32;
2680         m33 = this.m30*m1.m03 + this.m31*m1.m13 +
2681               this.m32*m1.m23 + this.m33*m1.m33;
2682
2683         this.m00 = m00; this.m01 = m01; this.m02 = m02; this.m03 = m03;
2684         this.m10 = m10; this.m11 = m11; this.m12 = m12; this.m13 = m13;
2685         this.m20 = m20; this.m21 = m21; this.m22 = m22; this.m23 = m23;
2686         this.m30 = m30; this.m31 = m31; this.m32 = m32; this.m33 = m33;
2687     }
2688
2689     /**
2690      * Sets the value of this matrix to the result of multiplying
2691      * the two argument matrices together.
2692      * @param m1 the first matrix
2693      * @param m2 the second matrix
2694      */
2695     public final void mul(Matrix4d m1, Matrix4d m2)
2696     {
2697         if (this != m1 && this != m2) {
2698             // code for mat mul
2699             this.m00 = m1.m00*m2.m00 + m1.m01*m2.m10 +
2700                        m1.m02*m2.m20 + m1.m03*m2.m30;
2701             this.m01 = m1.m00*m2.m01 + m1.m01*m2.m11 +
2702                        m1.m02*m2.m21 + m1.m03*m2.m31;
2703             this.m02 = m1.m00*m2.m02 + m1.m01*m2.m12 +
2704                        m1.m02*m2.m22 + m1.m03*m2.m32;
2705             this.m03 = m1.m00*m2.m03 + m1.m01*m2.m13 +
2706                        m1.m02*m2.m23 + m1.m03*m2.m33;
2707
2708             this.m10 = m1.m10*m2.m00 + m1.m11*m2.m10 +
2709                        m1.m12*m2.m20 + m1.m13*m2.m30;
2710             this.m11 = m1.m10*m2.m01 + m1.m11*m2.m11 +
2711                        m1.m12*m2.m21 + m1.m13*m2.m31;
2712             this.m12 = m1.m10*m2.m02 + m1.m11*m2.m12 +
2713                        m1.m12*m2.m22 + m1.m13*m2.m32;
2714             this.m13 = m1.m10*m2.m03 + m1.m11*m2.m13 +
2715                        m1.m12*m2.m23 + m1.m13*m2.m33;
2716
2717             this.m20 = m1.m20*m2.m00 + m1.m21*m2.m10 +
2718                        m1.m22*m2.m20 + m1.m23*m2.m30;
2719             this.m21 = m1.m20*m2.m01 + m1.m21*m2.m11 +
2720                        m1.m22*m2.m21 + m1.m23*m2.m31;
2721             this.m22 = m1.m20*m2.m02 + m1.m21*m2.m12 +
2722                        m1.m22*m2.m22 + m1.m23*m2.m32;
2723             this.m23 = m1.m20*m2.m03 + m1.m21*m2.m13 +
2724                        m1.m22*m2.m23 + m1.m23*m2.m33;
2725
2726             this.m30 = m1.m30*m2.m00 + m1.m31*m2.m10 +
2727                        m1.m32*m2.m20 + m1.m33*m2.m30;
2728             this.m31 = m1.m30*m2.m01 + m1.m31*m2.m11 +
2729                        m1.m32*m2.m21 + m1.m33*m2.m31;
2730             this.m32 = m1.m30*m2.m02 + m1.m31*m2.m12 +
2731                        m1.m32*m2.m22 + m1.m33*m2.m32;
2732             this.m33 = m1.m30*m2.m03 + m1.m31*m2.m13 +
2733                        m1.m32*m2.m23 + m1.m33*m2.m33;
2734         } else {
2735             double      m00, m01, m02, m03,
2736                         m10, m11, m12, m13,
2737                         m20, m21, m22, m23,
2738                         m30, m31, m32, m33;  // vars for temp result matrix
2739
2740             // code for mat mul
2741             m00 = m1.m00*m2.m00 + m1.m01*m2.m10 + m1.m02*m2.m20 + m1.m03*m2.m30;
2742             m01 = m1.m00*m2.m01 + m1.m01*m2.m11 + m1.m02*m2.m21 + m1.m03*m2.m31;
2743             m02 = m1.m00*m2.m02 + m1.m01*m2.m12 + m1.m02*m2.m22 + m1.m03*m2.m32;
2744             m03 = m1.m00*m2.m03 + m1.m01*m2.m13 + m1.m02*m2.m23 + m1.m03*m2.m33;
2745
2746             m10 = m1.m10*m2.m00 + m1.m11*m2.m10 + m1.m12*m2.m20 + m1.m13*m2.m30;
2747             m11 = m1.m10*m2.m01 + m1.m11*m2.m11 + m1.m12*m2.m21 + m1.m13*m2.m31;
2748             m12 = m1.m10*m2.m02 + m1.m11*m2.m12 + m1.m12*m2.m22 + m1.m13*m2.m32;
2749             m13 = m1.m10*m2.m03 + m1.m11*m2.m13 + m1.m12*m2.m23 + m1.m13*m2.m33;
2750
2751             m20 = m1.m20*m2.m00 + m1.m21*m2.m10 + m1.m22*m2.m20 + m1.m23*m2.m30;
2752             m21 = m1.m20*m2.m01 + m1.m21*m2.m11 + m1.m22*m2.m21 + m1.m23*m2.m31;
2753             m22 = m1.m20*m2.m02 + m1.m21*m2.m12 + m1.m22*m2.m22 + m1.m23*m2.m32;
2754             m23 = m1.m20*m2.m03 + m1.m21*m2.m13 + m1.m22*m2.m23 + m1.m23*m2.m33;
2755
2756             m30 = m1.m30*m2.m00 + m1.m31*m2.m10 + m1.m32*m2.m20 + m1.m33*m2.m30;
2757             m31 = m1.m30*m2.m01 + m1.m31*m2.m11 + m1.m32*m2.m21 + m1.m33*m2.m31;
2758             m32 = m1.m30*m2.m02 + m1.m31*m2.m12 + m1.m32*m2.m22 + m1.m33*m2.m32;
2759             m33 = m1.m30*m2.m03 + m1.m31*m2.m13 + m1.m32*m2.m23 + m1.m33*m2.m33;
2760
2761             this.m00 = m00; this.m01 = m01; this.m02 = m02; this.m03 = m03;
2762             this.m10 = m10; this.m11 = m11; this.m12 = m12; this.m13 = m13;
2763             this.m20 = m20; this.m21 = m21; this.m22 = m22; this.m23 = m23;
2764             this.m30 = m30; this.m31 = m31; this.m32 = m32; this.m33 = m33;
2765
2766         }
2767     }
2768
2769    /**
2770      *  Multiplies the transpose of matrix m1 times the transpose of matrix
2771      *  m2, and places the result into this.
2772      *  @param m1  the matrix on the left hand side of the multiplication
2773      *  @param m2  the matrix on the right hand side of the multiplication
2774      */
2775     public final void mulTransposeBoth(Matrix4d m1, Matrix4d m2)
2776     {
2777        if (this != m1 && this != m2) {
2778             this.m00 = m1.m00*m2.m00 + m1.m10*m2.m01 + m1.m20*m2.m02 + m1.m30*m2.m03;
2779             this.m01 = m1.m00*m2.m10 + m1.m10*m2.m11 + m1.m20*m2.m12 + m1.m30*m2.m13;
2780             this.m02 = m1.m00*m2.m20 + m1.m10*m2.m21 + m1.m20*m2.m22 + m1.m30*m2.m23;
2781             this.m03 = m1.m00*m2.m30 + m1.m10*m2.m31 + m1.m20*m2.m32 + m1.m30*m2.m33;
2782
2783             this.m10 = m1.m01*m2.m00 + m1.m11*m2.m01 + m1.m21*m2.m02 + m1.m31*m2.m03;
2784             this.m11 = m1.m01*m2.m10 + m1.m11*m2.m11 + m1.m21*m2.m12 + m1.m31*m2.m13;
2785             this.m12 = m1.m01*m2.m20 + m1.m11*m2.m21 + m1.m21*m2.m22 + m1.m31*m2.m23;
2786             this.m13 = m1.m01*m2.m30 + m1.m11*m2.m31 + m1.m21*m2.m32 + m1.m31*m2.m33;
2787
2788             this.m20 = m1.m02*m2.m00 + m1.m12*m2.m01 + m1.m22*m2.m02 + m1.m32*m2.m03;
2789             this.m21 = m1.m02*m2.m10 + m1.m12*m2.m11 + m1.m22*m2.m12 + m1.m32*m2.m13;
2790             this.m22 = m1.m02*m2.m20 + m1.m12*m2.m21 + m1.m22*m2.m22 + m1.m32*m2.m23;
2791             this.m23 = m1.m02*m2.m30 + m1.m12*m2.m31 + m1.m22*m2.m32 + m1.m32*m2.m33;
2792
2793             this.m30 = m1.m03*m2.m00 + m1.m13*m2.m01 + m1.m23*m2.m02 + m1.m33*m2.m03;
2794             this.m31 = m1.m03*m2.m10 + m1.m13*m2.m11 + m1.m23*m2.m12 + m1.m33*m2.m13;
2795             this.m32 = m1.m03*m2.m20 + m1.m13*m2.m21 + m1.m23*m2.m22 + m1.m33*m2.m23;
2796             this.m33 = m1.m03*m2.m30 + m1.m13*m2.m31 + m1.m23*m2.m32 + m1.m33*m2.m33;
2797         } else {
2798             double      m00, m01, m02, m03,
2799                         m10, m11, m12, m13,
2800                         m20, m21, m22, m23,  // vars for temp result matrix
2801                         m30, m31, m32, m33;
2802
2803             m00 = m1.m00*m2.m00 + m1.m10*m2.m01 + m1.m20*m2.m02 + m1.m30*m2.m03;
2804             m01 = m1.m00*m2.m10 + m1.m10*m2.m11 + m1.m20*m2.m12 + m1.m30*m2.m13;
2805             m02 = m1.m00*m2.m20 + m1.m10*m2.m21 + m1.m20*m2.m22 + m1.m30*m2.m23;
2806             m03 = m1.m00*m2.m30 + m1.m10*m2.m31 + m1.m20*m2.m32 + m1.m30*m2.m33;
2807
2808             m10 = m1.m01*m2.m00 + m1.m11*m2.m01 + m1.m21*m2.m02 + m1.m31*m2.m03;
2809             m11 = m1.m01*m2.m10 + m1.m11*m2.m11 + m1.m21*m2.m12 + m1.m31*m2.m13;
2810             m12 = m1.m01*m2.m20 + m1.m11*m2.m21 + m1.m21*m2.m22 + m1.m31*m2.m23;
2811             m13 = m1.m01*m2.m30 + m1.m11*m2.m31 + m1.m21*m2.m32 + m1.m31*m2.m33;
2812
2813             m20 = m1.m02*m2.m00 + m1.m12*m2.m01 + m1.m22*m2.m02 + m1.m32*m2.m03;
2814             m21 = m1.m02*m2.m10 + m1.m12*m2.m11 + m1.m22*m2.m12 + m1.m32*m2.m13;
2815             m22 = m1.m02*m2.m20 + m1.m12*m2.m21 + m1.m22*m2.m22 + m1.m32*m2.m23;
2816             m23 = m1.m02*m2.m30 + m1.m12*m2.m31 + m1.m22*m2.m32 + m1.m32*m2.m33;
2817
2818             m30 = m1.m03*m2.m00 + m1.m13*m2.m01 + m1.m23*m2.m02 + m1.m33*m2.m03;
2819             m31 = m1.m03*m2.m10 + m1.m13*m2.m11 + m1.m23*m2.m12 + m1.m33*m2.m13;
2820             m32 = m1.m03*m2.m20 + m1.m13*m2.m21 + m1.m23*m2.m22 + m1.m33*m2.m23;
2821             m33 = m1.m03*m2.m30 + m1.m13*m2.m31 + m1.m23*m2.m32 + m1.m33*m2.m33;
2822
2823             this.m00 = m00; this.m01 = m01; this.m02 = m02; this.m03 = m03;
2824             this.m10 = m10; this.m11 = m11; this.m12 = m12; this.m13 = m13;
2825             this.m20 = m20; this.m21 = m21; this.m22 = m22; this.m23 = m23;
2826             this.m30 = m30; this.m31 = m31; this.m32 = m32; this.m33 = m33;
2827         }
2828
2829     }
2830
2831
2832
2833    /**
2834      *  Multiplies matrix m1 times the transpose of matrix m2, and
2835      *  places the result into this.
2836      *  @param m1  the matrix on the left hand side of the multiplication
2837      *  @param m2  the matrix on the right hand side of the multiplication
2838      */
2839     public final void mulTransposeRight(Matrix4d m1, Matrix4d m2)
2840     {
2841     if (this != m1 && this != m2) {
2842       this.m00 = m1.m00*m2.m00 + m1.m01*m2.m01 + m1.m02*m2.m02 + m1.m03*m2.m03;
2843       this.m01 = m1.m00*m2.m10 + m1.m01*m2.m11 + m1.m02*m2.m12 + m1.m03*m2.m13;
2844       this.m02 = m1.m00*m2.m20 + m1.m01*m2.m21 + m1.m02*m2.m22 + m1.m03*m2.m23;
2845       this.m03 = m1.m00*m2.m30 + m1.m01*m2.m31 + m1.m02*m2.m32 + m1.m03*m2.m33;
2846
2847       this.m10 = m1.m10*m2.m00 + m1.m11*m2.m01 + m1.m12*m2.m02 + m1.m13*m2.m03;
2848       this.m11 = m1.m10*m2.m10 + m1.m11*m2.m11 + m1.m12*m2.m12 + m1.m13*m2.m13;
2849       this.m12 = m1.m10*m2.m20 + m1.m11*m2.m21 + m1.m12*m2.m22 + m1.m13*m2.m23;
2850       this.m13 = m1.m10*m2.m30 + m1.m11*m2.m31 + m1.m12*m2.m32 + m1.m13*m2.m33;
2851
2852       this.m20 = m1.m20*m2.m00 + m1.m21*m2.m01 + m1.m22*m2.m02 + m1.m23*m2.m03;
2853       this.m21 = m1.m20*m2.m10 + m1.m21*m2.m11 + m1.m22*m2.m12 + m1.m23*m2.m13;
2854       this.m22 = m1.m20*m2.m20 + m1.m21*m2.m21 + m1.m22*m2.m22 + m1.m23*m2.m23;
2855       this.m23 = m1.m20*m2.m30 + m1.m21*m2.m31 + m1.m22*m2.m32 + m1.m23*m2.m33;
2856
2857       this.m30 = m1.m30*m2.m00 + m1.m31*m2.m01 + m1.m32*m2.m02 + m1.m33*m2.m03;
2858       this.m31 = m1.m30*m2.m10 + m1.m31*m2.m11 + m1.m32*m2.m12 + m1.m33*m2.m13;
2859       this.m32 = m1.m30*m2.m20 + m1.m31*m2.m21 + m1.m32*m2.m22 + m1.m33*m2.m23;
2860       this.m33 = m1.m30*m2.m30 + m1.m31*m2.m31 + m1.m32*m2.m32 + m1.m33*m2.m33;
2861     } else {
2862             double      m00, m01, m02, m03,
2863                         m10, m11, m12, m13,
2864                         m20, m21, m22, m23,  // vars for temp result matrix
2865                         m30, m31, m32, m33;
2866
2867       m00 = m1.m00*m2.m00 + m1.m01*m2.m01 + m1.m02*m2.m02 + m1.m03*m2.m03;
2868       m01 = m1.m00*m2.m10 + m1.m01*m2.m11 + m1.m02*m2.m12 + m1.m03*m2.m13;
2869       m02 = m1.m00*m2.m20 + m1.m01*m2.m21 + m1.m02*m2.m22 + m1.m03*m2.m23;
2870       m03 = m1.m00*m2.m30 + m1.m01*m2.m31 + m1.m02*m2.m32 + m1.m03*m2.m33;
2871
2872       m10 = m1.m10*m2.m00 + m1.m11*m2.m01 + m1.m12*m2.m02 + m1.m13*m2.m03;
2873       m11 = m1.m10*m2.m10 + m1.m11*m2.m11 + m1.m12*m2.m12 + m1.m13*m2.m13;
2874       m12 = m1.m10*m2.m20 + m1.m11*m2.m21 + m1.m12*m2.m22 + m1.m13*m2.m23;
2875       m13 = m1.m10*m2.m30 + m1.m11*m2.m31 + m1.m12*m2.m32 + m1.m13*m2.m33;
2876
2877       m20 = m1.m20*m2.m00 + m1.m21*m2.m01 + m1.m22*m2.m02 + m1.m23*m2.m03;
2878       m21 = m1.m20*m2.m10 + m1.m21*m2.m11 + m1.m22*m2.m12 + m1.m23*m2.m13;
2879       m22 = m1.m20*m2.m20 + m1.m21*m2.m21 + m1.m22*m2.m22 + m1.m23*m2.m23;
2880       m23 = m1.m20*m2.m30 + m1.m21*m2.m31 + m1.m22*m2.m32 + m1.m23*m2.m33;
2881
2882       m30 = m1.m30*m2.m00 + m1.m31*m2.m01 + m1.m32*m2.m02 + m1.m33*m2.m03;
2883       m31 = m1.m30*m2.m10 + m1.m31*m2.m11 + m1.m32*m2.m12 + m1.m33*m2.m13;
2884       m32 = m1.m30*m2.m20 + m1.m31*m2.m21 + m1.m32*m2.m22 + m1.m33*m2.m23;
2885       m33 = m1.m30*m2.m30 + m1.m31*m2.m31 + m1.m32*m2.m32 + m1.m33*m2.m33;
2886
2887       this.m00 = m00; this.m01 = m01; this.m02 = m02; this.m03 = m03;
2888       this.m10 = m10; this.m11 = m11; this.m12 = m12; this.m13 = m13;
2889       this.m20 = m20; this.m21 = m21; this.m22 = m22; this.m23 = m23;
2890       this.m30 = m30; this.m31 = m31; this.m32 = m32; this.m33 = m33;
2891     }
2892 }
2893
2894
2895    /**
2896      *  Multiplies the transpose of matrix m1 times matrix m2, and
2897      *  places the result into this.
2898      *  @param m1  the matrix on the left hand side of the multiplication
2899      *  @param m2  the matrix on the right hand side of the multiplication
2900      */
2901     public final void mulTransposeLeft(Matrix4d m1, Matrix4d m2)
2902     {
2903     if (this != m1 && this != m2) {
2904       this.m00 = m1.m00*m2.m00 + m1.m10*m2.m10 + m1.m20*m2.m20 + m1.m30*m2.m30;
2905       this.m01 = m1.m00*m2.m01 + m1.m10*m2.m11 + m1.m20*m2.m21 + m1.m30*m2.m31;
2906       this.m02 = m1.m00*m2.m02 + m1.m10*m2.m12 + m1.m20*m2.m22 + m1.m30*m2.m32;
2907       this.m03 = m1.m00*m2.m03 + m1.m10*m2.m13 + m1.m20*m2.m23 + m1.m30*m2.m33;
2908
2909       this.m10 = m1.m01*m2.m00 + m1.m11*m2.m10 + m1.m21*m2.m20 + m1.m31*m2.m30;
2910       this.m11 = m1.m01*m2.m01 + m1.m11*m2.m11 + m1.m21*m2.m21 + m1.m31*m2.m31;
2911       this.m12 = m1.m01*m2.m02 + m1.m11*m2.m12 + m1.m21*m2.m22 + m1.m31*m2.m32;
2912       this.m13 = m1.m01*m2.m03 + m1.m11*m2.m13 + m1.m21*m2.m23 + m1.m31*m2.m33;
2913
2914       this.m20 = m1.m02*m2.m00 + m1.m12*m2.m10 + m1.m22*m2.m20 + m1.m32*m2.m30;
2915       this.m21 = m1.m02*m2.m01 + m1.m12*m2.m11 + m1.m22*m2.m21 + m1.m32*m2.m31;
2916       this.m22 = m1.m02*m2.m02 + m1.m12*m2.m12 + m1.m22*m2.m22 + m1.m32*m2.m32;
2917       this.m23 = m1.m02*m2.m03 + m1.m12*m2.m13 + m1.m22*m2.m23 + m1.m32*m2.m33;
2918
2919       this.m30 = m1.m03*m2.m00 + m1.m13*m2.m10 + m1.m23*m2.m20 + m1.m33*m2.m30;
2920       this.m31 = m1.m03*m2.m01 + m1.m13*m2.m11 + m1.m23*m2.m21 + m1.m33*m2.m31;
2921       this.m32 = m1.m03*m2.m02 + m1.m13*m2.m12 + m1.m23*m2.m22 + m1.m33*m2.m32;
2922       this.m33 = m1.m03*m2.m03 + m1.m13*m2.m13 + m1.m23*m2.m23 + m1.m33*m2.m33;
2923     } else {
2924             double      m00, m01, m02, m03,
2925                         m10, m11, m12, m13,
2926                         m20, m21, m22, m23,  // vars for temp result matrix
2927                         m30, m31, m32, m33;
2928
2929
2930
2931       m00 = m1.m00*m2.m00 + m1.m10*m2.m10 + m1.m20*m2.m20 + m1.m30*m2.m30;
2932       m01 = m1.m00*m2.m01 + m1.m10*m2.m11 + m1.m20*m2.m21 + m1.m30*m2.m31;
2933       m02 = m1.m00*m2.m02 + m1.m10*m2.m12 + m1.m20*m2.m22 + m1.m30*m2.m32;
2934       m03 = m1.m00*m2.m03 + m1.m10*m2.m13 + m1.m20*m2.m23 + m1.m30*m2.m33;
2935
2936       m10 = m1.m01*m2.m00 + m1.m11*m2.m10 + m1.m21*m2.m20 + m1.m31*m2.m30;
2937       m11 = m1.m01*m2.m01 + m1.m11*m2.m11 + m1.m21*m2.m21 + m1.m31*m2.m31;
2938       m12 = m1.m01*m2.m02 + m1.m11*m2.m12 + m1.m21*m2.m22 + m1.m31*m2.m32;
2939       m13 = m1.m01*m2.m03 + m1.m11*m2.m13 + m1.m21*m2.m23 + m1.m31*m2.m33;
2940
2941       m20 = m1.m02*m2.m00 + m1.m12*m2.m10 + m1.m22*m2.m20 + m1.m32*m2.m30;
2942       m21 = m1.m02*m2.m01 + m1.m12*m2.m11 + m1.m22*m2.m21 + m1.m32*m2.m31;
2943       m22 = m1.m02*m2.m02 + m1.m12*m2.m12 + m1.m22*m2.m22 + m1.m32*m2.m32;
2944       m23 = m1.m02*m2.m03 + m1.m12*m2.m13 + m1.m22*m2.m23 + m1.m32*m2.m33;
2945
2946       m30 = m1.m03*m2.m00 + m1.m13*m2.m10 + m1.m23*m2.m20 + m1.m33*m2.m30;
2947       m31 = m1.m03*m2.m01 + m1.m13*m2.m11 + m1.m23*m2.m21 + m1.m33*m2.m31;
2948       m32 = m1.m03*m2.m02 + m1.m13*m2.m12 + m1.m23*m2.m22 + m1.m33*m2.m32;
2949       m33 = m1.m03*m2.m03 + m1.m13*m2.m13 + m1.m23*m2.m23 + m1.m33*m2.m33;
2950
2951       this.m00 = m00; this.m01 = m01; this.m02 = m02; this.m03 = m03;
2952       this.m10 = m10; this.m11 = m11; this.m12 = m12; this.m13 = m13;
2953       this.m20 = m20; this.m21 = m21; this.m22 = m22; this.m23 = m23;
2954       this.m30 = m30; this.m31 = m31; this.m32 = m32; this.m33 = m33;
2955     }
2956
2957     }
2958
2959
2960    /**
2961      * Returns true if all of the data members of Matrix4d m1 are
2962      * equal to the corresponding data members in this Matrix4d.
2963      * @param m1  the matrix with which the comparison is made
2964      * @return  true or false
2965      */
2966     public boolean equals(Matrix4d m1)
2967     {
2968       try {
2969         return(this.m00 == m1.m00 && this.m01 == m1.m01 && this.m02 == m1.m02
2970             && this.m03 == m1.m03 && this.m10 == m1.m10 && this.m11 == m1.m11
2971             && this.m12 == m1.m12 && this.m13 == m1.m13 && this.m20 == m1.m20
2972             && this.m21 == m1.m21 && this.m22 == m1.m22 && this.m23 == m1.m23
2973             && this.m30 == m1.m30 && this.m31 == m1.m31 && this.m32 == m1.m32
2974             && this.m33 == m1.m33);
2975       }
2976       catch (NullPointerException e2) { return false; }
2977
2978     }
2979
2980    /**
2981      * Returns true if the Object t1 is of type Matrix4d and all of the
2982      * data members of t1 are equal to the corresponding data members in
2983      * this Matrix4d.
2984      * @param t1  the matrix with which the comparison is made
2985      * @return  true or false
2986      */
2987     public boolean equals(Object t1)
2988     {
2989         try {
2990            Matrix4d m2 = (Matrix4d) t1;
2991            return(this.m00 == m2.m00 && this.m01 == m2.m01 && this.m02 == m2.m02
2992              && this.m03 == m2.m03 && this.m10 == m2.m10 && this.m11 == m2.m11
2993              && this.m12 == m2.m12 && this.m13 == m2.m13 && this.m20 == m2.m20
2994              && this.m21 == m2.m21 && this.m22 == m2.m22 && this.m23 == m2.m23
2995              && this.m30 == m2.m30 && this.m31 == m2.m31 && this.m32 == m2.m32
2996              && this.m33 == m2.m33);
2997         }
2998         catch (ClassCastException   e1) { return false; }
2999         catch (NullPointerException e2) { return false; }
3000     }
3001
3002     /**
3003      * @deprecated Use epsilonEquals(Matrix4d,double) instead
3004      */
3005     public boolean epsilonEquals(Matrix4d m1, float epsilon) {
3006         return epsilonEquals(m1, (double)epsilon);
3007     }
3008
3009     /**
3010      * Returns true if the L-infinite distance between this matrix
3011      * and matrix m1 is less than or equal to the epsilon parameter,
3012      * otherwise returns false.  The L-infinite
3013      * distance is equal to
3014      * MAX[i=0,1,2,3 ; j=0,1,2,3 ; abs(this.m(i,j) - m1.m(i,j)]
3015      * @param m1  the matrix to be compared to this matrix
3016      * @param epsilon  the threshold value
3017      */
3018     public boolean epsilonEquals(Matrix4d m1, double epsilon) {
3019        double diff;
3020
3021        diff = m00 - m1.m00;
3022        if((diff<0?-diff:diff) > epsilon) return false;
3023
3024        diff = m01 - m1.m01;
3025        if((diff<0?-diff:diff) > epsilon) return false;
3026
3027        diff = m02 - m1.m02;
3028        if((diff<0?-diff:diff) > epsilon) return false;
3029
3030        diff = m03 - m1.m03;
3031        if((diff<0?-diff:diff) > epsilon) return false;
3032
3033        diff = m10 - m1.m10;
3034        if((diff<0?-diff:diff) > epsilon) return false;
3035
3036        diff = m11 - m1.m11;
3037        if((diff<0?-diff:diff) > epsilon) return false;
3038
3039        diff = m12 - m1.m12;
3040        if((diff<0?-diff:diff) > epsilon) return false;
3041
3042        diff = m13 - m1.m13;
3043        if((diff<0?-diff:diff) > epsilon) return false;
3044
3045        diff = m20 - m1.m20;
3046        if((diff<0?-diff:diff) > epsilon) return false;
3047
3048        diff = m21 - m1.m21;
3049        if((diff<0?-diff:diff) > epsilon) return false;
3050
3051        diff = m22 - m1.m22;
3052        if((diff<0?-diff:diff) > epsilon) return false;
3053
3054        diff = m23 - m1.m23;
3055        if((diff<0?-diff:diff) > epsilon) return false;
3056
3057        diff = m30 - m1.m30;
3058        if((diff<0?-diff:diff) > epsilon) return false;
3059
3060        diff = m31 - m1.m31;
3061        if((diff<0?-diff:diff) > epsilon) return false;
3062
3063        diff = m32 - m1.m32;
3064        if((diff<0?-diff:diff) > epsilon) return false;
3065
3066        diff = m33 - m1.m33;
3067        if((diff<0?-diff:diff) > epsilon) return false;
3068
3069        return true;
3070     }
3071
3072     /**
3073      * Returns a hash code value based on the data values in this
3074      * object.  Two different Matrix4d objects with identical data values
3075      * (i.e., Matrix4d.equals returns true) will return the same hash
3076      * code value.  Two objects with different data members may return the
3077      * same hash value, although this is not likely.
3078      * @return the integer hash code value
3079      */
3080     public int hashCode() {
3081         long bits = 1L;
3082         bits = 31L * bits + VecMathUtil.doubleToLongBits(m00);
3083         bits = 31L * bits + VecMathUtil.doubleToLongBits(m01);
3084         bits = 31L * bits + VecMathUtil.doubleToLongBits(m02);
3085         bits = 31L * bits + VecMathUtil.doubleToLongBits(m03);
3086         bits = 31L * bits + VecMathUtil.doubleToLongBits(m10);
3087         bits = 31L * bits + VecMathUtil.doubleToLongBits(m11);
3088         bits = 31L * bits + VecMathUtil.doubleToLongBits(m12);
3089         bits = 31L * bits + VecMathUtil.doubleToLongBits(m13);
3090         bits = 31L * bits + VecMathUtil.doubleToLongBits(m20);
3091         bits = 31L * bits + VecMathUtil.doubleToLongBits(m21);
3092         bits = 31L * bits + VecMathUtil.doubleToLongBits(m22);
3093         bits = 31L * bits + VecMathUtil.doubleToLongBits(m23);
3094         bits = 31L * bits + VecMathUtil.doubleToLongBits(m30);
3095         bits = 31L * bits + VecMathUtil.doubleToLongBits(m31);
3096         bits = 31L * bits + VecMathUtil.doubleToLongBits(m32);
3097         bits = 31L * bits + VecMathUtil.doubleToLongBits(m33);
3098         return (int) (bits ^ (bits >> 32));
3099     }
3100
3101
3102   /**
3103    * Transform the vector vec using this Matrix4d and place the
3104    * result into vecOut.
3105    * @param vec  the double precision vector to be transformed
3106    * @param vecOut  the vector into which the transformed values are placed
3107    */
3108     public final void transform(Tuple4d vec, Tuple4d vecOut)
3109     {
3110            double x,y,z,w;
3111            x = (m00*vec.x + m01*vec.y
3112                   + m02*vec.z + m03*vec.w);
3113            y = (m10*vec.x + m11*vec.y
3114                    + m12*vec.z + m13*vec.w);
3115            z = (m20*vec.x + m21*vec.y
3116                      + m22*vec.z + m23*vec.w);
3117            vecOut.w = (m30*vec.x + m31*vec.y
3118                       + m32*vec.z + m33*vec.w);
3119            vecOut.x = x;
3120            vecOut.y = y;
3121            vecOut.z = z;
3122     }
3123
3124   /**
3125    * Transform the vector vec using this Matrix4d and place the
3126    * result back into vec.
3127    * @param vec  the double precision vector to be transformed
3128    */
3129     public final void transform(Tuple4d vec)
3130     {
3131          double x,y,z;
3132
3133            x = (m00*vec.x + m01*vec.y
3134                      + m02*vec.z + m03*vec.w);
3135            y = (m10*vec.x + m11*vec.y
3136                      + m12*vec.z + m13*vec.w);
3137            z = (m20*vec.x + m21*vec.y
3138                      + m22*vec.z + m23*vec.w);
3139            vec.w = (m30*vec.x + m31*vec.y
3140                       + m32*vec.z + m33*vec.w);
3141            vec.x = x;
3142            vec.y = y;
3143            vec.z = z;
3144     }
3145
3146   /**
3147    * Transform the vector vec using this Matrix4d and place the
3148    * result into vecOut.
3149    * @param vec  the single precision vector to be transformed
3150    * @param vecOut  the vector into which the transformed values are placed
3151    */
3152     public final void transform(Tuple4f vec, Tuple4f vecOut)
3153     {
3154         float x,y,z;
3155            x = (float) (m00*vec.x + m01*vec.y
3156                      + m02*vec.z + m03*vec.w);
3157            y = (float) (m10*vec.x + m11*vec.y
3158                      + m12*vec.z + m13*vec.w);
3159            z = (float) (m20*vec.x + m21*vec.y
3160                      + m22*vec.z + m23*vec.w);
3161            vecOut.w = (float) (m30*vec.x + m31*vec.y
3162                       + m32*vec.z + m33*vec.w);
3163            vecOut.x = x;
3164            vecOut.y = y;
3165            vecOut.z = z;
3166     }
3167
3168   /**
3169    * Transform the vector vec using this Transform and place the
3170    * result back into vec.
3171    * @param vec  the single precision vector to be transformed
3172    */
3173     public final void transform(Tuple4f vec)
3174     {
3175          float x,y,z;
3176
3177            x = (float) (m00*vec.x + m01*vec.y
3178                      + m02*vec.z + m03*vec.w);
3179            y = (float) (m10*vec.x + m11*vec.y
3180                      + m12*vec.z + m13*vec.w);
3181            z = (float) (m20*vec.x + m21*vec.y
3182                      + m22*vec.z + m23*vec.w);
3183            vec.w = (float) (m30*vec.x + m31*vec.y
3184                       + m32*vec.z + m33*vec.w);
3185            vec.x = x;
3186            vec.y = y;
3187            vec.z = z;
3188     }
3189
3190
3191   /**
3192    * Transforms the point parameter with this Matrix4d and
3193    * places the result into pointOut.  The fourth element of the
3194    * point input parameter is assumed to be one.
3195    * @param point  the input point to be transformed.
3196    * @param pointOut  the transformed point
3197    */
3198     public final void transform(Point3d point, Point3d pointOut)
3199     {
3200         double x,y;
3201         x = m00*point.x + m01*point.y + m02*point.z + m03;
3202         y = m10*point.x + m11*point.y + m12*point.z + m13;
3203         pointOut.z = m20*point.x + m21*point.y + m22*point.z + m23;
3204         pointOut.x = x;
3205         pointOut.y = y;
3206
3207     }
3208
3209
3210   /**
3211    * Transforms the point parameter with this Matrix4d and
3212    * places the result back into point.  The fourth element of the
3213    * point input parameter is assumed to be one.
3214    * @param point  the input point to be transformed.
3215    */
3216     public final void transform(Point3d point)
3217     {
3218         double x, y;
3219         x = m00*point.x + m01*point.y + m02*point.z + m03;
3220         y = m10*point.x + m11*point.y + m12*point.z + m13;
3221         point.z =  m20*point.x + m21*point.y + m22*point.z + m23;
3222         point.x = x;
3223         point.y = y;
3224     }
3225
3226
3227   /**
3228    * Transforms the point parameter with this Matrix4d and
3229    * places the result into pointOut.  The fourth element of the
3230    * point input parameter is assumed to be one.
3231    * @param point  the input point to be transformed.
3232    * @param pointOut  the transformed point
3233    */
3234     public final void transform(Point3f point, Point3f pointOut)
3235     {
3236         float x,y;
3237
3238         x = (float) (m00*point.x + m01*point.y + m02*point.z + m03);
3239         y = (float) (m10*point.x + m11*point.y + m12*point.z + m13);
3240         pointOut.z = (float) (m20*point.x + m21*point.y + m22*point.z + m23);
3241         pointOut.x = x;
3242         pointOut.y = y;
3243     }
3244
3245
3246   /**
3247    * Transforms the point parameter with this Matrix4d and
3248    * places the result back into point.  The fourth element of the
3249    * point input parameter is assumed to be one.
3250    * @param point  the input point to be transformed.
3251    */
3252     public final void transform(Point3f point)
3253     {
3254         float  x, y;
3255         x = (float) (m00*point.x + m01*point.y + m02*point.z + m03);
3256         y = (float) (m10*point.x + m11*point.y + m12*point.z + m13);
3257         point.z =  (float) (m20*point.x + m21*point.y + m22*point.z + m23);
3258         point.x = x;
3259         point.y = y;
3260     }
3261
3262
3263   /**
3264    * Transforms the normal parameter by this Matrix4d and places the value
3265    * into normalOut.  The fourth element of the normal is assumed to be zero.
3266    * @param normal   the input normal to be transformed.
3267    * @param normalOut  the transformed normal
3268    */
3269     public final void transform(Vector3d normal, Vector3d normalOut)
3270     {
3271         double x,y;
3272         x =  m00*normal.x + m01*normal.y + m02*normal.z;
3273         y =  m10*normal.x + m11*normal.y + m12*normal.z;
3274         normalOut.z =  m20*normal.x + m21*normal.y + m22*normal.z;
3275         normalOut.x = x;
3276         normalOut.y = y;
3277     }
3278
3279
3280   /**
3281    * Transforms the normal parameter by this transform and places the value
3282    * back into normal.  The fourth element of the normal is assumed to be zero.
3283    * @param normal   the input normal to be transformed.
3284    */
3285     public final void transform(Vector3d normal)
3286     {
3287         double x, y;
3288
3289         x =  m00*normal.x + m01*normal.y + m02*normal.z;
3290         y =  m10*normal.x + m11*normal.y + m12*normal.z;
3291         normal.z =  m20*normal.x + m21*normal.y + m22*normal.z;
3292         normal.x = x;
3293         normal.y = y;
3294     }
3295
3296
3297   /**
3298    * Transforms the normal parameter by this Matrix4d and places the value
3299    * into normalOut.  The fourth element of the normal is assumed to be zero.
3300    * @param normal   the input normal to be transformed.
3301    * @param normalOut  the transformed normal
3302    */
3303     public final void transform(Vector3f normal, Vector3f normalOut)
3304     {
3305         float x,y;
3306         x =  (float) (m00*normal.x + m01*normal.y + m02*normal.z);
3307         y =  (float) (m10*normal.x + m11*normal.y + m12*normal.z);
3308         normalOut.z =  (float) (m20*normal.x + m21*normal.y + m22*normal.z);
3309         normalOut.x = x;
3310         normalOut.y = y;
3311     }
3312
3313
3314   /**
3315    * Transforms the normal parameter by this transform and places the value
3316    * back into normal.  The fourth element of the normal is assumed to be zero.
3317    * @param normal   the input normal to be transformed.
3318    */
3319     public final void transform(Vector3f normal)
3320     {
3321         float x, y;
3322
3323         x =  (float) (m00*normal.x + m01*normal.y + m02*normal.z);
3324         y =  (float) (m10*normal.x + m11*normal.y + m12*normal.z);
3325         normal.z =  (float) (m20*normal.x + m21*normal.y + m22*normal.z);
3326         normal.x = x;
3327         normal.y = y;
3328     }
3329
3330    /**
3331      * Sets the rotational component (upper 3x3) of this matrix to the
3332      * matrix values in the double precision Matrix3d argument; the other
3333      * elements of this matrix are unchanged; a singular value
3334      * decomposition is performed on this object's upper 3x3 matrix to
3335      * factor out the scale, then this object's upper 3x3 matrix components
3336      * are replaced by the passed rotation components,
3337      * and then the scale is reapplied to the rotational components.
3338      * @param m1   double precision 3x3 matrix
3339      */
3340     public final void setRotation( Matrix3d m1){
3341         double[]    tmp_rot = new double[9];  // scratch matrix
3342         double[]    tmp_scale = new double[3];  // scratch matrix
3343
3344         getScaleRotate( tmp_scale, tmp_rot );
3345
3346         m00 = m1.m00*tmp_scale[0];
3347         m01 = m1.m01*tmp_scale[1];
3348         m02 = m1.m02*tmp_scale[2];
3349
3350         m10 = m1.m10*tmp_scale[0];
3351         m11 = m1.m11*tmp_scale[1];
3352         m12 = m1.m12*tmp_scale[2];
3353
3354         m20 = m1.m20*tmp_scale[0];
3355         m21 = m1.m21*tmp_scale[1];
3356         m22 = m1.m22*tmp_scale[2];
3357
3358     }
3359
3360
3361    /**
3362      * Sets the rotational component (upper 3x3) of this matrix to the
3363      * matrix values in the single precision Matrix3f argument; the other
3364      * elements of this matrix are unchanged; a singular value
3365      * decomposition is performed on this object's upper 3x3 matrix to
3366      * factor out the scale, then this object's upper 3x3 matrix components
3367      * are replaced by the passed rotation components,
3368      * and then the scale is reapplied to the rotational components.
3369      * @param m1   single precision 3x3 matrix
3370      */
3371    public final void setRotation( Matrix3f m1)
3372    {
3373
3374        double[]    tmp_rot = new double[9];  // scratch matrix
3375        double[]    tmp_scale = new double[3];  // scratch matrix
3376        getScaleRotate( tmp_scale, tmp_rot );
3377
3378         m00 = m1.m00*tmp_scale[0];
3379         m01 = m1.m01*tmp_scale[1];
3380         m02 = m1.m02*tmp_scale[2];
3381
3382         m10 = m1.m10*tmp_scale[0];
3383         m11 = m1.m11*tmp_scale[1];
3384         m12 = m1.m12*tmp_scale[2];
3385
3386         m20 = m1.m20*tmp_scale[0];
3387         m21 = m1.m21*tmp_scale[1];
3388         m22 = m1.m22*tmp_scale[2];
3389    }
3390
3391    /**
3392      * Sets the rotational component (upper 3x3) of this matrix to the
3393      * matrix equivalent values of the quaternion argument; the other
3394      * elements of this matrix are unchanged; a singular value
3395      * decomposition is performed on this object's upper 3x3 matrix to
3396      * factor out the scale, then this object's upper 3x3 matrix components
3397      * are replaced by the matrix equivalent of the quaternion,
3398      * and then the scale is reapplied to the rotational components.
3399      * @param q1    the quaternion that specifies the rotation
3400      */
3401    public final void setRotation(Quat4f q1){
3402        double[]    tmp_rot = new double[9];  // scratch matrix
3403        double[]    tmp_scale = new double[3];  // scratch matrix
3404        getScaleRotate( tmp_scale, tmp_rot );
3405
3406         m00 = (1.0 - 2.0f*q1.y*q1.y - 2.0f*q1.z*q1.z)*tmp_scale[0];
3407         m10 = (2.0*(q1.x*q1.y + q1.w*q1.z))*tmp_scale[0];
3408         m20 = (2.0*(q1.x*q1.z - q1.w*q1.y))*tmp_scale[0];
3409
3410         m01 = (2.0*(q1.x*q1.y - q1.w*q1.z))*tmp_scale[1];
3411         m11 = (1.0 - 2.0f*q1.x*q1.x - 2.0f*q1.z*q1.z)*tmp_scale[1];
3412         m21 = (2.0*(q1.y*q1.z + q1.w*q1.x))*tmp_scale[1];
3413
3414         m02 = (2.0*(q1.x*q1.z + q1.w*q1.y))*tmp_scale[2];
3415         m12 = (2.0*(q1.y*q1.z - q1.w*q1.x))*tmp_scale[2];
3416         m22 = (1.0 - 2.0f*q1.x*q1.x - 2.0f*q1.y*q1.y)*tmp_scale[2];
3417
3418       }
3419
3420
3421    /**
3422      * Sets the rotational component (upper 3x3) of this matrix to the
3423      * matrix equivalent values of the quaternion argument; the other
3424      * elements of this matrix are unchanged; a singular value
3425      * decomposition is performed on this object's upper 3x3 matrix to
3426      * factor out the scale, then this object's upper 3x3 matrix components
3427      * are replaced by the matrix equivalent of the quaternion,
3428      * and then the scale is reapplied to the rotational components.
3429      * @param q1    the quaternion that specifies the rotation
3430      */
3431    public final void setRotation(Quat4d q1){
3432
3433        double[]    tmp_rot = new double[9];  // scratch matrix
3434        double[]    tmp_scale = new double[3];  // scratch matrix
3435        getScaleRotate( tmp_scale, tmp_rot );
3436
3437         m00 = (1.0 - 2.0f*q1.y*q1.y - 2.0f*q1.z*q1.z)*tmp_scale[0];
3438         m10 = (2.0*(q1.x*q1.y + q1.w*q1.z))*tmp_scale[0];
3439         m20 = (2.0*(q1.x*q1.z - q1.w*q1.y))*tmp_scale[0];
3440
3441         m01 = (2.0*(q1.x*q1.y - q1.w*q1.z))*tmp_scale[1];
3442         m11 = (1.0 - 2.0f*q1.x*q1.x - 2.0f*q1.z*q1.z)*tmp_scale[1];
3443         m21 = (2.0*(q1.y*q1.z + q1.w*q1.x))*tmp_scale[1];
3444
3445         m02 = (2.0*(q1.x*q1.z + q1.w*q1.y))*tmp_scale[2];
3446         m12 = (2.0*(q1.y*q1.z - q1.w*q1.x))*tmp_scale[2];
3447         m22 = (1.0 - 2.0f*q1.x*q1.x - 2.0f*q1.y*q1.y)*tmp_scale[2];
3448
3449       }
3450
3451    /**
3452      * Sets the rotational component (upper 3x3) of this matrix to the
3453      * matrix equivalent values of the axis-angle argument; the other
3454      * elements of this matrix are unchanged; a singular value
3455      * decomposition is performed on this object's upper 3x3 matrix to
3456      * factor out the scale, then this object's upper 3x3 matrix components
3457      * are replaced by the matrix equivalent of the axis-angle,
3458      * and then the scale is reapplied to the rotational components.
3459      * @param a1 the axis-angle to be converted (x, y, z, angle)
3460      */
3461     public final void setRotation(AxisAngle4d a1)
3462     {
3463         double[]    tmp_rot = new double[9];  // scratch matrix
3464         double[]    tmp_scale = new double[3];  // scratch matrix
3465
3466         getScaleRotate( tmp_scale, tmp_rot );
3467
3468         double mag = 1.0/Math.sqrt( a1.x*a1.x + a1.y*a1.y + a1.z*a1.z);
3469         double ax = a1.x*mag;
3470         double ay = a1.y*mag;
3471         double az = a1.z*mag;
3472
3473         double sinTheta = Math.sin(a1.angle);
3474         double cosTheta = Math.cos(a1.angle);
3475         double t = 1.0 - cosTheta;
3476
3477         double xz = a1.x * a1.z;
3478         double xy = a1.x * a1.y;
3479         double yz = a1.y * a1.z;
3480
3481         m00 = (t * ax * ax + cosTheta)*tmp_scale[0];
3482         m01 = (t * xy - sinTheta * az)*tmp_scale[1];
3483         m02 = (t * xz + sinTheta * ay)*tmp_scale[2];
3484
3485         m10 = (t * xy + sinTheta * az)*tmp_scale[0];
3486         m11 = (t * ay * ay + cosTheta)*tmp_scale[1];
3487         m12 = (t * yz - sinTheta * ax)*tmp_scale[2];
3488
3489         m20 = (t * xz - sinTheta * ay)*tmp_scale[0];
3490         m21 = (t * yz + sinTheta * ax)*tmp_scale[1];
3491         m22 = (t * az * az + cosTheta)*tmp_scale[2];
3492
3493     }
3494
3495   /**
3496     *  Sets this matrix to all zeros.
3497     */
3498    public final void setZero()
3499    {
3500         m00 = 0.0;
3501         m01 = 0.0;
3502         m02 = 0.0;
3503         m03 = 0.0;
3504         m10 = 0.0;
3505         m11 = 0.0;
3506         m12 = 0.0;
3507         m13 = 0.0;
3508         m20 = 0.0;
3509         m21 = 0.0;
3510         m22 = 0.0;
3511         m23 = 0.0;
3512         m30 = 0.0;
3513         m31 = 0.0;
3514         m32 = 0.0;
3515         m33 = 0.0;
3516    }
3517
3518    /**
3519      * Negates the value of this matrix: this = -this.
3520      */
3521     public final void negate()
3522     {
3523         m00 = -m00;
3524         m01 = -m01;
3525         m02 = -m02;
3526         m03 = -m03;
3527         m10 = -m10;
3528         m11 = -m11;
3529         m12 = -m12;
3530         m13 = -m13;
3531         m20 = -m20;
3532         m21 = -m21;
3533         m22 = -m22;
3534         m23 = -m23;
3535         m30 = -m30;
3536         m31 = -m31;
3537         m32 = -m32;
3538         m33 = -m33;
3539     }
3540
3541    /**
3542      *  Sets the value of this matrix equal to the negation of
3543      *  of the Matrix4d parameter.
3544      *  @param m1  the source matrix
3545      */
3546     public final void negate(Matrix4d m1)
3547     {
3548         this.m00 = -m1.m00;
3549         this.m01 = -m1.m01;
3550         this.m02 = -m1.m02;
3551         this.m03 = -m1.m03;
3552         this.m10 = -m1.m10;
3553         this.m11 = -m1.m11;
3554         this.m12 = -m1.m12;
3555         this.m13 = -m1.m13;
3556         this.m20 = -m1.m20;
3557         this.m21 = -m1.m21;
3558         this.m22 = -m1.m22;
3559         this.m23 = -m1.m23;
3560         this.m30 = -m1.m30;
3561         this.m31 = -m1.m31;
3562         this.m32 = -m1.m32;
3563         this.m33 = -m1.m33;
3564     }
3565     private final void getScaleRotate(double scales[], double rots[]) {
3566         double[]    tmp = new double[9];  // scratch matrix
3567         tmp[0] = m00;
3568         tmp[1] = m01;
3569         tmp[2] = m02;
3570
3571         tmp[3] = m10;
3572         tmp[4] = m11;
3573         tmp[5] = m12;
3574
3575         tmp[6] = m20;
3576         tmp[7] = m21;
3577         tmp[8] = m22;
3578
3579         Matrix3d.compute_svd( tmp, scales, rots);
3580
3581         return;
3582     }
3583
3584     /**
3585      * Creates a new object of the same class as this object.
3586      *
3587      * @return a clone of this instance.
3588      * @exception OutOfMemoryError if there is not enough memory.
3589      * @see java.lang.Cloneable
3590      * @since vecmath 1.3
3591      */
3592     public Object clone() {
3593         Matrix4d m1 = null;
3594         try {
3595             m1 = (Matrix4d)super.clone();
3596         } catch (CloneNotSupportedException e) {
3597             // this shouldn't happen, since we are Cloneable
3598             throw new InternalError();
3599         }
3600
3601         return m1;
3602     }
3603
3604     /**
3605          * Get the first matrix element in the first row.
3606          *
3607          * @return Returns the m00.
3608          *
3609          * @since vecmath 1.5
3610          */
3611         public final  double getM00() {
3612                 return m00;
3613         }
3614
3615         /**
3616          * Set the first matrix element in the first row.
3617          *
3618          * @param m00 The m00 to set.
3619          *
3620          *
3621          * @since vecmath 1.5
3622          */
3623         public final  void setM00(double m00) {
3624                 this.m00 = m00;
3625         }
3626
3627         /**
3628          * Get the second matrix element in the first row.
3629          *
3630          * @return Returns the m01.
3631          *
3632          * @since vecmath 1.5
3633          */
3634         public final  double getM01() {
3635                 return m01;
3636         }
3637
3638         /**
3639          * Set the second matrix element in the first row.
3640          *
3641          * @param m01 The m01 to set.
3642          *
3643          * @since vecmath 1.5
3644          */
3645         public  final void setM01(double m01) {
3646                 this.m01 = m01;
3647         }
3648
3649         /**
3650          * Get the third matrix element in the first row.
3651          *
3652          * @return Returns the m02.
3653          *
3654          * @since vecmath 1.5
3655          */
3656         public final double getM02() {
3657                 return m02;
3658         }
3659
3660         /**
3661          * Set the third matrix element in the first row.
3662          *
3663          * @param m02 The m02 to set.
3664          *
3665          * @since vecmath 1.5
3666          */
3667         public final  void setM02(double m02) {
3668                 this.m02 = m02;
3669         }
3670
3671         /**
3672          * Get first matrix element in the second row.
3673          *
3674          * @return Returns the m10.
3675          *
3676          * @since vecmath 1.5
3677          */
3678         public final  double getM10() {
3679                 return m10;
3680         }
3681
3682         /**
3683          * Set first matrix element in the second row.
3684          *
3685          * @param m10 The m10 to set.
3686          *
3687          * @since vecmath 1.5
3688          */
3689         public final  void setM10(double m10) {
3690                 this.m10 = m10;
3691         }
3692
3693         /**
3694          * Get second matrix element in the second row.
3695          *
3696          * @return Returns the m11.
3697          *
3698          * @since vecmath 1.5
3699          */
3700         public final  double getM11() {
3701                 return m11;
3702         }
3703
3704         /**
3705          * Set the second matrix element in the second row.
3706          *
3707          * @param m11 The m11 to set.
3708          *
3709          * @since vecmath 1.5
3710          */
3711         public final  void setM11(double m11) {
3712                 this.m11 = m11;
3713         }
3714
3715         /**
3716          * Get the third matrix element in the second row.
3717          *
3718          * @return Returns the m12.
3719          *
3720          * @since vecmath 1.5
3721          */
3722         public final  double getM12() {
3723                 return m12;
3724         }
3725
3726         /**
3727          * Set the third matrix element in the second row.
3728          *
3729          * @param m12 The m12 to set.
3730          *
3731          *
3732          * @since vecmath 1.5
3733          */
3734         public final  void setM12(double m12) {
3735                 this.m12 = m12;
3736         }
3737
3738         /**
3739          * Get the first matrix element in the third row.
3740          *
3741          * @return Returns the m20.
3742          *
3743          * @since vecmath 1.5
3744          */
3745         public final  double getM20() {
3746                 return m20;
3747         }
3748
3749         /**
3750          * Set the first matrix element in the third row.
3751          *
3752          * @param m20 The m20 to set.
3753          *
3754          * @since vecmath 1.5
3755          */
3756         public final void setM20(double m20) {
3757                 this.m20 = m20;
3758         }
3759
3760         /**
3761          * Get the second matrix element in the third row.
3762          *
3763          * @return Returns the m21.
3764          *
3765          * @since vecmath 1.5
3766          */
3767         public final double getM21() {
3768                 return m21;
3769         }
3770
3771         /**
3772          * Set the second matrix element in the third row.
3773          *
3774          * @param m21 The m21 to set.
3775          *
3776          * @since vecmath 1.5
3777          */
3778         public final void setM21(double m21) {
3779                 this.m21 = m21;
3780         }
3781
3782         /**
3783          * Get the third matrix element in the third row.
3784          *
3785          * @return Returns the m22.
3786          *
3787          * @since vecmath 1.5
3788          */
3789         public final double getM22() {
3790                 return m22;
3791         }
3792
3793         /**
3794          * Set the third matrix element in the third row.
3795          *
3796          * @param m22 The m22 to set.
3797          *
3798          * @since vecmath 1.5
3799          */
3800         public final void setM22(double m22) {
3801                 this.m22 = m22;
3802         }
3803
3804         /**
3805          * Get the fourth element of the first row.
3806          *
3807          * @return Returns the m03.
3808          *
3809          * @since vecmath 1.5
3810          */
3811         public final double getM03() {
3812                 return m03;
3813         }
3814
3815         /**
3816          * Set the fourth element of the first row.
3817          *
3818          * @param m03 The m03 to set.
3819          *
3820          * @since vecmath 1.5
3821          */
3822         public final void setM03(double m03) {
3823                 this.m03 = m03;
3824         }
3825
3826         /**
3827          * Get the fourth element of the second row.
3828          *
3829          * @return Returns the m13.
3830          *
3831          * @since vecmath 1.5
3832          */
3833         public final double getM13() {
3834                 return m13;
3835         }
3836
3837         /**
3838          * Set the fourth element of the second row.
3839          *
3840          * @param m13 The m13 to set.
3841          *
3842          * @since vecmath 1.5
3843          */
3844         public final void setM13(double m13) {
3845                 this.m13 = m13;
3846         }
3847
3848         /**
3849          * Get the fourth element of the third row.
3850          *
3851          * @return Returns the m23.
3852          *
3853          * @since vecmath 1.5
3854          */
3855         public final double getM23() {
3856                 return m23;
3857         }
3858
3859         /**
3860          * Set the fourth element of the third row.
3861          *
3862          * @param m23 The m23 to set.
3863          *
3864          * @since vecmath 1.5
3865          */
3866         public final void setM23(double m23) {
3867                 this.m23 = m23;
3868         }
3869
3870         /**
3871          * Get the first element of the fourth row.
3872          *
3873          * @return Returns the m30.
3874          *
3875          * @since vecmath 1.5
3876          */
3877         public final double getM30() {
3878                 return m30;
3879         }
3880
3881         /**
3882          * Set the first element of the fourth row.
3883          *
3884          * @param m30 The m30 to set.
3885          *
3886          * @since vecmath 1.5
3887          */
3888         public final void setM30(double m30) {
3889                 this.m30 = m30;
3890         }
3891
3892         /**
3893          * Get the second element of the fourth row.
3894          *
3895          * @return Returns the m31.
3896          *
3897          * @since vecmath 1.5
3898          */
3899         public final double getM31() {
3900                 return m31;
3901         }
3902
3903         /**
3904          * Set the second element of the fourth row.
3905          *
3906          * @param m31 The m31 to set.
3907          *
3908          * @since vecmath 1.5
3909          */
3910         public final void setM31(double m31) {
3911                 this.m31 = m31;
3912         }
3913
3914         /**
3915          * Get the third element of the fourth row.
3916          *
3917          * @return Returns the m32.
3918          *
3919          *
3920          * @since vecmath 1.5
3921          */
3922         public final double getM32() {
3923                 return m32;
3924         }
3925
3926         /**
3927          * Set the third element of the fourth row.
3928          *
3929          * @param m32 The m32 to set.
3930          *
3931          * @since vecmath 1.5
3932          */
3933         public final void setM32(double m32) {
3934                 this.m32 = m32;
3935         }
3936
3937         /**
3938          * Get the fourth element of the fourth row.
3939          *
3940          * @return Returns the m33.
3941          *
3942          * @since vecmath 1.5
3943          */
3944         public final double getM33() {
3945                 return m33;
3946         }
3947
3948         /**
3949          * Set the fourth element of the fourth row.
3950          *
3951          * @param m33 The m33 to set.
3952          *
3953          * @since vecmath 1.5
3954          */
3955         public final void setM33(double m33) {
3956                 this.m33 = m33;
3957         }
3958 }