--- trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp	2006/02/22 20:35:16	891
+++ trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp	2006/02/23 23:16:43	892
@@ -43,6 +43,7 @@
 #include "math/LU.hpp"
 #include "math/DynamicRectMatrix.hpp"
 #include "math/SquareMatrix3.hpp"
+#include "utils/OOPSEConstant.hpp"
 namespace oopse {
 /**
  * Reference:
@@ -50,7 +51,28 @@ namespace oopse {
  * Comparison of Different Modeling and Computational Procedures. 
  * Biophysical Journal, 75(6), 3044, 1999
  */
-bool HydrodynamicsModel::calcHydrodyanmicsProps(double eta) {
+
+HydrodynamicsModel::HydrodynamicsModel(StuntDouble* sd, const DynamicProperty& extraParams) : sd_(sd){
+    DynamicProperty::const_iterator iter;
+
+    iter = extraParams.find("Viscosity");
+    if (iter != extraParams.end()) {
+        boost::any param = iter->second;
+        viscosity_ = boost::any_cast<double>(param);
+    }else {
+        std::cout << "HydrodynamicsModel Error\n" ;
+    }
+
+    iter = extraParams.find("Temperature");
+    if (iter != extraParams.end()) {
+        boost::any param = iter->second;
+        temperature_ = boost::any_cast<double>(param);
+    }else {
+        std::cout << "HydrodynamicsModel Error\n" ;
+    }    
+}
+ 
+bool HydrodynamicsModel::calcHydrodyanmicsProps() {
     if (!createBeads(beads_)) {
         std::cout << "can not create beads" << std::endl;
         return false;
@@ -60,6 +82,10 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps(double
     DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
     DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
     Mat3x3d I;
+    I(0, 0) = 1.0;
+    I(1, 1) = 1.0;
+    I(2, 2) = 1.0;
+    
     for (std::size_t i = 0; i < nbeads; ++i) {
         for (std::size_t j = 0; j < nbeads; ++j) {
             Mat3x3d Tij;
@@ -70,10 +96,10 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps(double
                 double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
                 Mat3x3d tmpMat;
                 tmpMat = outProduct(beads_[i].pos, beads_[j].pos) / rij2;
-                double constant = 8.0 * NumericConstant::PI * eta * rij;
+                double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
                 Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
             }else {
-                double constant = 1.0 / (6.0 * NumericConstant::PI * eta * beads_[i].radius);
+                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
                 Tij(0, 0) = constant;
                 Tij(1, 1) = constant;
                 Tij(2, 2) = constant;
@@ -84,7 +110,6 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps(double
 
     //invert B Matrix
     invertMatrix(B, C);
-    
     //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
     std::vector<Mat3x3d> U;
     for (int i = 0; i < nbeads; ++i) {
@@ -115,35 +140,70 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps(double
     Xi.setSubMatrix(0, 0, Xitt);
     Xi.setSubMatrix(0, 3, Xitr.transpose());
     Xi.setSubMatrix(3, 0, Xitr);
-    Xi.setSubMatrix(3, 3, Xitt);
-    invertMatrix(Xi, Do);
+    Xi.setSubMatrix(3, 3, Xirr);
+    //invertMatrix(Xi, Do);
+    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
+    //Do *= kt;    
 
+
     Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
     Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
     Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
-    Do.getSubMatrix(0, 0 , Dott);
-    Do.getSubMatrix(3, 0, Dotr);
-    Do.getSubMatrix(3, 3, Dorr);
 
+    Mat3x3d XirrInv(0.0);
+    Mat3x3d XirrCopy;
+    XirrCopy = Xirr;
+    
+    Mat3x3d XittInv(0.0);
+    Mat3x3d XittCopy;
+    XittCopy = Xitt;
+    invertMatrix(XittCopy, XittInv);
+
+    Mat3x3d tmp;
+    Mat3x3d tmpInv;
+    tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
+
+    const static Mat3x3d zeroMat(0.0);
+    if (!invertMatrix(tmp, tmpInv)) {
+        tmpInv = zeroMat;
+    }
+
+    Dott = kt * tmpInv;
+    Dotr = -kt*XirrInv * Xitr * tmpInv;
+
+    tmp = Xirr - Xitr * XittInv * Xitr.transpose();
+    
+    if(!invertMatrix(tmp, tmpInv)) {
+        tmpInv = zeroMat;
+    }
+    Dorr = kt * tmpInv;
+
+    //Do.getSubMatrix(0, 0 , Dott);
+    //Do.getSubMatrix(3, 0, Dotr);
+    //Do.getSubMatrix(3, 3, Dorr);
+
     //calculate center of diffusion
-    Mat3x3d tmpMat;
-    tmpMat(0, 0) = Dorr(1, 1) + Dorr(2, 2);
-    tmpMat(0, 1) = - Dorr(0, 1);
-    tmpMat(0, 2) = -Dorr(0, 2);
-    tmpMat(1, 0) = -Dorr(0, 1);
-    tmpMat(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
-    tmpMat(1, 2) = -Dorr(1, 2);
-    tmpMat(2, 0) = -Dorr(0, 2);
-    tmpMat(2, 1) = -Dorr(1, 2);
-    tmpMat(2, 2) = Dorr(1, 1) + Dorr(0, 0);
+    tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
+    tmp(0, 1) = - Dorr(0, 1);
+    tmp(0, 2) = -Dorr(0, 2);
+    tmp(1, 0) = -Dorr(0, 1);
+    tmp(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
+    tmp(1, 2) = -Dorr(1, 2);
+    tmp(2, 0) = -Dorr(0, 2);
+    tmp(2, 1) = -Dorr(1, 2);
+    tmp(2, 2) = Dorr(1, 1) + Dorr(0, 0);
 
     Vector3d tmpVec;
     tmpVec[0] = Dotr(1, 2) - Dotr(2, 1);
     tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
     tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
-        
-    Vector3d rod = tmpMat.inverse() * tmpVec;
 
+    if(!invertMatrix(tmp, tmpInv)) {
+        tmpInv = zeroMat;
+    }
+    
+    Vector3d rod = tmpInv * tmpVec;
+
     //calculate Diffusion Tensor at center of diffusion
     Mat3x3d Uod;
     Uod.setupSkewMat(rod);
@@ -165,11 +225,49 @@ void HydrodynamicsModel::writeBeads(std::ostream& os) 
 }
 
 void HydrodynamicsModel::writeBeads(std::ostream& os) {
+    std::vector<BeadParam>::iterator iter;
+    os << beads_.size() << std::endl;
+    os << "Generated by Hydro" << std::endl;
+    for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
+        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
+    }
 
 }
 
 void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
+    os << "//viscosity = " << viscosity_ << std::endl;
+    os << "//temperature = " << temperature_<< std::endl;
+    std::vector<BeadParam>::iterator iter;
+    os << sd_->getType() << "\n";
 
+    os << "//diffusion center" << std::endl;
+    os << props_.diffCenter << std::endl;
+
+    os << "//translational diffusion tensor" << std::endl;
+    os << props_.transDiff << std::endl;
+
+    os << "//translational diffusion tensor" << std::endl;
+    os << props_.transRotDiff << std::endl;
+
+    os << "//rotational diffusion tensor" << std::endl;
+    os << props_.rotDiff << std::endl;
+    
+    /*
+    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
+
+    os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t"
+        << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t"
+        << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n";
+    
+    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
+        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
+        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
+
+    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
+        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
+        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
+        << std::endl;
+    */
 }
 
 }