Index: src/sans/simulation/geoshapespy/Make.mm
===================================================================
--- src/sans/simulation/geoshapespy/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,47 @@
+# -*- Makefile -*-
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+#
+# <LicenseText>
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+
+PROJECT = geoshapespy
+
+# directory structure
+
+BUILD_DIRS = \
+    libgeoshapespy \
+    geoshapespymodule \
+    geoshapespy \
+
+OTHER_DIRS = \
+    tests \
+    examples
+
+RECURSE_DIRS = $(BUILD_DIRS) $(OTHER_DIRS)
+
+#--------------------------------------------------------------------------
+#
+
+all:
+	BLD_ACTION="all" $(MM) recurse
+
+distclean::
+	BLD_ACTION="distclean" $(MM) recurse
+
+clean::
+	BLD_ACTION="clean" $(MM) recurse
+
+tidy::
+	BLD_ACTION="tidy" $(MM) recurse
+
+# version
+# $Id$
+
+# End of file
Index: src/sans/simulation/geoshapespy/examples/Make.mm
===================================================================
--- src/sans/simulation/geoshapespy/examples/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/examples/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,32 @@
+# -*- Makefile -*-
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+#
+# <LicenseText>
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+
+PROJECT = geoshapespy
+PACKAGE = examples
+
+#--------------------------------------------------------------------------
+#
+
+all: clean
+
+release: clean
+	cvs release .
+
+update: clean
+	cvs update .
+
+
+# version
+# $Id$
+
+# End of file
Index: src/sans/simulation/geoshapespy/geoshapespy/Make.mm
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespy/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespy/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,35 @@
+# -*- Makefile -*-
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+#
+# <LicenseText>
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+PROJECT = geoshapespy
+PACKAGE = geoshapespy
+
+#--------------------------------------------------------------------------
+#
+
+all: export
+
+
+#--------------------------------------------------------------------------
+#
+# export
+
+EXPORT_PYTHON_MODULES = \
+    __init__.py
+
+
+export:: export-python-modules
+
+# version
+# $Id$
+
+# End of file
Index: src/sans/simulation/geoshapespy/geoshapespy/__init__.py
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespy/__init__.py	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespy/__init__.py	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,21 @@
+#!/usr/bin/env python
+# 
+#  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# 
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+# 
+#  <LicenseText>
+# 
+#  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# 
+
+def copyright():
+    return "geoshapespy pyre module: Copyright (c) 1998-2005 Michael A.G. Aivazis";
+
+
+# version
+__id__ = "$Id$"
+
+#  End of file 
Index: src/sans/simulation/geoshapespy/geoshapespymodule/Make.mm
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespymodule/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespymodule/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,32 @@
+# -*- Makefile -*-
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+#
+# <LicenseText>
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+
+PROJECT = SANSsimulation
+PACKAGE = geoshapespymodule
+MODULE = geoshapespy
+
+include std-pythonmodule.def
+include local.def
+
+PROJ_CXX_SRCLIB = -lgeoshapespy
+
+PROJ_SRCS = \
+    bindings.cc \
+    exceptions.cc \
+    misc.cc
+
+
+# version
+# $Id$
+
+# End of file
Index: src/sans/simulation/geoshapespy/geoshapespymodule/bindings.cc
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespymodule/bindings.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespymodule/bindings.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,63 @@
+// -*- C++ -*-
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+//                               Michael A.G. Aivazis
+//                        California Institute of Technology
+//                        (C) 1998-2005  All Rights Reserved
+// 
+//  <LicenseText>
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+
+//#include <portinfo>
+#include <Python.h>
+
+#include "bindings.h"
+
+#include "misc.h"          // miscellaneous methods
+
+// the method table
+
+struct PyMethodDef pygeoshapespy_methods[] = {
+
+   //geoshapes methods: set_orientation; set_center
+   {pygeoshapespy_set_orientation__name__, pygeoshapespy_set_orientation,
+    METH_VARARGS, pygeoshapespy_set_orientation__doc__},
+
+   {pygeoshapespy_set_center__name__, pygeoshapespy_set_center,
+    METH_VARARGS, pygeoshapespy_set_center__doc__},
+
+   // new sphere
+   {pyanalmodelpy_new_sphere__name__, pyanalmodelpy_new_sphere,
+    METH_VARARGS, pyanalmodelpy_new_sphere__doc__},
+
+   // new cylinder
+   {pyanalmodelpy_new_cylinder__name__, pyanalmodelpy_new_cylinder,
+    METH_VARARGS, pyanalmodelpy_new_cylinder__doc__},
+
+   // new ellipsoid
+   {pyanalmodelpy_new_ellipsoid__name__, pyanalmodelpy_new_ellipsoid,
+    METH_VARARGS, pyanalmodelpy_new_ellipsoid__doc__},
+
+   // new hollowsphere
+   {pyanalmodelpy_new_hollowsphere__name__, pyanalmodelpy_new_hollowsphere,
+    METH_VARARGS, pyanalmodelpy_new_hollowsphere__doc__},
+
+   // new singlehelix
+   {pyanalmodelpy_new_singlehelix__name__, pyanalmodelpy_new_singlehelix,
+    METH_VARARGS, pyanalmodelpy_new_singlehelix__doc__},
+
+   {pygeoshapespy_copyright__name__, pygeoshapespy_copyright,
+    METH_VARARGS, pygeoshapespy_copyright__doc__},
+
+
+// Sentinel
+    {0, 0}
+};
+
+// version
+// $Id$
+
+// End of file
Index: src/sans/simulation/geoshapespy/geoshapespymodule/bindings.h
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespymodule/bindings.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespymodule/bindings.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,26 @@
+// -*- C++ -*-
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+//                               Michael A.G. Aivazis
+//                        California Institute of Technology
+//                        (C) 1998-2005  All Rights Reserved
+// 
+//  <LicenseText>
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+
+#if !defined(pygeoshapespy_bindings_h)
+#define pygeoshapespy_bindings_h
+
+// the method table
+
+extern struct PyMethodDef pygeoshapespy_methods[];
+
+#endif
+
+// version
+// $Id$
+
+// End of file
Index: src/sans/simulation/geoshapespy/geoshapespymodule/exceptions.cc
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespymodule/exceptions.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespymodule/exceptions.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,22 @@
+// -*- C++ -*-
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+//                               Michael A.G. Aivazis
+//                        California Institute of Technology
+//                        (C) 1998-2005  All Rights Reserved
+// 
+//  <LicenseText>
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+
+//#include <portinfo>
+#include <Python.h>
+
+PyObject *pygeoshapespy_runtimeError = 0;
+    
+// version
+// $Id$
+
+// End of file
Index: src/sans/simulation/geoshapespy/geoshapespymodule/exceptions.h
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespymodule/exceptions.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespymodule/exceptions.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,26 @@
+// -*- C++ -*-
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+//                               Michael A.G. Aivazis
+//                        California Institute of Technology
+//                        (C) 1998-2005  All Rights Reserved
+// 
+//  <LicenseText>
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+
+#if !defined(pygeoshapespy_exceptions_h)
+#define pygeoshapespy_exceptions_h
+
+// exceptions
+
+extern PyObject *pygeoshapespy_runtimeError;
+
+#endif
+
+// version
+// $Id$
+
+// End of file
Index: src/sans/simulation/geoshapespy/geoshapespymodule/geoshapespymodule.cc
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespymodule/geoshapespymodule.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespymodule/geoshapespymodule.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,44 @@
+// -*- C++ -*-
+
+#include <Python.h>
+
+
+#include "exceptions.h"
+#include "bindings.h"
+#include "myutil.h"
+
+
+char pygeoshapespy_module__doc__[] = "";
+
+// Initialization function for the module (*must* be called initgeoshapespy)
+extern "C"
+void
+initgeoshapespy()
+{
+    // create the module and add the functions
+    PyObject * m = Py_InitModule4(
+        "geoshapespy", pygeoshapespy_methods,
+        pygeoshapespy_module__doc__, 0, PYTHON_API_VERSION);
+
+    // get its dictionary
+    PyObject * d = PyModule_GetDict(m);
+
+    // check for errors
+    if (PyErr_Occurred()) {
+        Py_FatalError("can't initialize module geoshapespy");
+    }
+
+    // install the module exceptions
+    pygeoshapespy_runtimeError = PyErr_NewException("geoshapespy.runtime", 0, 0);
+    PyDict_SetItemString(d, "RuntimeException", pygeoshapespy_runtimeError);
+
+	// Seed the random number generator
+	seed_rnd();
+
+    return;
+}
+
+// version
+// $Id$
+
+// End of file
Index: src/sans/simulation/geoshapespy/geoshapespymodule/local.def
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespymodule/local.def	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespymodule/local.def	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,24 @@
+# -*- Makefile -*-
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+#
+# <LicenseText>
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+
+# C++
+
+    PROJ_CXX_INCLUDES = ../libgeoshapespy   \
+			../../iqPy/libiqPy  \
+			../../iqPy/libiqPy/tnt  \
+			../../pointsmodelpy/libpointsmodelpy
+
+# version
+# $Id$
+
+# End of file
Index: src/sans/simulation/geoshapespy/geoshapespymodule/misc.cc
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespymodule/misc.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespymodule/misc.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,190 @@
+// -*- C++ -*-
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+//                               Michael A.G. Aivazis
+//                        California Institute of Technology
+//                        (C) 1998-2005  All Rights Reserved
+// 
+//  <LicenseText>
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+
+//#include <portinfo>
+#include <Python.h>
+
+#include "misc.h"
+#include "sphere.h"
+#include "cylinder.h"
+#include "ellipsoid.h"
+#include "hollow_sphere.h"
+#include "single_helix.h"
+
+// copyright
+
+char pygeoshapespy_copyright__doc__[] = "";
+char pygeoshapespy_copyright__name__[] = "copyright";
+
+static char pygeoshapespy_copyright_note[] = 
+    "geoshapespy python module: Copyright (c) 1998-2005 Michael A.G. Aivazis";
+
+
+PyObject * pygeoshapespy_copyright(PyObject *, PyObject *)
+{
+    return Py_BuildValue("s", pygeoshapespy_copyright_note);
+}
+
+//GeoShape methods
+char pygeoshapespy_set_orientation__name__[] = "set_orientation";
+char pygeoshapespy_set_orientation__doc__[] = "Set the rotation angles";
+
+PyObject * pygeoshapespy_set_orientation(PyObject *, PyObject *args){
+  PyObject *pyshape = 0;
+  double angX=0,angY=0,angZ=0;
+  
+  int ok = PyArg_ParseTuple(args, "Oddd", &pyshape,&angX,&angY,&angZ); 
+  if(!ok) return NULL;
+
+  void *temp = PyCObject_AsVoidPtr(pyshape);
+
+  GeoShape *shape = static_cast<GeoShape *>(temp);
+
+  shape->SetOrientation(angX,angY,angZ);
+
+  return Py_BuildValue("i", 0);
+}
+
+char pygeoshapespy_set_center__name__[] = "set_center";
+char pygeoshapespy_set_center__doc__[] = "new center for points translation";
+
+PyObject * pygeoshapespy_set_center(PyObject *, PyObject *args){
+
+  PyObject *pyshape = 0;
+  double tranX=0,tranY=0,tranZ=0;
+  
+  int ok = PyArg_ParseTuple(args, "Oddd", &pyshape,&tranX,&tranY,&tranZ); 
+  if(!ok) return NULL;
+
+  void *temp = PyCObject_AsVoidPtr(pyshape);
+
+  GeoShape *shape = static_cast<GeoShape *>(temp);
+
+  shape->SetCenter(tranX,tranY,tranZ);
+
+  return Py_BuildValue("i", 0);
+}
+
+//Sphere constructor
+char pyanalmodelpy_new_sphere__name__[] = "new_sphere";
+char pyanalmodelpy_new_sphere__doc__[] = "sphere constructor";
+
+PyObject * pyanalmodelpy_new_sphere(PyObject *, PyObject *args){
+  double r;
+  int ok = PyArg_ParseTuple(args,"d",&r);
+  if(!ok) return 0;
+
+  Sphere *newsph = new Sphere(r);
+
+  return PyCObject_FromVoidPtr(newsph, PyDelSphere);
+}
+
+static void PyDelSphere(void *ptr){
+  Sphere * oldsph = static_cast<Sphere *>(ptr);
+  delete oldsph;
+
+  return;
+}
+
+//Cylinder constructor
+char pyanalmodelpy_new_cylinder__name__[] = "new_cylinder";
+char pyanalmodelpy_new_cylinder__doc__[] = "cylinder constructor";
+
+PyObject * pyanalmodelpy_new_cylinder(PyObject *, PyObject *args){
+  double r,h;
+  int ok = PyArg_ParseTuple(args,"dd",&r,&h);
+  if(!ok) return 0;
+
+  Cylinder *newcyl = new Cylinder(r,h);
+
+  return PyCObject_FromVoidPtr(newcyl, PyDelCylinder);
+}
+
+static void PyDelCylinder(void *ptr){
+  Cylinder * oldcyl = static_cast<Cylinder *>(ptr);
+  delete oldcyl;
+
+  return;
+}
+
+//Ellipsoid constructor
+char pyanalmodelpy_new_ellipsoid__name__[] = "new_ellipsoid";
+char pyanalmodelpy_new_ellipsoid__doc__[] = "ellipsoid constructor";
+
+PyObject * pyanalmodelpy_new_ellipsoid(PyObject *, PyObject *args){
+  double rx,ry,rz;
+  int ok = PyArg_ParseTuple(args,"ddd",&rx,&ry,&rz);
+  if(!ok) return 0;
+
+  Ellipsoid *newelli = new Ellipsoid(rx,ry,rz);
+
+  return PyCObject_FromVoidPtr(newelli, PyDelEllipsoid);
+}
+
+static void PyDelEllipsoid(void *ptr){
+  Ellipsoid * oldelli = static_cast<Ellipsoid *>(ptr);
+  delete oldelli;
+
+  return;
+}
+
+//Hollow Sphere constructor & methods
+char pyanalmodelpy_new_hollowsphere__name__[] = "new_hollowsphere";
+char pyanalmodelpy_new_hollowsphere__doc__[] = "";
+
+PyObject * pyanalmodelpy_new_hollowsphere(PyObject *, PyObject *args)
+{
+  double r, th;
+  int ok = PyArg_ParseTuple(args,"dd",&r, &th);
+  if(!ok) return 0;
+
+  HollowSphere *newhosph = new HollowSphere(r,th);
+
+  return PyCObject_FromVoidPtr(newhosph, PyDelHollowSphere);
+
+}
+
+static void PyDelHollowSphere(void *ptr)
+{
+  HollowSphere * oldhosph = static_cast<HollowSphere *>(ptr);
+  delete oldhosph;
+  return;
+}
+
+//Single Helix constructor & methods
+char pyanalmodelpy_new_singlehelix__name__[] = "new_singlehelix";
+char pyanalmodelpy_new_singlehelix__doc__[] = "";
+
+PyObject * pyanalmodelpy_new_singlehelix(PyObject *, PyObject *args)
+{
+  double hr,tr,pitch,turns;
+  int ok = PyArg_ParseTuple(args,"dddd",&hr,&tr,&pitch,&turns);
+  if(!ok) return 0;
+
+  SingleHelix *newsinhel = new SingleHelix(hr,tr,pitch,turns);
+
+  return PyCObject_FromVoidPtr(newsinhel, PyDelSingleHelix);
+
+}
+
+static void PyDelSingleHelix(void *ptr)
+{
+  SingleHelix * oldsinhel = static_cast<SingleHelix *>(ptr);
+  delete oldsinhel;
+  return;
+}
+    
+// version
+// $Id$
+
+// End of file
Index: src/sans/simulation/geoshapespy/geoshapespymodule/misc.h
===================================================================
--- src/sans/simulation/geoshapespy/geoshapespymodule/misc.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/geoshapespymodule/misc.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,79 @@
+// -*- C++ -*-
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+//                               Michael A.G. Aivazis
+//                        California Institute of Technology
+//                        (C) 1998-2005  All Rights Reserved
+// 
+//  <LicenseText>
+// 
+//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// 
+
+#if !defined(pygeoshapespy_misc_h)
+#define pygeoshapespy_misc_h
+
+// copyright
+extern char pygeoshapespy_copyright__name__[];
+extern char pygeoshapespy_copyright__doc__[];
+extern "C"
+PyObject * pygeoshapespy_copyright(PyObject *, PyObject *);
+
+//GeoShape methods
+extern char pygeoshapespy_set_orientation__name__[];
+extern char pygeoshapespy_set_orientation__doc__[];
+extern "C"
+PyObject * pygeoshapespy_set_orientation(PyObject *, PyObject *);
+
+extern char pygeoshapespy_set_center__name__[];
+extern char pygeoshapespy_set_center__doc__[];
+extern "C"
+PyObject * pygeoshapespy_set_center(PyObject *, PyObject *);
+
+//Sphere constructor & methods
+extern char pyanalmodelpy_new_sphere__name__[];
+extern char pyanalmodelpy_new_sphere__doc__[];
+extern "C"
+PyObject * pyanalmodelpy_new_sphere(PyObject *, PyObject *);
+
+static void PyDelSphere(void *);
+
+//Cylinder constructor & methods
+extern char pyanalmodelpy_new_cylinder__name__[];
+extern char pyanalmodelpy_new_cylinder__doc__[];
+extern "C"
+PyObject * pyanalmodelpy_new_cylinder(PyObject *, PyObject *);
+
+static void PyDelCylinder(void *);
+
+//Ellipsoid constructor & method
+extern char pyanalmodelpy_new_ellipsoid__name__[];
+extern char pyanalmodelpy_new_ellipsoid__doc__[];
+extern "C"
+PyObject * pyanalmodelpy_new_ellipsoid(PyObject *, PyObject *);
+
+static void PyDelEllipsoid(void *);
+
+//Hollow Sphere constructor & methods
+extern char pyanalmodelpy_new_hollowsphere__name__[];
+extern char pyanalmodelpy_new_hollowsphere__doc__[];
+extern "C"
+PyObject * pyanalmodelpy_new_hollowsphere(PyObject *, PyObject *);
+
+static void PyDelHollowSphere(void *);
+
+//Single Helix constructor & methods
+extern char pyanalmodelpy_new_singlehelix__name__[];
+extern char pyanalmodelpy_new_singlehelix__doc__[];
+extern "C"
+PyObject * pyanalmodelpy_new_singlehelix(PyObject *, PyObject *);
+
+static void PyDelSingleHelix(void *);
+
+#endif
+
+// version
+// $Id$
+
+// End of file
Index: src/sans/simulation/geoshapespy/libgeoshapespy/Make.mm
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,79 @@
+# -*- Makefile -*-
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+
+include local.def
+
+PROJECT = geoshapespy
+PACKAGE = libgeoshapespy
+
+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 = \
+    geo_shape.cc \
+    sphere.cc \
+    hollow_sphere.cc \
+    cylinder.cc   \
+    ellipsoid.cc  \
+    single_helix.cc \
+    Point3D.cc  \
+    myutil.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)
+
+# export
+export:: export-headers export-libraries
+
+endif
+
+EXPORT_HEADERS = \
+    geo_shape.h   \
+    sphere.h      \
+    hollow_sphere.h \
+    cylinder.h \
+    ellipsoid.h \
+    single_helix.h \
+    Point3D.h  \
+    myutil.h
+
+
+EXPORT_LIBS = $(PROJ_SAR)
+EXPORT_BINS = $(PROJ_DLL)
+
+
+# version
+# $Id$
+
+#
+# End of file
Index: src/sans/simulation/geoshapespy/libgeoshapespy/Point3D.cc
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/Point3D.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/Point3D.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,220 @@
+//Point3D.cpp  
+
+#include <cmath>
+#include <cstdlib>
+#include <stdexcept>
+#include "Point3D.h"
+
+using namespace std;
+
+Point3D::Point3D(double a, double b, double c, double sld)
+{
+  x = a;
+  y = b;
+  z = c;
+  sld_ =  sld;
+}
+
+ostream& operator<< (ostream& os, const Point3D& p)
+{
+  os << "(" << p.getX() << ", " << p.getY() << ", " << p.getZ() << "," << p.getSLD() << ")";
+
+  return os;
+}
+
+double Point3D::distanceToPoint(const Point3D &p) const
+{
+	double dx=x-p.x;
+	double dy=y-p.y;
+	double dz=z-p.z;
+	return sqrt(dx*dx+dy*dy+dz*dz);
+}
+
+double Point3D::distanceToLine(const Point3D &p1, const Point3D &p2, bool * pIsOutside /* 0 */) const
+{
+	double u = ((x-p1.x)*(p2.x-p1.x) + (y-p1.y)*(p2.y-p1.y) + (z-p1.z)*(p2.z-p1.z))/(p1.distanceToPoint(p2)*p1.distanceToPoint(p2));
+	
+	if(pIsOutside != 0) {
+		if ( u < 0 || u > 1)
+			*pIsOutside=true;
+		else 
+			*pIsOutside=false;
+	}
+
+	double interX=p1.x+u*(p2.x-p1.x);
+	double interY=p1.y+u*(p2.y-p1.y);
+	double interZ=p1.z+u*(p2.z-p1.z);
+	
+	return sqrt((x-interX)*(x-interX)+(y-interY)*(y-interY)+(z-interZ)*(z-interZ));
+}
+
+// p1 and p2 determine a line, they are two end points
+Point3D Point3D::getInterPoint(const Point3D &p1, const Point3D &p2, bool * pIsOutside /* 0 */) const
+{
+	double u = ((x-p1.x)*(p2.x-p1.x) + (y-p1.y)*(p2.y-p1.y) + (z-p1.z)*(p2.z-p1.z))/(p1.distanceToPoint(p2)*p1.distanceToPoint(p2));
+	
+	if(pIsOutside != 0) {
+
+		if ( u < 0 || u > 1)
+			*pIsOutside=true;
+		else 
+			*pIsOutside=false;
+	}
+
+	return Point3D(p1.x+u*(p2.x-p1.x), p1.y+u*(p2.y-p1.y), p1.z+u*(p2.z-p1.z));
+}
+
+void Point3D::set(double x1, double y1, double z1)
+{
+	x = x1;
+	y = y1;
+	z = z1;
+}
+
+double Point3D::norm() const
+{
+	return sqrt(x * x + y * y + z * z);
+}
+
+double Point3D::normalize()
+{
+	double v = norm();
+	
+	if(v != 0) {
+		x /= v;
+		y /= v;
+		z /= v;
+	}
+	
+	return v;
+}
+
+Point3D Point3D::normVector() const
+{
+	Point3D p;	
+	double v = norm();
+	
+	if(v != 0) {
+		p.x = x / v;
+		p.y = y / v;
+		p.z = z / v;
+	}
+
+	return p;
+}
+
+double Point3D::dotProduct(const Point3D &p) const
+{
+	return x * p.x + y * p.y + z * p.z;
+}
+
+Point3D Point3D::minus(const Point3D &p) const
+{
+	return Point3D(x - p.x, y - p.y, z - p.z);
+}
+
+Point3D Point3D::plus(const Point3D &p) const
+{
+	return Point3D(x + p.x, y + p.y, z + p.z);
+}
+
+Point3D& Point3D::operator=(const Point3D &p)
+{
+	x = p.x;
+	y = p.y;
+	z = p.z;
+	sld_ = p.sld_;
+
+	return *this;
+}
+
+void Point3D::scale(double s)
+{
+	x *= s;
+	y *= s;
+	z *= s;
+}
+
+Point3D Point3D::multiplyProduct(const Point3D &p)
+{
+	return Point3D(y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x);
+}
+
+void Point3D::Transform(const vector<double> &orien, const vector<double> &center){
+  if(orien[1]) RotateY(Degree2Radian(orien[1]));
+  if(orien[0]) RotateX(Degree2Radian(orien[0]));
+  if(orien[2]) RotateZ(Degree2Radian(orien[2]));
+
+  Translate(center[0],center[1],center[2]);
+}
+
+void Point3D::RotateX(const double ang_x)
+{
+  double sinA_ = sin(ang_x);
+  double cosA_ = cos(ang_x);
+
+  //x doesn't change;
+  double y_new= y*cosA_ - z*sinA_;
+  double z_new= y*sinA_ + z*cosA_;
+
+  y = y_new;
+  z = z_new;
+}
+
+void Point3D::RotateY(const double ang_y)
+{
+  double sinA_ = sin(ang_y);
+  double cosA_ = cos(ang_y);
+ 
+  double x_new = z*sinA_ + x*cosA_;
+  //y doesn't change
+  double z_new = z*cosA_ - x*sinA_;
+
+  x = x_new;
+  z = z_new;
+}
+
+void Point3D::RotateZ(const double ang_z)
+{
+  double sinA_ = sin(ang_z);
+  double cosA_ = cos(ang_z);
+
+  double x_new = x*cosA_ - y*sinA_;
+  double y_new = x*sinA_ + y*cosA_;
+  //z doesn't change
+
+  x = x_new;
+  y = y_new;
+}
+
+void Point3D::Translate(const double trans_x, const double trans_y, const double trans_z)
+{
+  double x_new = x + trans_x;
+  double y_new = y + trans_y;
+  double z_new = z + trans_z;
+
+  x = x_new;
+  y = y_new;
+  z = z_new;
+}
+
+double Point3D::Degree2Radian(const double degree)
+{
+  return degree/180*pi;
+}
+
+void Point3D::TransformMatrix(const vector<double> &rotmatrix, const vector<double> &center)
+{
+  if (rotmatrix.size() != 9 || center.size() != 3)
+    throw std::runtime_error("The size for rotation matrix vector has to be 9.");
+
+  double xold = x;
+  double yold = y;
+  double zold = z;
+
+  x = rotmatrix[0]*xold + rotmatrix[1]*yold + rotmatrix[2]*zold;
+  y = rotmatrix[3]*xold + rotmatrix[4]*yold + rotmatrix[5]*zold;
+  z = rotmatrix[6]*xold + rotmatrix[7]*yold + rotmatrix[8]*zold;
+
+  Translate(center[0],center[1],center[2]);
+}
Index: src/sans/simulation/geoshapespy/libgeoshapespy/Point3D.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/Point3D.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/Point3D.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,87 @@
+//class point3D    1/29/2004 Jing
+//properties: copy point coordinates
+//	      calculate point to point distance 
+//	      shortest distance from a point to a line 
+
+
+#ifndef __POINT3D_
+#define __POINT3D_
+
+#include <iostream>
+#include <vector>
+
+using namespace std;
+
+const double pi = 3.1415926;
+
+class Point3D
+{
+  public:
+        Point3D() {}
+
+	//assign values with SLD to a point
+	Point3D(double a, double b, double c, double sld = 0);
+
+	// output
+	friend std::ostream& operator<<(std::ostream&, const Point3D&);
+
+	//distance to a point
+  	double distanceToPoint(const Point3D &p) const;
+  	
+  	//distance to a line
+  	double distanceToLine(const Point3D &p1, const Point3D &p2, bool *pv = 0) const; 
+
+	// get the point lying on the axis from a point
+	Point3D getInterPoint(const Point3D &p1, const Point3D &p2, bool *pv = 0) const;
+
+	// normalization
+	Point3D normVector() const;
+
+	// get length
+	double norm() const;
+	
+	double normalize();
+	
+	// assignment operator
+	Point3D& operator=(const Point3D &p);
+	
+	// scale this point with s
+	void scale(double s);
+	
+	// multiplication product of two vectors
+	Point3D multiplyProduct(const Point3D &p);
+
+	// p0 - p
+	Point3D minus(const Point3D &p) const;
+
+	// p0 + p
+	Point3D plus(const Point3D &p) const;
+
+	// dot product
+	double dotProduct(const Point3D &p) const;
+
+	//if you do not care if the point falls into the range of line,you can directly use distanceToLine(p1,p2)
+	void set(double x1, double y1, double z1);
+
+	double getX() const { return x; }
+	double getY() const { return y; }
+	double getZ() const { return z; }
+	double getSLD() const { return sld_; }
+
+        //Transformation
+        void Transform(const vector<double> &orien, const vector<double> &center);
+        void TransformMatrix(const vector<double> &rotmatrix, const vector<double> &center);
+  private:
+  	double x, y, z;
+	double sld_;
+
+        void RotateX(const double ang_x);
+        void RotateY(const double ang_y);
+        void RotateZ(const double ang_z);
+        void Translate(const double trans_x, const double trans_y, const double trans_z);
+        double Degree2Radian(const double degree);
+};
+
+typedef std::vector<Point3D> Point3DVector;
+
+#endif
Index: src/sans/simulation/geoshapespy/libgeoshapespy/cylinder.cc
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/cylinder.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/cylinder.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,108 @@
+/** \file Cylinder.cc */
+#include <cmath>
+#include <cassert>
+#include "cylinder.h"
+#include "myutil.h"
+
+Cylinder::Cylinder()
+{
+  r_ = 0;
+}
+
+Cylinder::Cylinder(double radius,double length)
+{
+  r_ = radius;
+  l_ = length;
+  topcenter_ = Point3D(0,length/2,0);
+  botcenter_ = Point3D(0,-length/2,0);
+}
+
+Cylinder::~Cylinder()
+{
+}
+
+void Cylinder::SetRadius(double r)
+{
+  r_ = r;
+}
+
+void Cylinder::SetLength(double l)
+{
+  l_ = l;
+}
+
+double Cylinder::GetRadius()
+{
+  return r_;
+}
+
+double Cylinder::GetLength()
+{
+  return l_;
+}
+
+double Cylinder::GetMaxRadius()
+{
+  double maxr = sqrt(4*r_*r_ + l_*l_)/2;
+  return maxr;
+}
+
+ShapeType Cylinder::GetShapeType() const
+{
+  return CYLINDER;
+}
+
+double Cylinder::GetVolume()
+{
+  double V = pi * square(r_) * l_;
+  return V;
+}
+
+void Cylinder::GetFormFactor(IQ * iq)
+{
+  /** number of I for output, equal to the number of rows of array IQ*/
+  /** to be finished */
+}
+
+Point3D Cylinder::GetAPoint(double sld)
+{
+  /** cylinder long axis is along Y to match vtk actor */
+  static int max_try = 100;
+  for (int i = 0; i < max_try; ++i) {
+    double x = (ran1()-0.5) * 2 * r_;
+    double z = (ran1()-0.5) * 2 * r_;	
+    double y = (ran1()-0.5) * l_;
+
+    Point3D apoint(x,y,z,sld);
+    //check the cross section on xy plane within a sphere at (0,)
+    if (apoint.distanceToPoint(Point3D(0,y,0)) <= r_ ) 
+      return apoint;
+  }
+
+  std::cerr << "Max try "
+            << max_try
+            << " is reached while generating a point in cylinder" << std::endl;
+  return Point3D(0, 0, 0);
+}
+
+bool Cylinder::IsInside(const Point3D& point) const
+{
+  bool isOutside = false;
+  double distToLine = point.distanceToLine(GetTopCenter(),GetBotCenter(),&isOutside);
+
+  return (distToLine <= r_ && !isOutside);
+}
+
+Point3D Cylinder::GetTopCenter() const
+{
+  Point3D new_center(topcenter_);
+  new_center.Transform(GetOrientation(),GetCenter());
+  return new_center;
+}
+
+Point3D Cylinder::GetBotCenter() const
+{
+  Point3D new_center(botcenter_);
+  new_center.Transform(GetOrientation(),GetCenter());
+  return new_center;
+}
Index: src/sans/simulation/geoshapespy/libgeoshapespy/cylinder.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/cylinder.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/cylinder.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,63 @@
+/** \file cylinder.h class Cylinder:GeoShape */
+
+#ifndef Cylinder_H
+#define Cylinder_H
+
+#include "geo_shape.h"
+
+/** class Cylinder, subclass of GeoShape */
+class Cylinder : public GeoShape {
+ public:
+
+  /** initialize */
+  Cylinder();
+
+  /** constructor with radius initialization */
+  Cylinder(double radius,double length);
+
+  ~Cylinder();
+
+  /** set parameter radius */
+  void SetRadius(double r);
+
+  /** set parameter length */
+  void SetLength(double l);
+
+  /** get the radius */
+  double GetRadius();
+  
+  /** get the length */
+  double GetLength();
+
+  /** get the radius of the sphere to cover this shape */
+  double GetMaxRadius();
+
+  /** get the volume */
+  double GetVolume();
+
+  /** calculate the cylinder form factor, no scale, no background*/
+  void GetFormFactor(IQ * iq);
+
+  /** using a equation to check whether a point with XYZ lies
+      within the cylinder with center (0,0,0)
+  */
+  Point3D GetAPoint(double sld);
+
+  /** check whether a point is inside the cylinder at any position 
+      in the 3D space
+  */
+  bool IsInside(const Point3D& point) const;
+
+  ShapeType GetShapeType() const;
+
+ protected:
+  Point3D GetTopCenter() const;
+  Point3D GetBotCenter() const;
+
+ private:
+  double r_, l_;
+  Point3D topcenter_,botcenter_;
+
+};
+
+#endif
Index: src/sans/simulation/geoshapespy/libgeoshapespy/ellipsoid.cc
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/ellipsoid.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/ellipsoid.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -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;
+}
Index: src/sans/simulation/geoshapespy/libgeoshapespy/ellipsoid.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/ellipsoid.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/ellipsoid.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,65 @@
+/** \file cylinder.h class Cylinder:GeoShape */
+
+#ifndef Ellipsoid_H
+#define Ellipsoid_H
+
+#include "geo_shape.h"
+
+/** class Ellipsoid, subclass of GeoShape */
+class Ellipsoid : public GeoShape {
+ public:
+
+  /** initialize */
+  Ellipsoid();
+
+  /** constructor with radii initialization */
+  Ellipsoid(double rx,double ry, double rz);
+
+  ~Ellipsoid();
+
+  /** set parameter radii */
+  void SetRadii(double rx, double ry, double rz);
+
+  /** get the radii */
+  double GetRadiusX();
+  double GetRadiusY();
+  double GetRadiusZ();
+  
+  /** get the radius of the sphere to cover this shape */
+  double GetMaxRadius();
+
+  /** get the volume */
+  double GetVolume();
+
+  /** calculate the cylinder form factor, no scale, no background*/
+  void GetFormFactor(IQ * iq);
+
+  /** using a equation to check whether a point with XYZ lies
+      within the cylinder with center (0,0,0)
+  */
+  Point3D GetAPoint(double sld);
+
+  /** check whether a point is inside the cylinder at any position 
+      in the 3D space
+  */
+  bool IsInside(const Point3D& point) const;
+
+  ShapeType GetShapeType() const;
+
+ protected:
+  //get a point on three axis for function IsInside,using the edge point
+  Point3D GetXaxisPlusEdge() const;
+  Point3D GetXaxisMinusEdge() const;
+  Point3D GetYaxisPlusEdge() const;
+  Point3D GetYaxisMinusEdge() const;
+  Point3D GetZaxisPlusEdge() const;
+  Point3D GetZaxisMinusEdge() const;
+
+ private:
+  double rx_,ry_,rz_;
+  Point3D xedge_plus_,xedge_minus_;
+  Point3D yedge_plus_,yedge_minus_;
+  Point3D zedge_plus_,zedge_minus_;
+};
+
+#endif
Index: src/sans/simulation/geoshapespy/libgeoshapespy/geo_shape.cc
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/geo_shape.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/geo_shape.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,26 @@
+#include "geo_shape.h"
+
+void GeoShape::SetOrientation(double angX,double angY, double angZ){
+  orientation_[0]=angX;
+  orientation_[1]=angY;
+  orientation_[2]=angZ;
+}
+
+void GeoShape::SetCenter(double cenX,double cenY,double cenZ){
+  center_[0]=cenX;
+  center_[1]=cenY;
+  center_[2]=cenZ;
+}
+
+vector<double> GeoShape::GetOrientation() const {
+  return orientation_;
+}
+
+vector<double> GeoShape::GetCenter() const {
+  return center_;
+}
+
+Point3D GeoShape::GetCenterP() const {
+  Point3D p(center_[0],center_[1],center_[2]);
+  return p;
+}
Index: src/sans/simulation/geoshapespy/libgeoshapespy/geo_shape.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/geo_shape.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/geo_shape.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,66 @@
+/*! \file a abstract class GeoShape 
+*/
+#ifndef GEOSHAPE_H
+#define GEOSHAPE_H
+
+#include "iq.h"
+#include "Point3D.h"
+#include <vector>
+
+using namespace std;
+
+#define triple(x) ((x) * (x) * (x))
+#define square(x) ((x) * (x))
+
+enum ShapeType{ SPHERE, HOLLOWSPHERE, CYLINDER , ELLIPSOID,SINGLEHELIX};
+
+/**
+ *  class GeoShape, abstract class, parent class for Sphere, Cylinder ....
+ */
+
+class GeoShape{
+ public:
+
+  GeoShape(){
+    vector<double> init(3,0.0);
+    orientation_ = init;
+    center_ = init;
+  }
+  
+  virtual ~GeoShape() {}
+
+  /** calculate the form factor for a simple shape */
+  virtual void GetFormFactor(IQ * iq) = 0;
+
+  /** Get a point that is within the simple shape*/
+  virtual Point3D GetAPoint(double sld) = 0;
+
+  /** check whether a point is inside the shape*/
+  virtual bool IsInside(const Point3D& point) const = 0;
+
+  virtual double GetVolume() = 0;
+
+  virtual ShapeType GetShapeType() const = 0;
+
+  /** get the radius of the sphere to cover the shape*/
+  virtual double GetMaxRadius() = 0;
+
+  void SetOrientation(double angX, double angY, double angZ);
+
+  void SetCenter(double cenX, double cenY, double cenZ);
+
+  vector<double> GetOrientation() const;
+
+  vector<double> GetCenter() const;
+
+  Point3D GetCenterP() const;
+
+ private:
+  vector<double> orientation_;
+  
+ protected:
+  vector<double> center_;
+ 
+};
+
+#endif
Index: src/sans/simulation/geoshapespy/libgeoshapespy/hollow_sphere.cc
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/hollow_sphere.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/hollow_sphere.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,83 @@
+/** \file hollowsphere.cc  */
+#include <cmath>
+#include <cassert>
+#include "hollow_sphere.h"
+
+HollowSphere::HollowSphere()
+{
+  ro_ = 0;
+  th_ = 0;
+}
+
+HollowSphere::HollowSphere(double radius, double thickness)
+{
+  ro_ = radius;
+  th_ = thickness;
+}
+
+double HollowSphere::GetMaxRadius()
+{
+  double maxr = ro_;
+  return maxr;
+}
+
+void HollowSphere::GetFormFactor(IQ * iq)
+{
+  /** number of I for output, equal to the number of rows of array IQ*/
+  int numI = iq->iq_data.dim1();
+  double qmin = iq->GetQmin();
+  double qmax = iq->GetQmax();
+
+  assert(numI > 0);
+  assert(qmin > 0);
+  assert(qmax > 0);
+  assert( qmax > qmin );
+
+  double logMin = log10(qmin);
+  double z = logMin;
+  double logMax = log10(qmax);
+  double delta = (logMax - logMin) / (numI-1);
+
+  //not finished yet, need to find the part which is equal to 1
+  for(int i = 0; i < numI; ++i) {
+
+    /** temp for Q*/
+    double q = pow(z,10);
+
+    double ri_ = ro_ - th_ ;
+
+    double bes1 = 3.0 * (sin(q*ri_) - q*ri_*cos(q*ri_)) / triple(q) / triple(ri_);
+    double bes2 = 3.0 * (sin(q*ro_) - q*ro_*cos(q*ro_)) / triple(q) / triple(ro_);
+    double bes = (triple(ro_)*bes1 - triple(ri_)*bes2)/(triple(ro_) - triple(ri_));
+    /** double f is the temp for I, should be equal to one when q is 0*/
+    double f = bes * bes;
+
+    /** IQ[i][0] is Q,Q starts from qmin (non zero),q=0 handle separately IQ[i][1] is I */
+    iq->iq_data[i][0]= q;
+    iq->iq_data[i][1]= f;
+
+    z += delta;
+  }
+
+}
+
+Point3D HollowSphere::GetAPoint(double sld)
+{
+  return Point3D(0,0,0);
+}
+
+double HollowSphere::GetVolume()
+{
+  return 0;
+}
+
+bool HollowSphere::IsInside(const Point3D& point) const
+{
+  return true;
+}
+
+ShapeType HollowSphere::GetShapeType() const
+{
+  return HOLLOWSPHERE;
+}
+
Index: src/sans/simulation/geoshapespy/libgeoshapespy/hollow_sphere.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/hollow_sphere.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/hollow_sphere.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,43 @@
+/** \file hollowsphere.h class HollowSphere:GeoShape */
+
+#ifndef HOLLOWSPHERE_H
+#define HOLLOWSPHERE_H
+
+#include "geo_shape.h"
+
+/** class Sphere, subclass of GeoShape */
+class HollowSphere : public GeoShape {
+ public:
+
+  /** initialize */
+  HollowSphere();
+
+  /** constructor with radius&length initialization */
+  HollowSphere(double radius, double thickness);
+
+  /** get the radius of the sphere to cover this shape */
+  double GetMaxRadius();
+
+  /** calculate the sphere form factor, no scale, no background*/
+  void GetFormFactor(IQ * iq);
+
+  /** using a equation to check whether a point with XYZ lies
+      within the sphere with center (0,0,0)
+  */
+  Point3D GetAPoint(double sld);
+
+  /** check whether point is inside hollow sphere*/
+  bool IsInside(const Point3D& point) const;
+
+  double GetVolume();
+
+  //get the shape type, return sphere 
+  ShapeType GetShapeType() const;
+
+ private:
+  //outer radius, and the thickness of the shell
+  double ro_, th_;
+
+};
+
+#endif
Index: src/sans/simulation/geoshapespy/libgeoshapespy/local.def
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/local.def	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/local.def	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,30 @@
+# -*- 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  \
+		     ../../pointsmodelpy/libpointsmodelpy 
+
+    PROJ_CXX_FLAGS += $(CXX_SOFLAGS)
+    PROJ_LCXX_FLAGS += $(LCXX_SARFLAGS) -liqPy
+
+# version
+# $Id$
+
+# End of file
Index: src/sans/simulation/geoshapespy/libgeoshapespy/minmax.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/minmax.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/minmax.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,23 @@
+#ifndef GUARD_minmax_H
+#define GUARD_minmax_H
+#ifdef _MSC_VER
+// needed to cope with bug in MS library:
+// it fails to define min/max
+
+template <class T> inline T max(const T& a, const T& b)
+{
+
+      return (a > b) ? a : b;
+
+}
+
+template <class T> inline T min(const T& a, const T& b)
+{
+
+      return (a < b) ? a : b;
+
+}
+
+#endif
+
+#endif
Index: src/sans/simulation/geoshapespy/libgeoshapespy/myutil.cc
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/myutil.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/myutil.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,16 @@
+/** \file myutil.cc */
+#include <cstdlib>
+#include <ctime> 
+#include "myutil.h"
+
+using namespace std;
+
+void seed_rnd() {
+	srand((unsigned)time(NULL)); 
+}
+
+double ran1()
+{
+	return (double)rand() / (double) RAND_MAX;
+}
+
Index: src/sans/simulation/geoshapespy/libgeoshapespy/myutil.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/myutil.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/myutil.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,9 @@
+/** \file myutil.h */
+#ifndef _MYUTIL_H
+#define _MYUTIL_H
+
+/** return a random float number between [0,1] */
+double ran1();
+void seed_rnd();
+
+#endif
Index: src/sans/simulation/geoshapespy/libgeoshapespy/single_helix.cc
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/single_helix.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/single_helix.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,140 @@
+/** \file SingleHelix.cc */
+#include <cmath>
+#include <cassert>
+#include "single_helix.h"
+#include "myutil.h"
+
+SingleHelix::SingleHelix()
+{
+  hr_ = 0;
+  tr_ = 0;
+  pitch_ = 0;
+  turns_ = 0;
+}
+
+SingleHelix::SingleHelix(double helix_radius,double tube_radius,double pitch,double turns)
+{
+  hr_ = helix_radius;
+  tr_ = tube_radius;
+  pitch_ = pitch;
+  turns_ =  turns;
+}
+
+SingleHelix::~SingleHelix()
+{
+}
+
+void SingleHelix::SetHelixRadius(double hr)
+{
+  hr_ = hr;
+}
+
+void SingleHelix::SetTubeRadius(double tr)
+{
+  tr_ = tr;
+}
+
+void SingleHelix::SetPitch(double pitch)
+{
+  pitch_ = pitch;
+}
+
+void SingleHelix::SetTurns(double turns)
+{
+  turns_ = turns;
+}
+
+double SingleHelix::GetHelixRadius()
+{
+  return hr_;
+}
+
+double SingleHelix::GetTubeRadius()
+{
+  return tr_;
+}
+
+double SingleHelix::GetPitch()
+{
+  return pitch_;
+}
+
+double SingleHelix::GetTurns()
+{
+  return turns_;
+}
+
+double SingleHelix::GetMaxRadius()
+{
+  // put helix into a cylinder
+  double r_ = hr_ + tr_;
+  double l_ = pitch_*turns_ + 2*tr_;
+  double maxr = sqrt(4*r_*r_ + l_*l_)/2;
+  return maxr;
+}
+
+ShapeType SingleHelix::GetShapeType() const
+{
+  return SINGLEHELIX;
+}
+
+double SingleHelix::GetVolume()
+{
+  double V = pi*square(tr_)*sqrt(square(2*pi*hr_)+square(turns_*pitch_));
+  return V;
+}
+
+void SingleHelix::GetFormFactor(IQ * iq)
+{
+  /** number of I for output, equal to the number of rows of array IQ*/
+  /** to be finished */
+}
+
+Point3D SingleHelix::GetAPoint(double sld)
+{
+  int max_try = 100;
+  int i = 0;
+  double point1 = 0, point2 = 0, point3 = 0;
+ 
+  do{
+    i++;
+    if (i > max_try)
+	break;
+    zr_ = tr_*sqrt(1+square((pitch_/(2*pi*hr_))));
+    point1 = (ran1()-0.5)*2*tr_;
+    point2 = (ran1()-0.5)*2*zr_;
+    point3 = (ran1()-0.5)*4*pi*turns_;
+  } while (((square(point1/hr_)+square(point2/zr_))>1) || (point2+pitch_*point3/(2*pi)<0));
+
+  double x = (point1 + hr_)*cos(point3);
+  double y = (point1 + hr_)*sin(point3);
+  //"-pitch_*turns_/2" is just to corretly center the helix
+  double z = point2 + pitch_*point3/(2*pi)-pitch_*turns_/2;
+
+  Point3D apoint(x,y,z,sld);
+  return apoint;
+
+  
+  std::cerr << "Max try "
+            << max_try
+            << " is reached while generating a point in cylinder" << std::endl;
+  
+  return Point3D(0,0,0);
+}
+
+bool SingleHelix::IsInside(const Point3D& point) const
+{
+  double x = point.getX()-center_[0];
+  double y = point.getY()-center_[1];
+  double z = point.getZ()-center_[2];
+
+  double p3 = atan(y/x);
+  double p2 = z-(pitch_*p3)/(2*pi);
+  double p1 = x/cos(p3)-hr_;
+
+  if ((square(p1/tr_)+square(p2/zr_))>1 || p2+pitch_*p3/(2*pi))
+    return false;
+  else
+    return true;
+}
+
Index: src/sans/simulation/geoshapespy/libgeoshapespy/single_helix.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/single_helix.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/single_helix.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,75 @@
+/** \file cylinder.h class Cylinder:GeoShape */
+
+#ifndef SingleHelix_H
+#define SingleHelix_H
+
+#include "geo_shape.h"
+
+/** class SingleHelix, subclass of GeoShape */
+class SingleHelix : public GeoShape {
+ public:
+
+  /** initialize */
+  SingleHelix();
+
+  /** constructor with radius initialization */
+  SingleHelix(double helix_radius,double tube_radius,double pitch,double turns);
+
+  ~SingleHelix();
+
+  /** set parameter helix_radius */
+  void SetHelixRadius(double hr);
+
+  /** set parameter tube_radius */
+  void SetTubeRadius(double tr);
+
+  /** set parameter pitch */
+  void SetPitch(double pitch);
+
+  /** set parameter number of turns */
+  void SetTurns(double turns);
+
+  /** get the helix_radius */
+  double GetHelixRadius();
+
+  /** get the tube_radius */
+  double GetTubeRadius();
+  
+  /** get the pitch */
+  double GetPitch();
+
+  /** get the number of turns */
+  double GetTurns();
+
+  /** get the radius of the sphere to cover this shape */
+  double GetMaxRadius();
+
+  /** get the volume */
+  double GetVolume();
+
+  /** calculate the single helix form factor, no scale, no background*/
+  void GetFormFactor(IQ * iq);
+
+  /** using a equation to check whether a point with XYZ lies
+      within the cylinder with center (0,0,0)
+  */
+  Point3D GetAPoint(double sld);
+
+  /** check whether a point is inside the cylinder at any position 
+      in the 3D space
+  */
+  bool IsInside(const Point3D& point) const;
+
+  ShapeType GetShapeType() const;
+
+ protected:
+  //Point3D GetTopCenter() const;
+  //Point3D GetBotCenter() const;
+
+ private:
+  double hr_,tr_, pitch_, turns_,zr_;
+  Point3D topcenter_,botcenter_;
+
+};
+
+#endif
Index: src/sans/simulation/geoshapespy/libgeoshapespy/sphere.cc
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/sphere.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/sphere.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,106 @@
+/** \file Sphere.cc */
+#include <cmath>
+#include <cassert>
+#include "sphere.h"
+#include "myutil.h"
+
+Sphere::Sphere()
+{
+  r_ = 0;
+}
+
+Sphere::Sphere(double radius)
+{
+  r_ = radius;
+}
+
+Sphere::~Sphere()
+{
+}
+
+void Sphere::SetRadius(double r)
+{
+  r_ = r;
+}
+
+double Sphere::GetRadius()
+{
+  return r_;
+}
+
+double Sphere::GetMaxRadius()
+{
+  double maxr = r_;
+  return maxr;
+}
+
+ShapeType Sphere::GetShapeType() const
+{
+  return SPHERE;
+}
+
+double Sphere::GetVolume()
+{
+  double V = 4 * pi / 3 * triple(r_);
+  return V;
+}
+
+void Sphere::GetFormFactor(IQ * iq)
+{
+  /** number of I for output, equal to the number of rows of array IQ*/
+  int numI = iq->iq_data.dim1();
+  double qmin = iq->GetQmin();
+  double qmax = iq->GetQmax();
+
+  assert(numI > 0);
+  assert(qmin > 0);
+  assert(qmax > 0);
+  assert( qmax > qmin );
+  
+  double logMin = log10(qmin);
+  double z = logMin;
+  double logMax = log10(qmax);
+  double delta = (logMax - logMin) / (numI-1);
+
+  for(int i = 0; i < numI; ++i) {
+
+    /** temp for Q*/
+    double q = pow(z, 10);
+ 
+    double bes = 3.0 * (sin(q*r_) - q*r_*cos(q*r_)) / triple(q) / triple(r_);
+
+    /** double f is the temp for I, should be equal to one when q is 0*/
+    double f = bes * bes;
+
+    /** IQ[i][0] is Q,Q starts from qmin (non zero),q=0 handle separately IQ[i][1] is IQ*/
+    iq->iq_data[i][0]= q;
+    iq->iq_data[i][1]= f;
+
+    z += delta;
+  }  
+}
+
+Point3D Sphere::GetAPoint(double sld)
+{
+  static int max_try = 100;
+  for (int i = 0; i < max_try; ++i) {
+    double x = (ran1()-0.5) * 2 * r_;
+    double y = (ran1()-0.5) * 2 * r_;	
+    double z = (ran1()-0.5) * 2 * r_;
+
+    Point3D apoint(x,y,z,sld);
+    if (apoint.distanceToPoint(Point3D(0,0,0)) <=r_ ) //dist to origin give sphere shape
+      return apoint;
+  }
+
+  std::cerr << "Max try "
+            << max_try
+            << " is reached while generating a point in sphere" << std::endl;
+  return Point3D(0, 0, 0);
+}
+
+bool Sphere::IsInside(const Point3D& point) const
+{
+  return point.distanceToPoint(Point3D(center_[0], center_[1], center_[2])) <= r_;
+  cout << "distance = " << point.distanceToPoint(Point3D(center_[0], center_[1], center_[2])) << endl;
+}
Index: src/sans/simulation/geoshapespy/libgeoshapespy/sphere.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/sphere.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/sphere.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,52 @@
+/** \file sphere.h class Sphere:GeoShape */
+
+#ifndef SPHERE_H
+#define SPHERE_H
+
+#include "geo_shape.h"
+
+/** class Sphere, subclass of GeoShape */
+class Sphere : public GeoShape {
+ public:
+
+  /** initialize */
+  Sphere();
+
+  /** constructor with radius initialization */
+  Sphere(double radius);
+
+  ~Sphere();
+
+  /** set parameter radius */
+  void SetRadius(double r);
+
+  /** get the radius */
+  double GetRadius();
+
+  /** Get the radius of the sphere to cover this shape */
+  double GetMaxRadius();
+
+  /** get the volume */
+  double GetVolume();
+
+  /** calculate the sphere form factor, no scale, no background*/
+  void GetFormFactor(IQ * iq);
+
+  /** using a equation to check whether a point with XYZ lies
+      within the sphere with center (0,0,0)
+  */
+  Point3D GetAPoint(double sld);
+
+  /** check whether a point is inside the sphere at any position
+      in the 3D space
+  */
+  bool IsInside(const Point3D& point) const;
+
+  ShapeType GetShapeType() const;
+
+ private:
+  double r_;
+
+};
+
+#endif
Index: src/sans/simulation/geoshapespy/libgeoshapespy/transformation.cc
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/transformation.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/transformation.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,73 @@
+#include "transformation.h"
+#include <math.h>
+#include <stdexcept>
+
+void RotateX(const double ang_x,Point3D &a_point)
+{
+  double sinA_ = sin(ang_x);
+  double cosA_ = cos(ang_x);
+  double pointX_ = a_point.getX();
+  double pointY_ = a_point.getY();
+  double pointZ_ = a_point.getZ();
+
+  double x_new_ = pointX_;
+  double y_new_ = pointY_*cosA_ - pointZ_*sinA_;
+  double z_new_ = pointY_*sinA_ - pointZ_*cosA_;
+
+  a_point.set(x_new_, y_new_, z_new_);
+}
+
+void RotateY(const double ang_y,Point3D &a_point)
+{
+  double sinA_ = sin(ang_y);
+  double cosA_ = cos(ang_y);
+  double pointX_ = a_point.getX();
+  double pointY_ = a_point.getY();
+  double pointZ_ = a_point.getZ();
+
+  double x_new_ = pointX_*cosA_ + pointZ_*sinA_;
+  double y_new_ = pointY_;
+  double z_new_ = -pointX_*sinA_ + pointZ_*cosA_;
+
+  a_point.set(x_new_, y_new_, z_new_);
+}
+
+void RotateZ(const double ang_z,Point3D &a_point)
+{
+  double sinA_= sin(ang_z);
+  double cosA_ = cos(ang_z);
+  double pointX_ = a_point.getX();
+  double pointY_ = a_point.getY();
+  double pointZ_ = a_point.getZ();
+
+  double x_new_ = pointX_*cosA_ - pointY_*sinA_;
+  double y_new_ = pointX_*sinA_+ pointY_*cosA_;
+  double z_new_ = pointZ_;
+
+  a_point.set(x_new_, y_new_, z_new_);
+}
+
+void Translate(const double trans_x, const double trans_y, const double trans_z, Point3D &a_point)
+{
+  double x_new_ = a_point.getX() + trans_x;
+  double y_new_ = a_point.getY() + trans_y;
+  double z_new_ = a_point.getZ() + trans_z;
+
+  a_point.set(x_new_, y_new_, z_new_);
+}
+
+void RotateMatrix(const vector<double> &rotmatrix, Point3D &a_point)
+{
+  if (rotmatrix.size() != 9)
+    throw std::runtime_error("The size for rotation matrix vector has to be 9.");
+
+  double xold = a_point.getX();
+  double yold = a_point.getY();
+  double zold = a_point.getZ();
+
+  double x_new_ = rotmatrix[0]*xold + rotmatrix[1]*yold + rotmatrix[2]*zold;
+  double y_new_ = rotmatrix[3]*xold + rotmatrix[4]*yold + rotmatrix[5]*zold;
+  double z_new_ = rotmatrix[6]*xold + rotmatrix[7]*yold + rotmatrix[8]*zold;
+
+  a_point.set(x_new_, y_new_, z_new_);
+}
Index: src/sans/simulation/geoshapespy/libgeoshapespy/transformation.h
===================================================================
--- src/sans/simulation/geoshapespy/libgeoshapespy/transformation.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/libgeoshapespy/transformation.h	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,21 @@
+/** \file transformation.h: a few functions to do the transformation */
+
+#ifndef TRANSFORMATION_H
+//Transform operation on a Point3D
+//Not in use, similar transformation are implemented as Point3D class functions
+
+#define TRANSFORMATION_H
+
+#include <vector>
+#include "Point3D.h"
+
+void RotateX(const double ang_x,Point3D &);
+void RotateY(const double ang_y,Point3D &);
+void RotateZ(const double ang_z,Point3D &);
+void Translate(const double trans_x, const double trans_y, const double trans_z, Point3D &);
+
+//rotate a point by given a rotation 3X3 matrix
+//vec[0]=R00, vec[1]=R01.......vec[8]=R22
+void RotateMatrix(const vector<double> &, Point3D &);
+
+#endif
Index: src/sans/simulation/geoshapespy/tests/Make.mm
===================================================================
--- src/sans/simulation/geoshapespy/tests/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/tests/Make.mm	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,59 @@
+# -*- Makefile -*-
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+#
+# <LicenseText>
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+PROJECT = geoshapespy
+PACKAGE = tests
+
+PROJ_CLEAN += $(PROJ_CPPTESTS)
+
+PROJ_CXX_INCLUDES += ../../iqPy/libiqPy   \
+                     ../../iqPy/libiqPy/tnt \
+                     ../libgeoshapespy   \
+		     ../../pointsmodelpy/libpointsmodelpy
+
+PROJ_PYTESTS = 
+PROJ_CPPTESTS = testpoint
+
+PROJ_TESTS = $(PROJ_PYTESTS) $(PROJ_CPPTESTS)
+PROJ_LIBRARIES = -L$(BLD_LIBDIR) -lgeoshapespy -liqPy -lpointsmodelpy
+
+
+#--------------------------------------------------------------------------
+#
+
+all: $(PROJ_TESTS)
+
+test:
+	for test in $(PROJ_TESTS) ; do $${test}; done
+
+release: tidy
+	cvs release .
+
+update: clean
+	cvs update .
+
+#--------------------------------------------------------------------------
+#
+
+testshapes: testshapes.cc $(BLD_LIBDIR)/libgeoshapespy.$(EXT_SAR)
+	$(CXX) $(CXXFLAGS) $(LCXXFLAGS) -o $@ testshapes.cc $(PROJ_LIBRARIES)
+
+testpoint: testPoint.cc $(BLD_LIBDIR)/libgeoshapespy.$(EXT_SAR)
+	$(CXX) $(CXXFLAGS) $(LCXXFLAGS) -o $@ testPoint.cc $(PROJ_LIBRARIES)
+
+testtransformation:  testtransformation.cc $(BLD_LIBDIR)/libgeoshapespy.$(EXT_SAR)
+	$(CXX) $(CXXFLAGS) $(LCXXFLAGS) -o $@ testtransformation.cc $(PROJ_LIBRARIES)
+
+# version
+# $Id$
+
+# End of file
Index: src/sans/simulation/geoshapespy/tests/testPoint.cc
===================================================================
--- src/sans/simulation/geoshapespy/tests/testPoint.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/tests/testPoint.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,40 @@
+#include <iostream>
+#include <vector>
+#include "Point3D.h"
+
+
+using namespace std;
+
+int main(){
+
+  cout << "test 1:initialize a point, set orientation, set center"
+       << " and perform transform" <<endl;
+  Point3D apoint(3,4,5);
+
+  vector<double> orient(3),center(3);
+  for (int i = 0; i!=3; ++i){
+    center[i] = 0;
+  }
+  orient[0]=10;
+  orient[1]=0;
+  orient[2]=30;
+
+  apoint.Transform(orient,center);
+
+  cout << apoint.getX() <<" "<<apoint.getY()<<" "<<apoint.getZ()<<endl;
+  cout << orient[0] <<" "<<orient[1]<<" "<<orient[2]<<endl;
+
+  cout << "test 2: initialize a point, and perform transformation \n"
+       << "through a rotation matrix and a translation" <<endl;
+
+  vector<double> rotmatrix(9,1);
+  for (size_t i = 0; i != rotmatrix.size(); ++i){
+    cout << rotmatrix[i] <<" ";
+  }
+
+  Point3D pp(1,1,1);
+  pp.TransformMatrix(rotmatrix, center);
+
+  cout << pp.getX() <<" "<<pp.getY()<<" "<<pp.getZ()<<endl;
+  return 0;
+}
Index: src/sans/simulation/geoshapespy/tests/testorientation.cc
===================================================================
--- src/sans/simulation/geoshapespy/tests/testorientation.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/tests/testorientation.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,13 @@
+#include "sphere.h"
+#include <iostream>
+
+int main(){
+  Sphere sph1(10);
+  sph1.SetOrientation(10,20,30);
+  sph1.SetCenter(1,5,10);
+
+  vector<double> v1 = sph1.GetOrientation();
+  cout << v1[0] <<endl;
+
+  return 0;
+}
Index: src/sans/simulation/geoshapespy/tests/testshapes.cc
===================================================================
--- src/sans/simulation/geoshapespy/tests/testshapes.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/tests/testshapes.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,171 @@
+#include <cassert>
+#include <iostream>
+#include <vector>
+#include "sphere.h"
+#include "cylinder.h"
+#include "hollow_sphere.h"
+#include "ellipsoid.h"
+#include "single_helix.h"
+#include "iq.h"
+
+using namespace std;
+
+void TestGetFormFactor_sphere() {
+  Sphere sphere(1.0);
+
+  IQ iq1(10,0.001, 0.3);
+  sphere.GetFormFactor(&iq1);
+
+  for (int i = 0; i< iq1.iq_data.dim1(); i++)
+    cout << iq1.iq_data[i][0]<< " " << iq1.iq_data[i][1] <<endl;
+}
+
+void TestGetShapeType_sphere() {
+  Sphere sphere;
+  //if (assert(sphere.GetShapeType() == SPHERE))
+  //  cout << " it is a sphere" << endl;
+  cout << sphere.GetShapeType() << endl;
+}
+
+void TestGetShapeType(GeoShape &gs){
+  cout <<gs.GetShapeType() << endl;
+}
+
+void TestGetVolume_sphere() {
+  Sphere s(1);
+  cout << "volume is " << s.GetVolume();
+}
+
+void TestGetAPoint_sphere(){
+  Sphere asphere(10);
+  cout << asphere.GetAPoint(1) << endl;
+}
+
+void TestGetFormFactor_HS(){
+
+  HollowSphere hs(10, 2);
+
+  IQ iq1(10, 0.001, 0.3);
+  hs.GetFormFactor(&iq1);
+
+  for (int i = 0; i< iq1.iq_data.dim1(); i++)
+    cout << iq1.iq_data[i][0]<< " " << iq1.iq_data[i][1] <<endl;
+
+
+}
+
+void TestGetShapeType_HS(){
+  HollowSphere hs1;
+  assert(hs1.GetShapeType() == HOLLOWSPHERE);
+}
+
+void TestIsInside_sphere(){
+  cout << "Testing IsInside(): ";
+  Sphere s1(10);
+  Point3D pp(20,30,30);
+  Point3D pp2(1,1,1);
+
+  cout << s1.IsInside(pp) <<" " << s1.IsInside(pp2) << " ";
+  cout << "Finished testing IsInside()" <<endl;
+}
+
+void TestCylinder(){
+  Cylinder c1(10,20);
+  Point3D p1 = c1.GetAPoint(1);
+  Point3D p4 = c1.GetAPoint(1);
+  cout << "Got two points inside the cylinder" << endl;
+  cout << "the distance is: " << p1.distanceToPoint(p4);
+  cout << p1 << endl;
+  cout << "is this point inside?:" << c1.IsInside(p1) << endl;
+  
+  Point3D p2(20,20,20);
+  cout << "this point should be outside:" << c1.IsInside(p2) <<endl;
+
+  c1.SetOrientation(10,10,10);
+  c1.SetCenter(20,20,20);
+
+  vector<double> ori;
+  ori.push_back(10);
+  ori.push_back(10);
+  ori.push_back(20);
+
+  vector<double> cen;
+  cen.push_back(20);
+  cen.push_back(20);
+  cen.push_back(20);
+
+  p1.Transform(ori,cen);
+  cout << "is the p1 still inside after orientation?" << c1.IsInside(p1) << endl;
+  cout << "p1 p2 distance: " <<p1.distanceToPoint(p4) << endl;
+
+  Point3D p3(0,0,0);
+  cout <<" this origin point should be outside after transformation:" << c1.IsInside(p3) <<endl;
+
+  Point3D pp1(6.8,-2.11,11);
+  cout <<"is pp1 inside? "<< c1.IsInside(pp1) << endl;
+
+}
+
+void TestEllipsoid(){
+  Ellipsoid e1(5,10,15);
+  cout << "generated an ellipsoid:"<<e1.GetRadiusX()<<" ";
+  cout <<e1.GetRadiusY()<<" " << e1.GetRadiusZ() << endl;
+  cout << "maximum radius is: " << e1.GetMaxRadius() << endl;
+  cout << "volume is : " << e1.GetVolume() << endl;
+  cout << "Get a point: " << e1.GetAPoint(1.0) << endl;
+
+  Point3D p1(1,2,3);
+  Point3D p2(6,11,16);
+  Point3D p3(100,100,100);
+  cout <<"is " << p1 << " inside? " << e1.IsInside(p1) << endl;
+  cout <<"is " << p2 << " inside? " << e1.IsInside(p2) << endl;
+
+  e1.SetCenter(10,10,10);
+  Point3D p4(11,12,13);
+  cout <<"is " << p2 << " inside after moved? " << e1.IsInside(p2) << endl;
+
+  e1.SetOrientation(90,90,90);
+  cout <<"is " << p3 << " inside after rotated? " << e1.IsInside(p3) << endl;
+
+}
+
+void TestSingleHelix(){
+
+  cout <<" Testing single helix" << endl;
+  SingleHelix sh(10,3,30,1);
+  cout << "Get a point" << sh.GetAPoint(1.0) << endl;  
+
+  Point3D p1(1,1,-50); //cannot be in sh(10,3,30,1)
+  cout <<"is " << p1 <<" inside?" << sh.IsInside(p1)<< endl;
+  Point3D p2(-3.15922, -11.6193, 20.5302); //should be inside
+  cout <<"is " << p2 <<" inside?" << sh.IsInside(p2)<< endl;
+
+}
+
+int main(){
+ 
+  //cout << "testing sphere :\n" ;
+  //TestGetAPoint_sphere();
+  //TestGetShapeType_sphere();
+  //testing getting shape type from parent
+  //Sphere ss(1);
+  //TestGetShapeType(ss);
+  //TestGetVolume_sphere();
+  //TestIsInside_sphere();
+
+  //cout << "testing hollowsphere: \n" ;	
+  //TestGetFormFactor_HS();
+  //TestGetShapeType_HS();
+  //  TestGetVolume();
+  
+  //cout << "testing cylinder:\n";
+  //TestCylinder();
+  
+  //cout << "testing ellipsoid :\n";
+  //TestEllipsoid();
+
+  TestSingleHelix();
+  printf("PASS.\n");
+
+  return 0;
+}
Index: src/sans/simulation/geoshapespy/tests/testshapes.py
===================================================================
--- src/sans/simulation/geoshapespy/tests/testshapes.py	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/tests/testshapes.py	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+# 
+#  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# 
+#                               Michael A.G. Aivazis
+#                        California Institute of Technology
+#                        (C) 1998-2005  All Rights Reserved
+# 
+#  <LicenseText>
+# 
+#  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# 
+
+if __name__ == "__main__":
+
+    from SANSsimulation import geoshapespy
+
+    print
+    print "module information:"
+    print "    file:", geoshapespy.__file__
+    print "    doc:", geoshapespy.__doc__
+    print "    contents:", dir(geoshapespy)
+
+    sp = geoshapespy.new_sphere(10)
+#    geoshapespy.set_orientation(sp,10,20,10)
+    cy = geoshapespy.new_cylinder(2,6)
+
+    el = geoshapespy.new_ellipsoid(25,15,10)
+
+    hs = geoshapespy.new_hollowsphere(10,2)
+
+    sh = geoshapespy.new_singlehelix(10,2,30,2)
+
+# version
+__id__ = "$Id$"
+
+#  End of file 
Index: src/sans/simulation/geoshapespy/tests/testsphere.cc
===================================================================
--- src/sans/simulation/geoshapespy/tests/testsphere.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
+++ src/sans/simulation/geoshapespy/tests/testsphere.cc	(revision aa639ea269b76d78cd120b41705f70280ef4e03e)
@@ -0,0 +1,36 @@
+#include <cassert>
+#include <iostream>
+#include "sphere.h"
+#include "iq.h"
+
+using namespace std;
+
+void TestGetAnalyticalIQ() {
+  Sphere sphere(1.0);
+
+  IQ iq1(10,0.001, 0.3);
+  sphere.GetFormFactor(&iq1);
+
+  for (int i = 0; i< iq1.iq_data.dim1(); i++)
+    cout << iq1.iq_data[i][0]<< " " << iq1.iq_data[i][1] <<endl;
+}
+
+void TestGetShapeType() {
+  Sphere sphere;
+  assert(sphere.GetShapeType() == SPHERE);
+}
+
+void TestGetVolume() {
+  Sphere s(1);
+  cout << "volume is " << s.GetVolume();
+}
+
+int main(){
+  TestGetAnalyticalIQ();
+  TestGetShapeType();
+  TestGetVolume();
+
+  printf("PASS.\n");
+
+  return 0;
+}
