Index: realspace_modeling/geoshapespy/libgeoshapespy/ellipsoid.cc
===================================================================
--- realspace_modeling/geoshapespy/libgeoshapespy/ellipsoid.cc	(revision a778da591a2680b9cec2b94403f040cefd931b7f)
+++ realspace_modeling/geoshapespy/libgeoshapespy/ellipsoid.cc	(revision a778da591a2680b9cec2b94403f040cefd931b7f)
@@ -0,0 +1,165 @@
+/** \file Ellipsoid.cc */
+#include <cmath>
+#include <cassert>
+#include <algorithm>
+//MS library does not define min max in algorithm
+//#include "minmax.h"
+#include "ellipsoid.h"
+#include "myutil.h"
+
+using namespace std;
+
+Ellipsoid::Ellipsoid()
+{
+  rx_ = 0;
+  ry_ = 0;
+  rz_ = 0;
+}
+
+Ellipsoid::Ellipsoid(double rx, double ry, double rz)
+{
+  rx_ = rx;
+  ry_ = ry;
+  rz_ = rz;
+  xedge_plus_ = Point3D(rx,0,0);
+  xedge_minus_ = Point3D(-rx,0,0);
+  yedge_plus_ = Point3D(0,ry,0);
+  yedge_minus_ = Point3D(0,-ry,0);
+  zedge_plus_ = Point3D(0,0,rz);
+  zedge_minus_ = Point3D(0,0,-rz);
+}
+
+Ellipsoid::~Ellipsoid()
+{
+}
+
+void Ellipsoid::SetRadii(double rx, double ry, double rz)
+{
+  rx_ = rx;
+  ry_ = ry;
+  rz_ = rz;
+}
+
+double Ellipsoid::GetRadiusX()
+{
+  return rx_;
+}
+
+double Ellipsoid::GetRadiusY()
+{
+  return ry_;
+}
+
+double Ellipsoid::GetRadiusZ()
+{
+  return rz_;
+}
+
+double Ellipsoid::GetMaxRadius()
+{
+  double maxr = max(max(rx_,ry_),max(ry_,rz_));
+  return maxr;
+}
+
+ShapeType Ellipsoid::GetShapeType() const
+{
+  return ELLIPSOID;
+}
+
+double Ellipsoid::GetVolume()
+{
+  double V = (4./3.)*pi*rx_*ry_*rz_;
+  return V;
+}
+
+void Ellipsoid::GetFormFactor(IQ * iq)
+{
+  /** number of I for output, equal to the number of rows of array IQ*/
+  /** to be finished */
+}
+
+Point3D Ellipsoid::GetAPoint(double sld)
+{
+  static int max_try = 100;
+  for (int i = 0; i < max_try; ++i) {
+    double x = (ran1()-0.5) * 2 * rx_;
+    double y = (ran1()-0.5) * 2 * ry_;	
+    double z = (ran1()-0.5) * 2 * rz_;
+
+    if ((square(x/rx_) + square(y/ry_) + square(z/rz_)) <= 1){
+      Point3D apoint(x,y,z,sld);
+      return apoint;
+    }
+  }
+
+  std::cerr << "Max try "
+            << max_try
+            << " is reached while generating a point in ellipsoid" << std::endl;
+  return Point3D(0, 0, 0);
+}
+
+bool Ellipsoid::IsInside(const Point3D& point) const
+{
+  //x, y, z axis are internal axis
+  bool isOutsideX = false;
+  Point3D pointOnX = point.getInterPoint(GetXaxisPlusEdge(),GetXaxisMinusEdge(),&isOutsideX);
+  bool isOutsideY = false;
+  Point3D pointOnY = point.getInterPoint(GetYaxisPlusEdge(),GetYaxisMinusEdge(),&isOutsideY);
+  bool isOutsideZ = false;
+  Point3D pointOnZ = point.getInterPoint(GetZaxisPlusEdge(),GetZaxisMinusEdge(),&isOutsideZ);
+
+  if (isOutsideX || isOutsideY || isOutsideZ){
+    //one is outside axis is true -> the point is not inside
+    return false;    
+  }
+  else{
+    Point3D pcenter = GetCenterP();
+    double distX = pointOnX.distanceToPoint(pcenter);
+    double distY = pointOnY.distanceToPoint(pcenter);
+    double distZ = pointOnZ.distanceToPoint(pcenter);
+    return ((square(distX/rx_)+square(distY/ry_)+square(distZ/rz_)) <= 1);
+  }
+
+}
+
+Point3D Ellipsoid::GetXaxisPlusEdge() const
+{
+  Point3D new_edge(xedge_plus_);
+  new_edge.Transform(GetOrientation(),GetCenter());
+  return new_edge;
+}
+
+Point3D Ellipsoid::GetXaxisMinusEdge() const
+{
+  Point3D new_edge(xedge_minus_);
+  new_edge.Transform(GetOrientation(),GetCenter());
+  return new_edge;
+}
+
+Point3D Ellipsoid::GetYaxisPlusEdge() const
+{
+  Point3D new_edge(yedge_plus_);
+  new_edge.Transform(GetOrientation(),GetCenter());
+  return new_edge;
+}
+
+Point3D Ellipsoid::GetYaxisMinusEdge() const
+{
+  Point3D new_edge(yedge_minus_);
+  new_edge.Transform(GetOrientation(),GetCenter());
+  return new_edge;
+}
+
+Point3D Ellipsoid::GetZaxisPlusEdge() const
+{
+  Point3D new_edge(zedge_plus_);
+  new_edge.Transform(GetOrientation(),GetCenter());
+  return new_edge;
+}
+
+Point3D Ellipsoid::GetZaxisMinusEdge() const
+{
+  Point3D new_edge(zedge_minus_);
+  new_edge.Transform(GetOrientation(),GetCenter());
+  return new_edge;
+}
