--- trunk/src/applications/hydrodynamics/AnalyticalModel.cpp	2006/05/16 02:06:37	956
+++ trunk/src/applications/hydrodynamics/AnalyticalModel.cpp	2006/05/24 18:31:12	972
@@ -38,101 +38,45 @@
  * University of Notre Dame has been advised of the possibility of
  * such damages.
  */
-#include "hydrodynamics/AnalyticalModel.hpp"
+#include "applications/hydrodynamics/AnalyticalModel.hpp"
 #include "hydrodynamics/Sphere.hpp"
 #include "hydrodynamics/Ellipsoid.hpp"
-#include "hydrodynamics/CompositeShape.hpp"
+#include "applications/hydrodynamics/CompositeShape.hpp"
 #include "math/LU.hpp"
 namespace oopse {
   
-  bool AnalyticalModel::calcHydroProps(Sphere* sphere, double viscosity, double temperature) {
+  bool AnalyticalModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
     
-    double radius = sphere->getRadius(); 
     HydroProps props;
-    props.center =V3Zero;
-    double Xitt  = 6.0 * NumericConstant::PI * viscosity * radius;
-    double Xirr = 8.0 * NumericConstant::PI * viscosity * radius * radius * radius;
-    props.Xi(0, 0) = Xitt;
-    props.Xi(1, 1) = Xitt;
-    props.Xi(2, 2) = Xitt;
-    props.Xi(3, 3) = Xirr;
-    props.Xi(4, 4) = Xirr;
-    props.Xi(5, 5) = Xirr;
-    
-    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
-    props.Xi *= convertConstant;
-    Mat6x6d XiCopy = props.Xi;
-    invertMatrix(XiCopy, props.D);
-    double kt = OOPSEConstant::kB * temperature;
-    props.D *= kt;
-    props.Xi *= OOPSEConstant::kb * temperature;
-    
-    setCR(props);
-    setCD(props);
-    
-    return true;    
-}
-
-  /**
-   * Reference:
-   * (2) F. Perrin , J. Phys. Radium, [7] 5, 497-511, 1934
-   * (3) F. Perrin, J. Phys. Radium, [7] 7, 1-11, 1936
-   */        
-  bool AnalyticalModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) {
-    
-    double rMajor = ellipsoid->getRMajor();
-    double rMinor = ellipsoid->getRMinor();
-    
-    double a = rMinor;
-    double b = rMajor;
-    double a2 = a * a;
-    double b2 = b* b;
-    
-    double p = a /b;
-    double S;
-    if (p > 1.0) { //prolate
-      S = 2.0/sqrt(a2 - b2) * log((a + sqrt(a2-b2))/b);
-    } else { //oblate
-      S = 2.0/sqrt(b2 - a2) * atan(sqrt(b2-a2)/a);
-    }
-    
-    double P = 1.0/(a2 - b2) * (S - 2.0/a);
-    double Q = 0.5/(a2-b2) * (2.0*a/b2 - S);
-    
-    double transMinor = 16.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-b2)*S -2.0*a);
-    double transMajor = 32.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-3.0*b2)*S +2.0*a);
-    double rotMinor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2 - b2) * b2 /(2.0*a -b2*S);
-    double rotMajor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2*a2 - b2*b2)/((2.0*a2-b2)*S-2.0*a);
-    
-    
-    HydroProps props;
-    
-    props.Xi(0,0) = transMajor;
-    props.Xi(1,1) = transMajor;
-    props.Xi(2,2) = transMinor;
-    props.Xi(3,3) = rotMajor;
-    props.Xi(4,4) = rotMajor;
-    props.Xi(5,5) = rotMinor;
-    
-    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
-    props.Xi *= convertConstant;    
-    
-    Mat6x6d XiCopy = props.Xi;
-    invertMatrix(XiCopy, props.D);
-    double kt = OOPSEConstant::kB * temperature;
-    props.D *= kt;
-    props.Xi *= OOPSEConstant::kb * temperature;
-    
-    setCR(props);
-    setCD(props);
-    
-    return true;
+    Sphere* sphere = dynamic_cast<Sphere*>(shape);        
+    if (sphere != NULL) {
+      props = sphere->getHydroProps(viscosity, temperature);
+      setCR(props);
+      setCD(props);
+      return true;
+    } else {
+      Ellipsoid* ellipsoid = dynamic_cast<Ellipsoid*>(shape);        
+      if (ellipsoid != NULL) {
+        props = ellipsoid->getHydroProps(viscosity, temperature); 
+        setCR(props);
+        setCD(props);
+        return true;
+      } else {
+        CompositeShape* composite = dynamic_cast<CompositeShape*>(shape);
+        if (composite != NULL) {
+          return false;
+        } else {
+          sprintf( painCave.errMsg,
+                   "Could not figure out what kind of shape this is!\n");
+          painCave.severity = OOPSE_ERROR;
+          painCave.isFatal = 1;
+          simError();    
+          return false;
+        }
+      }
+    }          
   }
-  
-  bool AnalyticalModel::calcHydroProps(CompositeShape* compositeShape, double viscosity, double temperature) {
-    return false;
-  }
-  
+    
   void AnalyticalModel::writeBeads(std::ostream& os) {
     os << "1\n";
     os << "Generated by Hydro\n";