Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/Make.mm
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/Make.mm	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/Make.mm	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,69 @@
+# -*- Makefile -*-
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+
+include local.def
+
+PROJECT = pointsmodelpy
+PACKAGE = libpointsmodelpy
+
+PROJ_SAR = $(BLD_LIBDIR)/$(PACKAGE).$(EXT_SAR)
+PROJ_DLL = $(BLD_BINDIR)/$(PACKAGE).$(EXT_SO)
+PROJ_TMPDIR = $(BLD_TMPDIR)/$(PROJECT)/$(PACKAGE)
+PROJ_CLEAN += $(PROJ_SAR) $(PROJ_DLL)
+
+PROJ_SRCS = \
+	lores_model.cc	\
+	points_model.cc	\
+        pdb_model.cc  \
+        complex_model.cc
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# build the library
+
+all: $(PROJ_SAR) export
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ifeq (Win32, ${findstring Win32, $(PLATFORM_ID)})
+
+# build the shared object
+$(PROJ_SAR): product_dirs $(PROJ_OBJS)
+	$(CXX) $(LCXXFLAGS) -o $(PROJ_DLL) \
+	-Wl,--out-implib=$(PROJ_SAR) $(PROJ_OBJS)
+
+# export
+export:: export-headers export-libraries export-binaries
+
+else
+
+# build the shared object
+$(PROJ_SAR): product_dirs $(PROJ_OBJS)
+	$(CXX) $(LCXXFLAGS) -o $(PROJ_SAR) $(PROJ_OBJS) -L$(BLD_LIBDIR)
+
+# export
+export:: export-headers export-libraries
+
+endif
+
+EXPORT_HEADERS = \
+	lores_model.h	\
+	points_model.h	\
+        complex_model.h
+
+EXPORT_LIBS = $(PROJ_SAR)
+EXPORT_BINS = $(PROJ_DLL)
+
+
+# version
+# $Id$
+
+#
+# End of file
Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/complex_model.cc
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/complex_model.cc	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/complex_model.cc	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,97 @@
+/** \file complex_model.cc */
+
+#include "complex_model.h"
+#include <algorithm>
+#include <stdexcept>
+
+using namespace std;
+
+ComplexModel::ComplexModel(){
+
+}
+
+void ComplexModel::Add(PointsModel *pm)
+{
+  models_.push_back(pm);
+}
+
+double ComplexModel::GetDimBound()
+{
+  //Get the vector of centers of pointsmodel instances
+  //and a vector of individual boundary
+  vector<PointsModel *>::iterator itr;
+  vector<Point3D> veccenter;
+  vector<double> bounds;
+
+  for(itr = models_.begin(); itr != models_.end(); ++itr){
+    if ((*itr)->GetDimBound() != 0){
+      Point3D pp((*itr)->GetCenter()[0],(*itr)->GetCenter()[1],(*itr)->GetCenter()[2]);
+      veccenter.push_back(pp);
+      bounds.push_back((*itr)->GetDimBound());
+    }
+  }
+
+  //max bound
+  double maxbound = *max_element(bounds.begin(),bounds.end());
+
+  //max distance
+  vector<double> vecdist;
+  size_t num = veccenter.size();
+
+  if (num > 1){
+    for (size_t i = 0; i != num; ++i){
+      for (size_t j= 1; j != num; ++j){
+        double dist = veccenter[i].distanceToPoint(veccenter[j]);
+        vecdist.push_back(dist);
+      }
+    }
+  }
+  else{
+    vecdist.push_back(maxbound);
+  }
+
+  double maxcenterdist = *max_element(vecdist.begin(),vecdist.end());
+
+  double finalbound = maxbound + maxcenterdist;
+
+  return finalbound;
+}
+
+vector<double> ComplexModel::GetCenter()
+{
+  double sumx = 0, sumy = 0, sumz = 0;
+  size_t num = 0;
+
+  vector<PointsModel *>::iterator itr;
+  for(itr = models_.begin(); itr != models_.end(); ++itr){
+    sumx += (*itr)->GetCenter()[0];
+    sumy += (*itr)->GetCenter()[1];
+    sumz += (*itr)->GetCenter()[2];
+    ++num;
+  }
+
+  vector<double> v(3);
+  v[0] = sumx/num;
+  v[1] = sumy/num;
+  v[2] = sumz/num;
+  center_ = v;
+  return center_;
+}
+
+int ComplexModel::GetPoints(Point3DVector &vp)
+{
+  if (vp.size() != 0){
+    throw runtime_error("GetPoints(Point3DVector &VP):VP has to be empty"); 
+  }
+
+  vector<PointsModel *>::iterator itr;
+
+  for(itr = models_.begin(); itr != models_.end(); ++itr){
+    vector<Point3D> temp;
+    (*itr)->GetPoints(temp);    	
+    if (temp.size() != 0){
+      vp.insert(vp.end(),temp.begin(),temp.end());
+    }
+  }
+  return vp.size();
+}
Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/complex_model.h
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/complex_model.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/complex_model.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,36 @@
+/** \file complex_model.h child class of PointsModel */
+
+#ifndef COMPLEX_MODEL_H
+#define COMPLEX_MODEL_H
+
+#include "points_model.h"
+#include <string>
+#include <map>
+
+/**
+ *  Class ComplexModel : container class for LORESModel & PDBModel
+ *  The main functionality is to merge points from instances of 
+ *  LORESModel & PDBModel
+ */
+
+class ComplexModel : public PointsModel {
+ public:
+  ComplexModel();
+  
+  //add PointsModel instance
+  void Add(PointsModel *);
+
+  //Parse all coordinates from ATOM section 
+  //of the PDB file into vector of points
+  int GetPoints(Point3DVector &);
+
+  //Get distance boundary for the pointsmodel instances
+  double GetDimBound();
+
+  vector<double> GetCenter();
+
+ private:
+  vector<PointsModel *> models_;
+};
+
+#endif
Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/local.def
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/local.def	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/local.def	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,32 @@
+# -*- Makefile -*-
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+#
+# <LicenseText>
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+
+#
+# Local (project) definitions
+#
+
+# C++
+
+PROJ_CXX_INCLUDES += ../../iqPy/libiqPy   \
+                     ../../iqPy/libiqPy/tnt \
+		     ../../geoshapespy/libgeoshapespy   \
+		     ../../analmodelpy/libanalmodelpy   \
+		     ./ 
+
+    PROJ_CXX_FLAGS += $(CXX_SOFLAGS)
+    PROJ_LCXX_FLAGS += $(LCXX_SARFLAGS) -liqPy -lgeoshapespy
+
+# version
+# $Id$
+
+# End of file
Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/lores_model.cc
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/lores_model.cc	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/lores_model.cc	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,237 @@
+/** \file lores_model.cc */
+
+#include <algorithm>
+#include <stdexcept>
+#include "lores_model.h"
+#include "sphere.h"
+#include "hollow_sphere.h"
+#include "cylinder.h"
+#include "ellipsoid.h"
+#include "single_helix.h"
+#include "Point3D.h"
+
+LORESModel::LORESModel(double density)
+{
+  density_ = density;
+}
+
+GeoShape* LORESModel::GetGeoShape(GeoShape& geo_shape)
+{
+  GeoShape* shape = NULL;
+
+  switch (geo_shape.GetShapeType()){
+    case SPHERE:
+      shape = new Sphere(static_cast<const Sphere&>(geo_shape));
+      break;
+    case HOLLOWSPHERE:
+      shape = new HollowSphere(static_cast<const HollowSphere&>(geo_shape));
+      break;
+    case CYLINDER:
+      shape = new Cylinder(static_cast<const Cylinder&>(geo_shape));
+      break;
+    case ELLIPSOID:
+      shape = new Ellipsoid(static_cast<const Ellipsoid&>(geo_shape));
+      break;
+    case SINGLEHELIX:
+      shape = new SingleHelix(static_cast<const SingleHelix&>(geo_shape));
+      break;
+  
+  }
+
+  return shape;
+}
+
+LORESModel::~LORESModel()
+{
+  for (RealSpaceShapeCollection::iterator it = shapes_.begin();
+       it != shapes_.end(); ++it) {
+    delete *it;
+  }
+}
+
+void LORESModel::SetDensity(double density)
+{
+  density_ = density;
+}
+
+double LORESModel::GetDensity()
+{
+  return density_;
+}
+
+// Add a shape into LORES Model
+void LORESModel::Add(GeoShape& geo_shape,
+                     double sld) {
+  GeoShape* shape = GetGeoShape(geo_shape);
+  assert(shape != NULL);
+
+  RealSpaceShape* real_shape = new RealSpaceShape(shape);
+  FillPoints(real_shape, sld);
+  shapes_.push_back(real_shape);
+}
+
+// Delete ith shape in shapes_ list
+// If we have 3 shapes in our model, the index starts
+// from 0, which means we need to call Delete(0) to delete
+// the first shape, and call Delete(1) to delete the second
+// shape, etc..
+void LORESModel::Delete(size_t i) {
+  if (i >= shapes_.size()) {
+    std::cerr << "Delete shapes out of scope" << std::endl;
+    return;
+  }
+
+  RealSpaceShape* real_shape = shapes_[i];
+  if (i + 1 < shapes_.size()) {
+    // if it is not the last shape, we have to distribute its points
+    // to the shapes after i-th shape.
+    const Point3DVector& points = real_shape->points;
+    for (Point3DVector::const_iterator pit = points.begin();
+         pit != points.end(); ++pit)
+      DistributePoint(*pit, i + 1);
+  }
+
+  shapes_.erase(shapes_.begin() + i);
+  delete real_shape;
+}
+
+// Get the points in the realspaceshapecollection
+int LORESModel::GetPoints(Point3DVector &pp)
+{
+  if (pp.size() != 0){
+    throw runtime_error("GetPoints(Point3DVector &VP):VP has to be empty"); 
+  }
+
+  if (shapes_.size() != 0){
+    for (size_t j = 0; j <shapes_.size(); ++j ){
+      pp.insert(pp.end(),shapes_[j]->points.begin(),shapes_[j]->points.end());
+    }
+  }
+  return pp.size();
+}
+
+//Write points to a file, mainly for testing right now
+void LORESModel::WritePoints2File(Point3DVector &vp){
+  ofstream outfile("test.coor");
+  for(size_t i=0; i<vp.size(); ++i){
+    outfile<<vp[i].getX()<<" "<<vp[i].getY()<<" "<<vp[i].getZ()<<" "<<vp[i].getSLD()<<endl;
+  }
+}
+
+double LORESModel::GetDimBound()
+{
+  if(shapes_.size() == 0)
+    return 0;
+
+  //find maximum distance between centers of shapes
+  //get the vector of centers
+  vector<Point3D> vec_center;
+  for (size_t m = 0; m < shapes_.size(); ++m){
+    assert(shapes_[m]->shape != NULL);
+    vector<double> center(3);
+    center = shapes_[m]->shape->GetCenter();
+    Point3D p_center;
+    p_center.set(center[0],center[1],center[2]);
+    vec_center.push_back(p_center);
+  }
+  size_t vecsize = vec_center.size();
+
+  //get the maximum distance among centers
+  double max_cen_dist;
+  if (vecsize == 1){
+    max_cen_dist = 0;
+  }
+  else{
+    vector<double> vecdist;
+    for (size_t m1=0; m1<vecsize -1; ++m1){
+      for (size_t m2 = m1+1; m2<vecsize; ++m2){
+	double dist = vec_center[m1].distanceToPoint(vec_center[m2]);
+	vecdist.push_back(dist);
+      }
+    }
+    max_cen_dist = *max_element(vecdist.begin(), vecdist.end());
+  }
+  //find the maximum radius for individual shape
+  vector<double> maxradii;
+  for (size_t n = 0; n < shapes_.size(); ++n){
+    assert(shapes_[n]->shape != NULL);
+    double maxradius = shapes_[n]->shape->GetMaxRadius();
+    maxradii.push_back(maxradius);
+  }
+  double max_maxradius = *max_element(maxradii.begin(),maxradii.end());
+  
+  return 2*(max_cen_dist/2 + max_maxradius);
+}
+
+
+// Distribute points of the shape we are going to delete
+void LORESModel::DistributePoint(const Point3D& point, size_t i)
+{
+  for (size_t k = i; k < shapes_.size(); ++k) {
+    assert(shapes_[k]->shape != NULL);
+    if (shapes_[k]->shape->IsInside(point)) {
+      shapes_[k]->points.push_back(point);
+      return;
+    }
+  }
+}
+
+void LORESModel::FillPoints(RealSpaceShape* real_shape, double sld)
+{
+  assert(real_shape != NULL);
+
+  GeoShape* shape = real_shape->shape;
+  assert(shape != NULL);
+
+  int npoints = static_cast<int>(density_ * shape->GetVolume());
+
+  for (int i = 0; i < npoints; ++i){
+    Point3D apoint = shape->GetAPoint(sld);
+    apoint.Transform(shape->GetOrientation(),shape->GetCenter());
+    if (!IsInside(apoint)){
+      real_shape->points.push_back(apoint);
+    }
+  }
+
+}
+
+bool LORESModel::IsInside(const Point3D& point)
+{
+  for (RealSpaceShapeCollection::const_iterator it = shapes_.begin();
+       it != shapes_.end(); ++it) {
+    const GeoShape* shape = (*it)->shape;
+    assert(shape != NULL);
+
+    if (shape->IsInside(point)) {
+       return true;
+    }
+  }
+
+  return false;
+}
+
+vector<double> LORESModel::GetCenter()
+{
+  //get the vector of centers from the list of shapes
+  size_t numshapes = 0;
+  double sumx = 0, sumy = 0, sumz = 0;
+  for (size_t m = 0; m < shapes_.size(); ++m){
+    assert(shapes_[m]->shape != NULL);
+    vector<double> center(3);
+    center = shapes_[m]->shape->GetCenter();
+
+    sumx += center[0];
+    sumy += center[1];
+    sumz += center[2];
+
+    ++numshapes;
+  }
+
+  vector<double> shapescenter(3);
+  shapescenter[0]= sumx/numshapes;
+  shapescenter[1]= sumy/numshapes;
+  shapescenter[2]= sumz/numshapes;
+
+  center_ = shapescenter;
+  return center_;
+}
Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/lores_model.h
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/lores_model.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/lores_model.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,73 @@
+/** \file lores_model.h child class of PointsModel*/
+
+#ifndef LORESMODEL_H
+#define LORESMODEL_H
+
+#include <deque>
+#include "points_model.h"
+#include "geo_shape.h"
+//temporary
+#include <fstream>
+
+/**
+ *  Class LORESModel, low resolution models
+ */
+
+class LORESModel : public PointsModel {
+ private:
+  typedef vector<Point3D> PointsVector;
+
+  struct RealSpaceShape {
+    GeoShape* shape;
+    PointsVector points;
+
+    RealSpaceShape(GeoShape* s = NULL) : shape(s) {
+    }
+
+    ~RealSpaceShape() {
+      if (shape) delete shape;
+    }
+  };
+
+  typedef deque<RealSpaceShape*> RealSpaceShapeCollection;
+
+ public:
+  LORESModel(double density = 1.0);
+  ~LORESModel();
+
+  // Change density
+  void SetDensity(double density);
+  double GetDensity();
+
+  // Add new shape
+  void Add(GeoShape& geo_shape, double sld = 1);
+
+  // Delete ith shape at shapes_
+  void Delete(size_t i);
+
+  int GetPoints(Point3DVector &);
+  //Write points to a file, mainly for testing right now
+  void WritePoints2File(Point3DVector &);
+
+  //get the maximum possible dimension
+  double GetDimBound();
+
+  //will be used in determining the maximum distance for
+  //P(r) calculation for a complex model (merge several 
+  //pointsmodel instance together
+  vector<double> GetCenter();
+
+ protected:
+  GeoShape* GetGeoShape(GeoShape& geo_shape);
+  void FillPoints(RealSpaceShape* real_shape, double sld);
+  bool IsInside(const Point3D& point);
+  void DistributePoint(const Point3D& point, size_t i);
+
+ private:
+  RealSpaceShapeCollection shapes_;
+  int npoints_;
+  double density_;
+  
+};
+
+#endif
Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/pdb_model.cc
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/pdb_model.cc	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/pdb_model.cc	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,123 @@
+/** \file pdb_model.cc */
+
+#include "pdb_model.h"
+#include <fstream>
+#include <string>
+#include <cstdlib>
+#include <cmath>
+#include <stdexcept>
+
+using namespace std;
+
+typedef map<string,double>::const_iterator mymapitr;
+
+PDBModel::PDBModel(){
+  res_sld_["ALA"] = 1.645;  res_sld_["ARG"] = 3.466;
+  res_sld_["ASN"] = 3.456;  res_sld_["ASP"] = 3.845;
+  res_sld_["CYS"] = 1.930;  res_sld_["GLU"] = 3.762;
+  res_sld_["GLN"] = 3.373;  res_sld_["GLY"] = 1.728;
+  res_sld_["HIS"] = 4.959;  res_sld_["LLE"] = 1.396;
+  res_sld_["LEU"] = 1.396;  res_sld_["LYS"] = 1.586;
+  res_sld_["MET"] = 1.763;  res_sld_["PHE"] = 4.139;
+  res_sld_["PRO"] = 2.227;  res_sld_["SER"] = 2.225;
+  res_sld_["THR"] = 2.142;  res_sld_["TRP"] = 6.035;
+  res_sld_["TYR"] = 4.719;  res_sld_["VAL"] = 1.479;
+
+  x_max_ = -10000; x_min_ = 10000;
+  y_max_ = -10000; y_min_ = 10000;
+  z_max_ = -10000; z_min_ = 10000;
+}
+
+void PDBModel::AddPDB(const string &pdbfile)
+{
+  pdbnames_.push_back(pdbfile);
+}
+
+void PDBModel::AddPDB(const char* pdbfile)
+{
+  string sname(pdbfile);
+  pdbnames_.push_back(sname);
+}
+
+int PDBModel::GetPoints(Point3DVector &vp)
+{
+  if (vp.size() != 0){
+    throw runtime_error("PDBModel::GetPoints(Point3DVector &VP):VP has to be empty"); 
+  }
+  vector<string>::iterator itr;
+  for (itr = pdbnames_.begin(); itr!=pdbnames_.end(); ++itr){
+    ifstream infile(itr->c_str());
+    if (!infile){
+      cerr << "error: unable to open input file: " 
+	   << infile << endl;
+    }
+    string line;
+    while (getline(infile,line)){
+      size_t len = line.size();
+      string header;
+      
+      //for a line with size >=4, retrieve the header
+      //if the header == "ATOM""
+      //make sure the length of the line >=54, then
+      //retrieve the corresponding coordinates,convert str to double
+      //then assign them to a point3d,which will be append to vp with a SLD
+      if (len >= 4){
+	for (size_t i=0; i !=4; ++i) header += line[i];
+	if (header.compare("ATOM") == 0 && len >53){
+	  string strx,stry,strz;
+	  double x = 0, y = 0, z = 0;
+	  for (size_t j = 30; j!=38; ++j){
+	    strx += line[j];
+	    stry += line[j+8];
+	    strz += line[j+16];
+	  }
+	  x = atof(strx.c_str());
+	  y = atof(stry.c_str());
+	  z = atof(strz.c_str());
+	  
+	  //find residue type, and assign SLD
+	  string resname;
+	  for (size_t k = 17; k!=20; ++k) resname +=line[k];
+	  mymapitr itr = res_sld_.find("ALA");
+	  
+	  //apoint(x,y,z,sld)
+	  Point3D apoint(x,y,z,itr->second);
+	  vp.push_back(apoint);
+	  
+	  //save x y z' max & min to calculate the size boundary
+	  x_max_ = x > x_max_ ? x : x_max_;
+	  x_min_ = x < x_min_ ? x : x_min_;
+	  y_max_ = y > y_max_ ? y : y_max_;
+	  y_min_ = y < y_min_ ? y : y_min_;
+	  z_max_ = z > z_max_ ? z : z_max_;
+	  z_min_ = z < z_min_ ? z : z_min_;
+	}
+      }
+    }
+
+    infile.close();
+  }
+  return vp.size();
+}
+
+double PDBModel::GetDimBound()
+{
+  if (pdbnames_.size() == 0 )
+    return 0;
+
+  //cout <<x_max_<<" "<<x_min_<<" "<<y_max_<<" "<<y_min_<<" "<<z_max_<<" "<<z_min_<<endl;
+  double size = sqrt((x_max_-x_min_)*(x_max_-x_min_)+(y_max_-y_min_)*(y_max_-y_min_)+(z_max_-z_min_)*(z_max_-z_min_));
+  return size;
+}
+
+vector<double> PDBModel::GetCenter()
+{
+
+  vector<double> cen(3);
+  cen[0]= (x_max_+x_min_)/2;
+  cen[1]= (y_max_+y_min_)/2;
+  cen[2]= (z_max_+z_min_)/2;
+  center_ = cen;
+
+  return center_;
+}
Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/pdb_model.h
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/pdb_model.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/pdb_model.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,43 @@
+/** \file pdb_model.h child class of PointsModel */
+
+#ifndef PDB_MODEL_H
+#define PDB_MODEL_H
+
+#include "points_model.h"
+#include <string>
+#include <map>
+
+/**
+ *  Class PDBModel, PDB models
+ */
+
+class PDBModel : public PointsModel {
+ public:
+  PDBModel();
+  
+  //add pdb file
+  void AddPDB(const string &);
+  void AddPDB(const char*);
+
+  //Parse all coordinates from ATOM section 
+  //of the PDB file into vector of points
+  int GetPoints(Point3DVector &);
+
+  //will be used in determining the maximum distance for
+  //P(r) calculation for a complex model (merge several 
+  //pointsmodel instance together
+  vector<double> GetCenter();
+
+  //get the maximum possible dimension= dist between XYZmax XYZmin
+  //May improve to directly get the dimension from unit cell
+  //if the technique is crystallography
+  double GetDimBound();
+
+ private:
+  vector<string> pdbnames_;
+  map<string,double> res_sld_;
+  double x_max_,x_min_,y_max_,y_min_,z_max_,z_min_;
+};
+
+#endif
+
Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/points_model.cc
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/points_model.cc	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/points_model.cc	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,461 @@
+/** \file points_model.cc */
+
+#include <vector>
+#include <algorithm>
+#include <fstream>
+#include <stdio.h>
+//#include <exception>
+#include <stdexcept>
+#include "points_model.h"
+#include "Point3D.h"
+
+PointsModel::PointsModel()
+{
+  r_grids_num_ = 2000;
+  rmax_ = 0;
+  cormax_ = 0;
+  rstep_ = 0;
+}
+
+void PointsModel::CalculateIQ(IQ *iq)
+{
+  //fourier transform of the returned Array2D<double> from ddFunction()
+  int nIpoints = iq->GetNumI();
+  double qstep = (iq->GetQmax()) / (nIpoints-1);
+  vector<double> fint(nIpoints, 0);
+
+  //I(0) is calculated seperately
+  int num_rstep = pr_.dim1();
+
+  for (int k = 1; k<nIpoints; k++){
+
+    double q = k * qstep;
+
+    double r =0;
+    double debeye = 0;
+    double fadd =0;
+
+
+    for (int i = 1; i < num_rstep; ++i){
+      r = i*rstep_;  //r should start from 1* rstep
+      double qr = q*r;
+      debeye = sin(qr)/qr;
+      fadd = pr_[i][1]*debeye;
+      fint[k] = fint[k] + fadd;
+    }
+  }
+
+  //I(0)
+  double Izero = 0;
+  for (int i = 0; i < num_rstep; ++i)
+    Izero += pr_[i][1];
+  fint[0] = Izero;
+
+  //assign I(Q) with normalization
+  for(int j = 0; j < nIpoints; ++j){
+    (*iq).iq_data[j][0] = j * qstep;
+    (*iq).iq_data[j][1] = fint[j];
+    // remove normalization Izero;
+  }
+}
+
+//return I with a single q value
+double PointsModel::CalculateIQ(double q)
+{
+  //fourier transform of the returned Array2D<double> from ddFunction()
+  int num_rstep = pr_.dim1();
+
+  double r =0;
+  double debeye = 0;
+  double fadd = 0;
+  double Irelative = 0;
+
+  //I(0) is calculated seperately
+  if (q == 0){
+    //I(0)
+    double Izero = 0;
+    for (int i = 0; i < num_rstep; ++i)
+      Izero += pr_[i][1];
+    Irelative = Izero;
+  }
+  else {
+    for (int i = 1; i < num_rstep; ++i){
+      r = i*rstep_;  //r should start from 1* rstep
+      double qr = q*r;
+      debeye = sin(qr)/qr;
+      fadd = pr_[i][1]*debeye;
+      Irelative = Irelative + fadd;
+    }
+  }
+  return Irelative;
+}
+
+double PointsModel::CalculateIQError(double q)
+{
+  //fourier transform of the returned Array2D<double> from ddFunction()
+  int num_rstep = pr_.dim1();
+
+  double r =0;
+  double debeye = 0;
+  double fadd = 0;
+  double Irelative = 0;
+
+  //I(0) is calculated seperately
+  for (int i = 1; i < num_rstep; ++i){
+	  r = i*rstep_;  //r should start from 1* rstep
+	  double qr = q*r;
+	  debeye = sin(qr)/qr;
+	  fadd = fabs(pr_[i][2])*debeye*debeye 
+	  	+ rstep_*rstep_/4.0/r/r*(cos(qr)*cos(qr) + debeye*debeye);
+	  Irelative = Irelative + fadd;
+  }
+  return sqrt(Irelative);
+}
+
+//pass in a vector of points, and calculate the P(r)
+double PointsModel::DistDistribution(const vector<Point3D> &vp)
+{
+  //get r axis:0,rstep,2rstep,3rstep......d_bound
+  int sizeofpr = r_grids_num_ + 1; //+1 just for overflow prevention
+
+  double d_bound = GetDimBound();
+  rstep_ = CalculateRstep(r_grids_num_,d_bound);
+
+  Array2D<double> pr(sizeofpr, 3); //third column is left for error for the future
+  pr = 0;
+
+  for (int i = 1; i != sizeofpr; ++i)
+    pr[i][0] = pr[i-1][0] + rstep_ ;  //column 1:  distance
+
+  int size = vp.size();
+
+  for (int i1 = 0; i1 < size - 1; ++i1) {
+    for (int i2 = i1 + 1; i2 < size; ++i2) {
+      //dist_.push_back(vp[i1].distanceToPoint(vp[i2]));
+      //product_sld_.push_back(vp[i1].getSLD() * vp[i2].getSLD());
+      double a_dist = vp[i1].distanceToPoint(vp[i2]);
+      double its_sld = vp[i1].getSLD() * vp[i2].getSLD();
+
+      //save maximum distance
+      if (a_dist>rmax_) {
+        rmax_ = a_dist;
+      }
+      //insert into pr array
+      int l = int(floor(a_dist/rstep_));
+
+      //cout << "i1,i2,l,a_dist"<<vp[i1]<<" "<<vp[i2]<<" "<<l<<" "<<a_dist<<endl;     
+      //overflow check
+      if (l >= sizeofpr) {
+        cerr << "one distance is out of range: " << l <<endl;
+        //throw "Out of range";
+      }
+      else {
+        pr[l][1] += its_sld;  //column 2 intermediate: sum of SLD of the points with specific distance     
+        // Estimate uncertainty (squared)
+        pr[l][2] += its_sld*its_sld; 
+        //keep maxium Pr absolute number, in order to normalize
+        //if (pr[l][1] > cormax_)	  cormax_ = pr[l][1];
+      }
+    }
+  }
+  
+  //normalize Pr
+  for (int j = 0; j != sizeofpr; ++j){           //final column2 for P(r)
+    //pr[j][1] = pr[j][1]/cormax_;
+
+    // 'Size' is the number of space points, without double counting (excluding 
+    // overlapping regions between shapes). The volume of the combined shape
+    // is given by  V = size * (sum of all sub-volumes) / (Total number of points) 
+    //              V = size / (lores_density)
+    
+    // - To transform the integral to a sum, we need to give a weight
+    // to each entry equal to the average space volume of a point (w = V/N = 1/lores_density).
+    // The final output, I(q), should therefore be multiplied by V*V/N*N.
+    // Since we will be interested in P(r)/V, we only need to multiply by 1/N*(V/N) = 1/N/lores_density.
+    // We don't have access to lores_density from this class; we will therefore apply
+    // this correction externally.
+    //
+    // - Since the loop goes through half the points, multiply by 2.
+    // TODO: have access to lores_density from this class.
+    // 
+    pr[j][1] = 2.0*pr[j][1]/size;
+    pr[j][2] = 4.0*pr[j][2]/size/size;
+  }
+  pr_ = pr;
+
+  return rmax_;
+}
+
+Array2D<double> PointsModel::GetPr()
+{
+  return pr_;
+}
+
+
+double PointsModel::CalculateRstep(int num_grids, double rmax)
+{
+  assert(num_grids > 0);
+
+  double rstep;
+  rstep = rmax / num_grids;
+
+  return rstep;
+}
+
+void PointsModel::OutputPR(const string &fpr){
+  ofstream outfile(fpr.c_str());
+  if (!outfile) {
+    cerr << "error: unable to open output file: "
+	 << outfile << endl;
+    exit(1);
+  }
+  
+  double sum = 0.0;
+  double r_stepsize = 1.0;
+  if (pr_.dim1()>2) r_stepsize = pr_[1][0] - pr_[0][0];
+  
+  for (int i = 0;  i < pr_.dim1(); ++i){
+	  sum += pr_[i][1]*r_stepsize;
+  }
+  
+  for (int i = 0;  i < pr_.dim1(); ++i){
+	  if (pr_[i][1]==0) continue;
+	  outfile << pr_[i][0] << "       " << (pr_[i][1]/sum) << endl;
+  }
+}
+
+void PointsModel::OutputPDB(const vector<Point3D> &vp,const char *fpr){
+  FILE *outfile=NULL;
+  outfile = fopen(fpr,"w+");
+  if (!outfile) {
+    cerr << "error: unable to open output file: "
+	 << outfile << endl;
+    exit(1);
+  }
+  int size = vp.size();
+  int index = 0;
+  for (int i = 0;  i < size; ++i){
+    ++index;
+    fprintf(outfile,"ATOM%7d  C%24.3lf%8.3lf%8.3lf%6.3lf\n", \
+           index,vp[i].getX(),vp[i].getY(),vp[i].getZ(),vp[i].getSLD());
+  }
+  fclose(outfile);
+}
+
+PointsModel::~PointsModel()
+{
+}
+
+void PointsModel::DistDistributionXY(const vector<Point3D> &vp)
+{
+  //the max box get from 3D should be more than enough for 2D,but doesn't hurt
+  double d_bound = GetDimBound();
+
+  //using 1A for rstep, so the total bins is the max distance for the object
+  int sizeofpr = ceil(d_bound) + 1; //+1 just for overflow prevention
+  rstep_ = 1;
+
+  Array2D<double> pr_xy(sizeofpr,sizeofpr); //2D histogram
+
+  //the max frequency in the correlation histogram
+  double cormax_xy_ = 0;
+
+  //initialization
+  pr_xy = 0;
+
+  for (int i = 1; i != sizeofpr; ++i){
+    pr_xy[i][0] = pr_xy[i-1][0] + rstep_ ;  //column 1:  distance
+  }
+
+  int size = vp.size();
+
+  for (int i1 = 0; i1 < size - 1; ++i1) {
+    for (int i2 = i1 + 1; i2 < size; ++i2) {
+      int jx = int(floor(fabs(vp[i1].getX()-vp[i2].getX())/rstep_));
+      int jy = int(floor(fabs(vp[i1].getY()-vp[i2].getY())/rstep_));
+      //the sld for the pair of points
+      double its_sld = vp[i1].getSLD()*vp[i2].getSLD();
+
+      //overflow check
+      if ((jx >= sizeofpr) || (jy >= sizeofpr))
+      {
+        cerr << "one distance is out of range: " <<endl;
+        //throw "Out of range";
+      }
+      else{
+	pr_xy[jx][jy] += its_sld;
+        if (pr_xy[jx][jy] > cormax_xy_ ) cormax_xy_ = pr_xy[jx][jy];
+      }
+    }
+  }
+
+  //normalize Pr_xy
+  for (int m = 0; m != sizeofpr; ++m){           //final column2 for P(r)
+    for (int n = 0; n != sizeofpr; ++n){
+      pr_xy[m][n] = pr_xy[m][n]/cormax_xy_;
+      //cout << "m n:"<<m<<" "<<n<<" "<<pr_xy[m][n]<<endl;
+    }
+  }
+
+  pr_xy_ = pr_xy;
+}
+
+void PointsModel::OutputPR_XY(const std::string &fpr)
+{
+  ofstream outfile(fpr.c_str());
+  if (!outfile) {
+    cerr << "error: unable to open output file: "
+	 << outfile << endl;
+    exit(1);
+  }
+  int size = pr_xy_.dim1();
+  //pr_xy_ is a N x N array
+  for (int i = 0;  i != size; ++i){
+    for (int j = 0; j != size; ++j)
+    {
+      outfile << i << "    " << j <<"    "<< pr_xy_[i][j] << endl;
+    }
+  }
+}
+
+void PointsModel::CalculateIQ_2D(IQ *iq,double phi)
+{
+  int nIpoints = iq->GetNumI();
+  double qstep = (iq->GetQmax()) / (nIpoints-1);
+  vector<double> fint(nIpoints, 0);
+  double Izero = 0;
+
+  //number of bins on x and y axis
+  int size_r = pr_xy_.dim1();
+  //rstep is set to one, otherwise should be cos(phi)*rstep
+  double cosphi = cos(phi);
+  double sinphi = sin(phi);
+
+  for(int k = 1; k != nIpoints; ++k){
+    double q = k * qstep;
+    double tmp = cos(q*(cosphi+sinphi));
+
+    for(int i=0; i!=size_r; ++i){
+      for(int j = 0; j!=size_r; ++j){
+        fint[k] += pr_xy_[i][j]*tmp;
+      }
+    }
+  }
+
+  for(int i=0; i!=size_r; ++i){
+    for(int j = 0; j!=size_r; ++j){
+      Izero += pr_xy_[i][j];
+    }
+  }
+  fint[0] = Izero;
+
+  //assign I(Q) with normalization
+  for(int j = 0; j < nIpoints; ++j){
+    (*iq).iq_data[j][0] = j * qstep;
+    (*iq).iq_data[j][1] = fint[j] / Izero;
+  }
+}
+
+vector<double> PointsModel::GetCenter()
+{
+  vector<double> vp(3,0);
+  return vp;
+}
+
+double PointsModel::CalculateIQ_2D(double qx, double qy)
+{
+  //for each (Qx,Qy) on 2D detector, calculate I
+  double q = sqrt(qx*qx+qy*qy);
+  double I = 0;
+
+  double cosphi = qx/q;
+  double sinphi = qy/q;
+  double tmp = cos(q*(cosphi+sinphi));
+
+  //loop through P(r) on xy plane
+  int size_r = pr_xy_.dim1();
+  for(int i=-size_r+1; i!=size_r; ++i){
+    for(int j = -size_r+1; j!=size_r; ++j){
+      //rstep is set to one, left out from calculation
+      I += pr_xy_[abs(i)][abs(j)]*cos(q*(cosphi*i+sinphi*j));
+    }
+  }
+
+  //return I, without normalization
+  return I;
+}
+
+/*
+ * 2D simulation for oriented systems
+ * The beam direction is assumed to be in the z direction.
+ * 
+ * @param points: vector of space points
+ * @param qx: qx [A-1]
+ * @param qy: qy [A-1]
+ * @return: I(qx, qy) for the system described by the space points [cm-1]
+ * 
+ */
+double PointsModel::CalculateIQ_2D(const vector<Point3D>&points, double qx, double qy){
+	/*
+	 * TODO: the vector of points should really be part of the class
+	 * 	     This is a design flaw inherited from the original programmer.
+	 */
+	
+	int size = points.size();
+
+	double cos_term = 0;
+	double sin_term = 0;
+	for (int i = 0; i < size; i++) {
+	  	//the sld for the pair of points
+	
+		double phase = qx*points[i].getX() + qy*points[i].getY();
+			  	
+		cos_term += cos(phase) * points[i].getSLD();
+		sin_term += sin(phase) * points[i].getSLD();
+	
+	}   			
+
+	// P(q) = 1/V I(q) = (V/N)^2 (1/V) (cos_term^2 + sin_term^2) 
+	// We divide by N here and we will multiply by the density later.
+
+	return (cos_term*cos_term + sin_term*sin_term)/size;	
+}
+
+double PointsModel::CalculateIQ_2D_Error(const vector<Point3D>&points, double qx, double qy){
+	
+	int size = points.size();
+
+	double delta_x, delta_y;
+	double q_t2 = qx*qx + qy*qy;
+	double cos_term = 0;
+	double sin_term = 0;
+	double cos_err = 0;
+	double sin_err = 0;
+	
+	// Estimate the error on the position of each point
+	// in x or y as V^(1/3)/N
+	
+	for (int i = 0; i < size; i++) {
+		
+		
+		//the sld for the pair of points
+	
+		double phase = qx*points[i].getX() + qy*points[i].getY();
+		double sld_fac = points[i].getSLD() * points[i].getSLD();
+		
+		cos_term += cos(phase) * points[i].getSLD();
+		sin_term += sin(phase) * points[i].getSLD();
+
+		sin_err += cos(phase) * cos(phase) * sld_fac;
+		cos_err += sin(phase) * sin(phase) * sld_fac;
+	
+	}   			
+
+	// P(q) = 1/V I(q) = (V/N)^2 (1/V) (cos_term^2 + sin_term^2) 
+	// We divide by N here and we will multiply by the density later.
+
+	// We will need to multiply this error by V^(1/3)/N.
+	// We don't have access to V from within this class.
+	return 2*sqrt(cos_term*cos_term*cos_err*cos_err + sin_term*sin_term*sin_err*sin_err)/size;	
+}
+
Index: src/sas/simulation/pointsmodelpy/libpointsmodelpy/points_model.h
===================================================================
--- src/sas/simulation/pointsmodelpy/libpointsmodelpy/points_model.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/pointsmodelpy/libpointsmodelpy/points_model.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,61 @@
+/** \file points_model.h child class of SANSmodel*/
+
+#ifndef POINTSMODEL_H
+#define POINTSMODEL_H
+
+#include <vector>
+#include "tnt/tnt.h"
+#include <string>
+#include "Point3D.h"
+#include "iq.h"
+#include "sans_model.h"
+
+using namespace std;
+
+class PointsModel : public SANSModel{
+ public:
+  PointsModel();
+  virtual ~PointsModel();
+
+  void CalculateIQ(IQ *iq);
+  double CalculateIQ(double q);
+  double CalculateIQError(double q);
+  
+  // Old lengthy 2D simulation (unchecked) 
+  void CalculateIQ_2D(IQ *iq,double phi);
+  double CalculateIQ_2D(double qx, double qy);
+  
+  // Fast 2D simulation
+  double CalculateIQ_2D(const vector<Point3D>&, double qx, double qy);
+  double CalculateIQ_2D_Error(const vector<Point3D>&, double qx, double qy);
+  
+  //given a set of points, calculate distance correlation
+  //function, and return the max dist
+  double DistDistribution(const vector<Point3D>&);
+  void DistDistributionXY(const vector<Point3D>&);
+
+  Array2D<double> GetPr();
+  void OutputPR(const std::string &);
+  void OutputPR_XY(const std::string &);
+
+  virtual int GetPoints(Point3DVector &) = 0;
+  void OutputPDB(const vector<Point3D>&,const char*);
+
+  //will be used in calculating P(r), the maximum distance
+  virtual double GetDimBound() = 0;
+  //will be used to determin the maximum distance for
+  //several pointsmodel instances
+  virtual vector<double> GetCenter();
+
+ protected:
+  double CalculateRstep(int num_points, double rmax);
+  double rmax_, rstep_,cormax_,cormax_xy_;
+  int r_grids_num_;
+  vector<double> center_;
+
+ private:
+  Array2D<double> pr_,pr_xy_;
+  vector<double> product_sld_;
+};
+
+#endif
