Index: sansdataloader/src/sans/dataloader/data_info.py
===================================================================
--- sansdataloader/src/sans/dataloader/data_info.py	(revision e38e3350180e57eace423b477cb12450a3d10b75)
+++ sansdataloader/src/sans/dataloader/data_info.py	(revision f60a8c20e2f0018d19feff6fc41e188d2939a7e4)
@@ -1,21 +1,18 @@
 """
-    Module that contains classes to hold information read from 
+    Module that contains classes to hold information read from
     reduced data files.
     
-    A good description of the data members can be found in 
+    A good description of the data members can be found in
     the CanSAS 1D XML data format:
     
     http://www.smallangles.net/wgwiki/index.php/cansas1d_documentation
 """
-
-
+#####################################################################
 #This software was developed by the University of Tennessee as part of the
 #Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
-#project funded by the US National Science Foundation. 
-#If you use DANSE applications to do scientific research that leads to 
-#publication, we ask that you acknowledge the use of the software with the 
-#following sentence:
-#This work benefited from DANSE software developed under NSF award DMR-0520547.
+#project funded by the US National Science Foundation.
+#See the license text in license.txt
 #copyright 2008, University of Tennessee
+######################################################################
 
 
@@ -28,4 +25,5 @@
 import numpy
 import math
+
 
 class plottable_1D:
@@ -54,8 +52,8 @@
         self.y = numpy.asarray(y)
         if dx is not None:
-             self.dx = numpy.asarray(dx)
-        if dy is not None: 
+            self.dx = numpy.asarray(dx)
+        if dy is not None:
             self.dy = numpy.asarray(dy)
-        if dxl is not None: 
+        if dxl is not None:
             self.dxl = numpy.asarray(dxl)
         if dxw is not None: 
@@ -75,4 +73,5 @@
         self._yaxis = label
         self._yunit = unit
+
 
 class plottable_2D:
@@ -85,11 +84,11 @@
     ymax = None
     data = None
-    qx_data     = None
-    qy_data     = None
-    q_data      = None
-    err_data    = None
-    dqx_data    = None
-    dqy_data    = None
-    mask        = None
+    qx_data = None
+    qy_data = None
+    q_data = None
+    err_data = None
+    dqx_data = None
+    dqy_data = None
+    mask = None
     
     # Units
@@ -102,6 +101,6 @@
     
     def __init__(self, data=None, err_data=None, qx_data=None,
-                  qy_data=None, q_data=None, mask=None,
-                   dqx_data=None, dqy_data=None):
+                 qy_data=None, q_data=None, mask=None,
+                 dqx_data=None, dqy_data=None):
         self.data = numpy.asarray(data)
         self.qx_data = numpy.asarray(qx_data)
@@ -110,6 +109,8 @@
         self.mask = numpy.asarray(mask)
         self.err_data = numpy.asarray(err_data)
-        if dqx_data is not None: self.dqx_data = numpy.asarray(dqx_data) 
-        if dqy_data is not None: self.dqy_data = numpy.asarray(dqy_data) 
+        if dqx_data is not None:
+            self.dqx_data = numpy.asarray(dqx_data) 
+        if dqy_data is not None:
+            self.dqy_data = numpy.asarray(dqy_data) 
                
     def xaxis(self, label, unit):
@@ -174,6 +175,6 @@
     distance = None
     distance_unit = 'mm'
-    ## Offset of this detector position in X, Y, 
-    #(and Z if necessary) [Vector] [mm] 
+    ## Offset of this detector position in X, Y,
+    #(and Z if necessary) [Vector] [mm]
     offset = None
     offset_unit = 'm'
@@ -182,5 +183,5 @@
     orientation = None
     orientation_unit = 'degree'
-    ## Center of the beam on the detector in X and Y 
+    ## Center of the beam on the detector in X and Y
     #(and Z if necessary) [Vector] [mm]
     beam_center = None
@@ -204,5 +205,4 @@
         self.pixel_size  = Vector()
         
-    
     def __str__(self):
         _str  = "Detector:\n"
@@ -222,4 +222,5 @@
         return _str
 
+
 class Aperture:
     ## Name
@@ -239,4 +240,5 @@
         self.size = Vector()
     
+    
 class Collimation:
     """
@@ -265,8 +267,9 @@
         return _str
 
+
 class Source:
     """
     Class to hold source information
-    """  
+    """
     ## Name
     name = None
@@ -296,5 +299,4 @@
         self.beam_size = Vector()
         
-    
     def __str__(self):
         _str  = "Source:\n"
@@ -314,5 +316,5 @@
     
     
-""" 
+"""
 Definitions of radiation types
 """
@@ -321,4 +323,5 @@
 MUON     = 'muon'
 ELECTRON = 'electron'
+    
     
 class Sample:
@@ -371,4 +374,5 @@
         return _str
   
+  
 class Process:
     """
@@ -405,5 +409,5 @@
     the data itself and any other meta data.
     """
-    ## Title 
+    ## Title
     title      = ''
     ## Run number
@@ -460,5 +464,5 @@
         self.meta_data  = {}
         ## Loading errors
-        self.errors = []        
+        self.errors = []
         
     def append_empty_process(self):
@@ -497,5 +501,5 @@
     # but should be implemented for each data class inherited from DataInfo
     # that holds actual data (ex.: Data1D)
-    def _perform_operation(self, other, operation): 
+    def _perform_operation(self, other, operation):
         """
         Private method to perform operation. Not implemented for DataInfo,
@@ -619,5 +623,6 @@
         def operation(a, b):
             return b/a
-        return self._perform_operation(other, operation)            
+        return self._perform_operation(other, operation)
+            
             
 class Data1D(plottable_1D, DataInfo):
@@ -631,5 +636,4 @@
         DataInfo.__init__(self)
         plottable_1D.__init__(self, x, y, dx, dy)
-        
         
     def __str__(self):
@@ -654,5 +658,5 @@
         
         """
-        def _check(v):            
+        def _check(v):
             if (v.__class__ == list or v.__class__ == numpy.ndarray) \
                 and len(v) > 0 and min(v) > 0:
@@ -675,8 +679,8 @@
         
         if clone is None or not issubclass(clone.__class__, Data1D):
-            x  = numpy.zeros(length) 
-            dx = numpy.zeros(length) 
-            y  = numpy.zeros(length) 
-            dy = numpy.zeros(length) 
+            x  = numpy.zeros(length)
+            dx = numpy.zeros(length)
+            y  = numpy.zeros(length)
+            dy = numpy.zeros(length)
             clone = Data1D(x, y, dx=dx, dy=dy)
         
@@ -684,13 +688,13 @@
         clone.run         = self.run
         clone.filename    = self.filename
-        clone.instrument   = self.instrument
-        clone.notes       = deepcopy(self.notes) 
-        clone.process     = deepcopy(self.process) 
-        clone.detector    = deepcopy(self.detector) 
-        clone.sample      = deepcopy(self.sample) 
-        clone.source      = deepcopy(self.source) 
-        clone.collimation = deepcopy(self.collimation) 
-        clone.meta_data   = deepcopy(self.meta_data) 
-        clone.errors      = deepcopy(self.errors) 
+        clone.instrument  = self.instrument
+        clone.notes       = deepcopy(self.notes)
+        clone.process     = deepcopy(self.process)
+        clone.detector    = deepcopy(self.detector)
+        clone.sample      = deepcopy(self.sample)
+        clone.source      = deepcopy(self.source)
+        clone.collimation = deepcopy(self.collimation)
+        clone.meta_data   = deepcopy(self.meta_data)
+        clone.errors      = deepcopy(self.errors)
         
         return clone
@@ -716,5 +720,5 @@
             if len(self.x) != len(other.x) or \
                 len(self.y) != len(other.y):
-                msg =  "Unable to perform operation: data length are not equal"
+                msg = "Unable to perform operation: data length are not equal"
                 raise ValueError, msg
             
@@ -734,5 +738,5 @@
         dy = self.dy
         if self.dy == None or (len(self.dy) != len(self.y)):
-            dy = numpy.zeros(len(self.y))            
+            dy = numpy.zeros(len(self.y))
             
         return dy, dy_other
@@ -761,4 +765,5 @@
         return result
         
+        
 class Data2D(plottable_2D, DataInfo):
     """
@@ -777,18 +782,17 @@
     y_bins = None
     
-    
     def __init__(self, data=None, err_data=None, qx_data=None,
-                  qy_data=None, q_data=None, mask=None, 
-                  dqx_data=None, dqy_data=None):
+                 qy_data=None, q_data=None, mask=None,
+                 dqx_data=None, dqy_data=None):
         self.y_bins = []
         self.x_bins = []
         DataInfo.__init__(self)
         plottable_2D.__init__(self, data, err_data, qx_data,
-                              qy_data, q_data,mask, dqx_data, dqy_data)
+                              qy_data, q_data, mask, dqx_data, dqy_data)
         if len(self.detector) > 0:
             raise RuntimeError, "Data2D: Detector bank already filled at init"
 
     def __str__(self):
-        _str =  "%s\n" % DataInfo.__str__(self)
+        _str = "%s\n" % DataInfo.__str__(self)
         
         _str += "Data:\n"
@@ -814,7 +818,7 @@
         from copy import deepcopy
         
-        if clone is None or not issubclass(clone.__class__, Data2D):  
-            data     = numpy.zeros(length) 
-            err_data = numpy.zeros(length) 
+        if clone is None or not issubclass(clone.__class__, Data2D):
+            data = numpy.zeros(length)
+            err_data = numpy.zeros(length)
             qx_data = numpy.zeros(length)
             qy_data = numpy.zeros(length)
@@ -824,5 +828,5 @@
             dqy_data = None
             clone = Data2D(data, err_data, qx_data, qy_data,
-                            q_data,mask, dqx_data=dqx_data, dqy_data=dqy_data)
+                            q_data, mask, dqx_data=dqx_data, dqy_data=dqy_data)
 
         clone.title       = self.title
@@ -830,15 +834,14 @@
         clone.filename    = self.filename
         clone.instrument  = self.instrument
-        clone.notes       = deepcopy(self.notes) 
-        clone.process     = deepcopy(self.process) 
-        clone.detector    = deepcopy(self.detector) 
-        clone.sample      = deepcopy(self.sample) 
-        clone.source      = deepcopy(self.source) 
-        clone.collimation = deepcopy(self.collimation) 
-        clone.meta_data   = deepcopy(self.meta_data) 
-        clone.errors      = deepcopy(self.errors) 
+        clone.notes       = deepcopy(self.notes)
+        clone.process     = deepcopy(self.process)
+        clone.detector    = deepcopy(self.detector)
+        clone.sample      = deepcopy(self.sample)
+        clone.source      = deepcopy(self.source)
+        clone.collimation = deepcopy(self.collimation)
+        clone.meta_data   = deepcopy(self.meta_data)
+        clone.errors      = deepcopy(self.errors)
         
         return clone
-  
   
     def _validity_check(self, other):
@@ -865,5 +868,5 @@
                
             # Check that the scales match
-            #TODO: matching scales?     
+            #TODO: matching scales?
             
             # Check that the other data set has errors, otherwise
@@ -885,5 +888,4 @@
         return err, err_other
   
-  
     def _perform_operation(self, other, operation):
         """
@@ -897,5 +899,5 @@
         dy, dy_other = self._validity_check(other)
     
-        result = self.clone_without_data([numpy.size(self.data, 0), 
+        result = self.clone_without_data([numpy.size(self.data, 0),
                                           numpy.size(self.data, 1)])
         
@@ -917,5 +919,2 @@
                 result.err_data[i][j] = math.sqrt(math.fabs(output.variance))
         return result
-    
-  
- 
Index: sansdataloader/src/sans/dataloader/loader.py
===================================================================
--- sansdataloader/src/sans/dataloader/loader.py	(revision 371cb85b0c40eb940e7b6c183b34791575781567)
+++ sansdataloader/src/sans/dataloader/loader.py	(revision f60a8c20e2f0018d19feff6fc41e188d2939a7e4)
@@ -1,27 +1,25 @@
-
+"""
+    File handler to support different file extensions.
+    Uses reflectometry's registry utility.
+    
+    The default readers are found in the 'readers' sub-module
+    and registered by default at initialization time.
+    
+    To add a new default reader, one must register it in
+    the register_readers method found in readers/__init__.py.
+    
+    A utility method (find_plugins) is available to inspect
+    a directory (for instance, a user plug-in directory) and
+    look for new readers/writers.
+"""
 #####################################################################
 #This software was developed by the University of Tennessee as part of the
 #Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
-#project funded by the US National Science Foundation. 
+#project funded by the US National Science Foundation.
 #See the license text in license.txt
 #copyright 2008, University of Tennessee
 ######################################################################
 
-"""
-    File handler to support different file extensions.
-    Uses reflectometry's registry utility.
-    
-    The default readers are found in the 'readers' sub-module
-    and registered by default at initialization time.
-    
-    To add a new default reader, one must register it in
-    the register_readers method found in readers/__init__.py. 
-    
-    A utility method (find_plugins) is available to inspect 
-    a directory (for instance, a user plug-in directory) and
-    look for new readers/writers.
-"""
-
-import os 
+import os
 import sys
 import logging
@@ -55,5 +53,5 @@
         readers.read_associations(self)
         
-        #TODO: remove the following line when ready to switch to 
+        #TODO: remove the following line when ready to switch to
         #the new default readers
         #readers.register_readers(self._identify_plugin)
@@ -62,5 +60,4 @@
         #self.find_plugins('plugins')
 
-        
     def load(self, path, format=None):
         """
@@ -68,10 +65,10 @@
 
         :param path: file path
-        :param format: explicit extension, to force the use 
+        :param format: explicit extension, to force the use
             of a particular reader
 
         Defaults to the ascii (multi-column) reader
         if no reader was registered for the file's
-        extension.   
+        extension.
         """
         try:
@@ -107,5 +104,5 @@
         dir = temp_path
         # Check whether the directory exists
-        if not os.path.isdir(dir): 
+        if not os.path.isdir(dir):
             msg = "DataLoader couldn't locate DataLoader plugin folder."
             msg += """ "%s" does not exist""" % dir
@@ -157,9 +154,9 @@
                         logging.error(msg)
                      
-        return readers_found 
+        return readers_found
     
     def associate_file_type(self, ext, module):
         """
-        Look into a module to find whether it contains a 
+        Look into a module to find whether it contains a
         Reader class. If so, APPEND it to readers and (potentially)
         to the list of writers for the given extension
@@ -226,5 +223,5 @@
                 type_name = loader.type_name
                 
-                wcard = "%s files (*%s)|*%s" % (type_name, ext.lower(), 
+                wcard = "%s files (*%s)|*%s" % (type_name, ext.lower(),
                                                 ext.lower())
                 if wcard not in self.wildcards:
@@ -237,5 +234,4 @@
         return reader_found
 
-    
     def _identify_plugin(self, module):
         """
@@ -268,5 +264,5 @@
                                                      ext.lower())
                     if wcard not in self.wildcards:
-                    	self.wildcards.append(wcard)
+                        self.wildcards.append(wcard)
                             
                 # Check whether writing is supported
@@ -275,5 +271,5 @@
                         if ext not in self.writers:
                             self.writers[ext] = []
-                        self.writers[ext].insert(0,loader.write)
+                        self.writers[ext].insert(0, loader.write)
                         
             except:
@@ -288,5 +284,5 @@
         
         :Raises ValueError: if file type is not known.
-        """        
+        """
         # Find matching extensions
         extlist = [ext for ext in self.extensions() if path.endswith(ext)]
@@ -301,5 +297,6 @@
             result = []
             for L in writers:
-                if L not in result: result.append(L)
+                if L not in result:
+                    result.append(L)
             writers = L
         # Raise an error if there are no matching extensions
@@ -316,5 +313,5 @@
         Raises KeyError if format is not available.
         
-        May raise a writer-defined exception if writer fails.        
+        May raise a writer-defined exception if writer fails.
         """
         if format is None:
@@ -326,7 +323,7 @@
                 return fn(path, data)
             except:
-                pass # give other loaders a chance to succeed
+                pass  # give other loaders a chance to succeed
         # If we get here it is because all loaders failed
-        raise # reraises last exception
+        raise  # reraises last exception
 
         
@@ -340,5 +337,5 @@
     def associate_file_type(self, ext, module):
         """
-        Look into a module to find whether it contains a 
+        Look into a module to find whether it contains a
         Reader class. If so, append it to readers and (potentially)
         to the list of writers for the given extension
@@ -394,16 +391,2 @@
     def get_wildcards(self):
         return self.__registry.wildcards
-        
-if __name__ == "__main__": 
-    logging.basicConfig(level=logging.INFO,
-                        format='%(asctime)s %(levelname)s %(message)s',
-                        filename='loader.log',
-                        filemode='w')
-    l = Loader()
-    test_data = l.load('test/cansas1d.xml')
-    l.save('test_file.xml', test_data, '.xml')
-    
-    print l.get_wildcards()
-        
-        
-    
Index: sansdataloader/src/sans/dataloader/manipulations.py
===================================================================
--- sansdataloader/src/sans/dataloader/manipulations.py	(revision 47045e660992f773f4d578ea70535571ce4090cd)
+++ sansdataloader/src/sans/dataloader/manipulations.py	(revision f60a8c20e2f0018d19feff6fc41e188d2939a7e4)
@@ -1,16 +1,15 @@
-
+"""
+Data manipulations for 2D data sets.
+Using the meta data information, various types of averaging
+are performed in Q-space
+"""
 #####################################################################
 #This software was developed by the University of Tennessee as part of the
 #Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
-#project funded by the US National Science Foundation. 
+#project funded by the US National Science Foundation.
 #See the license text in license.txt
 #copyright 2008, University of Tennessee
 ######################################################################
 
-"""
-Data manipulations for 2D data sets.
-Using the meta data information, various types of averaging
-are performed in Q-space 
-"""
 #TODO: copy the meta data from the 2D object to the resulting 1D object
 import math
@@ -33,4 +32,5 @@
     theta = 0.5 * math.atan(plane_dist/det_dist)
     return (4.0 * math.pi/wavelength) * math.sin(theta)
+
 
 def get_q_compo(dx, dy, det_dist, wavelength, compo=None):
@@ -55,4 +55,5 @@
     return out
 
+
 def flip_phi(phi):
     """
@@ -63,10 +64,11 @@
     Pi = math.pi
     if phi < 0:
-        phi_out = phi  + (2 * Pi)
+        phi_out = phi + (2 * Pi)
     elif phi > (2 * Pi):
-        phi_out = phi  - (2 * Pi)
+        phi_out = phi - (2 * Pi)
     else:
-        phi_out = phi 
+        phi_out = phi
     return phi_out
+
 
 def reader2D_converter(data2d=None):
@@ -93,10 +95,11 @@
     qy_data = new_y.flatten()
     q_data = numpy.sqrt(qx_data*qx_data + qy_data*qy_data)
-    if data2d.err_data == None or numpy.any(data2d.err_data <= 0): 
+    if data2d.err_data == None or numpy.any(data2d.err_data <= 0):
         new_err_data = numpy.sqrt(numpy.abs(new_data))
     else:
         new_err_data = data2d.err_data.flatten()
-    mask    = numpy.ones(len(new_data), dtype=bool)
-
+    mask = numpy.ones(len(new_data), dtype=bool)
+
+    #TODO: make sense of the following two lines...
     output = Data2D()
     output = data2d
@@ -109,4 +112,5 @@
 
     return output
+
 
 class _Slab(object):
@@ -149,9 +153,9 @@
             raise RuntimeError, msg
         
-        # Get data 
+        # Get data
         data = data2D.data[numpy.isfinite(data2D.data)]
         q_data = data2D.q_data[numpy.isfinite(data2D.data)]
         err_data = data2D.err_data[numpy.isfinite(data2D.data)]
-        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
+        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)]
         qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
              
@@ -159,24 +163,27 @@
         if maj == 'x':
             if self.fold:
-                x_min = 0     
-            else:  x_min = self.x_min
-            nbins = int(math.ceil((self.x_max - x_min)/self.bin_width))
+                x_min = 0
+            else:
+                x_min = self.x_min
+            nbins = int(math.ceil((self.x_max - x_min) / self.bin_width))
             qbins = self.bin_width * numpy.arange(nbins) + x_min
         elif maj == 'y':
-            if self.fold: y_min = 0     
-            else:  y_min = self.y_min
+            if self.fold:
+                y_min = 0
+            else:
+                y_min = self.y_min
             nbins = int(math.ceil((self.y_max - y_min)/self.bin_width))
-            qbins = self.bin_width * numpy.arange(nbins) + y_min          
+            qbins = self.bin_width * numpy.arange(nbins) + y_min
         else:
             raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj)
                                 
-        x  = numpy.zeros(nbins)
-        y  = numpy.zeros(nbins)
+        x = numpy.zeros(nbins)
+        y = numpy.zeros(nbins)
         err_y = numpy.zeros(nbins)
         y_counts = numpy.zeros(nbins)
 
-        # Average pixelsize in q space                                                
-        for npts in range(len(data)): 
-            # default frac                                               
+        # Average pixelsize in q space
+        for npts in range(len(data)):
+            # default frac
             frac_x = 0
             frac_y = 0
@@ -191,14 +198,14 @@
                 continue
             # binning: find axis of q
-            if maj == 'x': 
+            if maj == 'x':
                 q_value = qx_data[npts]
-                min = x_min 
-            if maj == 'y': 
-                q_value = qy_data[npts] 
+                min = x_min
+            if maj == 'y':
+                q_value = qy_data[npts]
                 min = y_min 
             if self.fold and q_value < 0:
-                q_value = -q_value   
+                q_value = -q_value
             # bin
-            i_q = int(math.ceil((q_value - min)/self.bin_width)) - 1
+            i_q = int(math.ceil((q_value - min) / self.bin_width)) - 1
             
             # skip outside of max bins
@@ -206,29 +213,27 @@
                 continue
             
-            # give it full weight
-            #frac = 1 
-            
             #TODO: find better definition of x[i_q] based on q_data
-            x[i_q] += frac * q_value#min + (i_q + 1) * self.bin_width / 2.0
+            x[i_q] += frac * q_value  # min + (i_q + 1) * self.bin_width / 2.0
             y[i_q] += frac * data[npts]
             
             if err_data == None or err_data[npts] == 0.0:
-                if data[npts] < 0: data[npts] = -data[npts]
+                if data[npts] < 0:
+                    data[npts] = -data[npts]
                 err_y[i_q] += frac * frac * data[npts]
             else:
                 err_y[i_q] += frac * frac * err_data[npts] * err_data[npts]
-            y_counts[i_q]  += frac
+            y_counts[i_q] += frac
                 
-        # Average the sums       
+        # Average the sums
         for n in range(nbins):
             err_y[n] = math.sqrt(err_y[n])
           
         err_y = err_y / y_counts
-        y    = y / y_counts
-        x    = x / y_counts
-        idx = (numpy.isfinite(y) & numpy.isfinite(x)) 
+        y = y / y_counts
+        x = x / y_counts
+        idx = (numpy.isfinite(y) & numpy.isfinite(x))
         
         if not idx.any(): 
-            msg = "Average Error: No points inside ROI to average..." 
+            msg = "Average Error: No points inside ROI to average..."
             raise ValueError, msg
         #elif len(y[idx])!= nbins:
@@ -237,4 +242,5 @@
         return Data1D(x=x[idx], y=y[idx], dy=err_y[idx])
         
+        
 class SlabY(_Slab):
     """
@@ -251,4 +257,5 @@
         return self._avg(data2D, 'y')
         
+        
 class SlabX(_Slab):
     """
@@ -264,6 +271,7 @@
         
         """
-        return self._avg(data2D, 'x') 
-        
+        return self._avg(data2D, 'x')
+
+
 class Boxsum(object):
     """
@@ -282,5 +290,5 @@
     def __call__(self, data2D):
         """
-        Perform the sum in the region of interest 
+        Perform the sum in the region of interest
          
         :param data2D: Data2D object
@@ -293,5 +301,5 @@
         # Average the sums
         counts = 0 if y_counts == 0 else y
-        error  = 0 if y_counts == 0 else math.sqrt(err_y)
+        error = 0 if y_counts == 0 else math.sqrt(err_y)
         
         return counts, error
@@ -299,9 +307,9 @@
     def _sum(self, data2D):
         """
-        Perform the sum in the region of interest 
-        
-        :param data2D: Data2D object
-        
-        :return: number of counts, 
+        Perform the sum in the region of interest
+        
+        :param data2D: Data2D object
+        
+        :return: number of counts,
             error on number of counts, number of entries summed
         
@@ -311,20 +319,20 @@
             msg += "of detectors: %g" % len(data2D.detector)
             raise RuntimeError, msg
-        # Get data 
+        # Get data
         data = data2D.data[numpy.isfinite(data2D.data)]
         q_data = data2D.q_data[numpy.isfinite(data2D.data)]
         err_data = data2D.err_data[numpy.isfinite(data2D.data)]
-        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
+        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)]
         qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
    
-        y  = 0.0
+        y = 0.0
         err_y = 0.0
         y_counts = 0.0
 
-        # Average pixelsize in q space                                                
-        for npts in range(len(data)): 
-            # default frac                                               
-            frac_x = 0 
-            frac_y = 0                 
+        # Average pixelsize in q space
+        for npts in range(len(data)):
+            # default frac
+            frac_x = 0
+            frac_y = 0
             
             # get min and max at each points
@@ -337,5 +345,5 @@
             if self.y_min <= qy and self.y_max > qy:
                 frac_y = 1
-            #Find the fraction along each directions            
+            #Find the fraction along each directions
             frac = frac_x * frac_y
             if frac == 0:
@@ -348,9 +356,8 @@
             else:
                 err_y += frac * frac * err_data[npts] * err_data[npts]
-            y_counts += frac     
+            y_counts += frac
         return y, err_y, y_counts
 
 
-      
 class Boxavg(Boxsum):
     """
@@ -363,5 +370,5 @@
     def __call__(self, data2D):
         """
-        Perform the sum in the region of interest 
+        Perform the sum in the region of interest
          
         :param data2D: Data2D object
@@ -373,12 +380,13 @@
         
         # Average the sums
-        counts = 0 if y_counts == 0 else y/y_counts
-        error  = 0 if y_counts == 0 else math.sqrt(err_y)/y_counts
+        counts = 0 if y_counts == 0 else y / y_counts
+        error = 0 if y_counts == 0 else math.sqrt(err_y) / y_counts
         
         return counts, error
         
+        
 def get_pixel_fraction_square(x, xmin, xmax):
     """
-    Return the fraction of the length 
+    Return the fraction of the length
     from xmin to x.::
    
@@ -437,15 +445,15 @@
         # Get the dq for resolution averaging
         if data2D.dqx_data != None and data2D.dqy_data != None:
-            # The pinholes and det. pix contribution present 
+            # The pinholes and det. pix contribution present
             # in both direction of the 2D which must be subtracted when
             # converting to 1D: dq_overlap should calculated ideally at
-            # q = 0. Note This method works on only pinhole geometry. 
+            # q = 0. Note This method works on only pinhole geometry.
             # Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average.
             z_max = max(data2D.q_data)
             z_min = min(data2D.q_data)
             x_max = data2D.dqx_data[data2D.q_data[z_max]]
-            x_min = data2D.dqx_data[data2D.q_data[z_min]]    
+            x_min = data2D.dqx_data[data2D.q_data[z_min]]
             y_max = data2D.dqy_data[data2D.q_data[z_max]]
-            y_min = data2D.dqy_data[data2D.q_data[z_min]]          
+            y_min = data2D.dqy_data[data2D.q_data[z_min]]
             # Find qdx at q = 0
             dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min)
@@ -453,5 +461,5 @@
             if dq_overlap_x > min(data2D.dqx_data):
                 dq_overlap_x = min(data2D.dqx_data)
-            dq_overlap_x *= dq_overlap_x 
+            dq_overlap_x *= dq_overlap_x
             # Find qdx at q = 0
             dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min)
@@ -462,5 +470,5 @@
             dq_overlap_y *= dq_overlap_y
 
-            dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y)/2.0)
+            dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0)
             # Final protection of dq
             if dq_overlap < 0:
@@ -470,5 +478,5 @@
             # def; dqx_data = dq_r dqy_data = dq_phi
             # Convert dq 2D to 1D here
-            dqx = dqx_data * dqx_data            
+            dqx = dqx_data * dqx_data
             dqy = dqy_data * dqy_data
             dq_data = numpy.add(dqx, dqy)
@@ -484,34 +492,34 @@
         qbins = self.bin_width * numpy.arange(nbins) + self.r_min
 
-        x  = numpy.zeros(nbins)
-        y  = numpy.zeros(nbins)
+        x = numpy.zeros(nbins)
+        y = numpy.zeros(nbins)
         err_y = numpy.zeros(nbins)
         err_x = numpy.zeros(nbins)
         y_counts = numpy.zeros(nbins)
 
-        for npt in range(len(data)):   
+        for npt in range(len(data)):
             
             if ismask and not mask_data[npt]:
-                continue   
+                continue
             
             frac = 0
             
             # q-value at the pixel (j,i)
-            q_value = q_data[npt] 
-            data_n = data[npt]                
+            q_value = q_data[npt]
+            data_n = data[npt]
             
             ## No need to calculate the frac when all data are within range
             if self.r_min >= self.r_max:
-                raise ValueError, "Limit Error: min > max" 
+                raise ValueError, "Limit Error: min > max"
             
             if self.r_min <= q_value and q_value <= self.r_max:
-                frac = 1                       
+                frac = 1
             if frac == 0:
                 continue 
-            i_q = int(math.floor((q_value - self.r_min) / self.bin_width))    
-
-            # Take care of the edge case at phi = 2pi. 
-            if i_q == nbins:  
-                i_q = nbins -1              
+            i_q = int(math.floor((q_value - self.r_min) / self.bin_width))
+
+            # Take care of the edge case at phi = 2pi.
+            if i_q == nbins:
+                i_q = nbins - 1
             y[i_q] += frac * data_n
             # Take dqs from data to get the q_average
@@ -524,16 +532,17 @@
                 err_y[i_q] += frac * frac * err_data[npt] * err_data[npt]
             if dq_data != None:
-                # To be consistent with dq calculation in 1d reduction, 
-                # we need just the averages (not quadratures) because 
-                # it should not depend on the number of the q points 
+                # To be consistent with dq calculation in 1d reduction,
+                # we need just the averages (not quadratures) because
+                # it should not depend on the number of the q points
                 # in the qr bins.
                 err_x[i_q] += frac * dq_data[npt]
             else:
                 err_x = None
-            y_counts[i_q]  += frac
-        
-        # Average the sums       
+            y_counts[i_q] += frac
+        
+        # Average the sums
         for n in range(nbins):
-            if err_y[n] < 0: err_y[n] = -err_y[n]
+            if err_y[n] < 0:
+                err_y[n] = -err_y[n]
             err_y[n] = math.sqrt(err_y[n])
             #if err_x != None:
@@ -541,8 +550,8 @@
             
         err_y = err_y / y_counts
-        err_y[err_y==0] = numpy.average(err_y)
-        y    = y / y_counts
-        x    = x / y_counts
-        idx = (numpy.isfinite(y)) & (numpy.isfinite(x)) 
+        err_y[err_y == 0] = numpy.average(err_y)
+        y = y / y_counts
+        x = x / y_counts
+        idx = (numpy.isfinite(y)) & (numpy.isfinite(x))
         
         if err_x != None:
@@ -551,6 +560,6 @@
             d_x = None
 
-        if not idx.any(): 
-            msg = "Average Error: No points inside ROI to average..." 
+        if not idx.any():
+            msg = "Average Error: No points inside ROI to average..."
             raise ValueError, msg
         
@@ -567,5 +576,5 @@
     around the ring as a function of phi.
     
-    Phi_min and phi_max should be defined between 0 and 2*pi 
+    Phi_min and phi_max should be defined between 0 and 2*pi
     in anti-clockwise starting from the x- axis on the left-hand side
     """
@@ -601,5 +610,5 @@
         q_data = data2D.q_data[numpy.isfinite(data2D.data)]
         err_data = data2D.err_data[numpy.isfinite(data2D.data)]
-        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
+        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)]
         qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
         
@@ -612,9 +621,9 @@
         phi_err    = numpy.zeros(self.nbins_phi)
         
-        for npt in range(len(data)):          
+        for npt in range(len(data)):
             frac = 0
             # q-value at the point (npt)
             q_value = q_data[npt]
-            data_n = data[npt]    
+            data_n = data[npt]
                         
             # phi-value at the point (npt)
@@ -622,5 +631,5 @@
             
             if self.r_min <= q_value and q_value <= self.r_max:
-                frac = 1                       
+                frac = 1
             if frac == 0:
                 continue
@@ -628,7 +637,7 @@
             i_phi = int(math.floor((self.nbins_phi) * phi_value / (2 * Pi)))
             
-            # Take care of the edge case at phi = 2pi. 
-            if i_phi == self.nbins_phi:  
-                i_phi =  self.nbins_phi - 1                            
+            # Take care of the edge case at phi = 2pi.
+            if i_phi == self.nbins_phi:
+                i_phi =  self.nbins_phi - 1
             phi_bins[i_phi] += frac * data[npt]
             
@@ -646,8 +655,8 @@
             phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i + 0.5)
             
-        idx = (numpy.isfinite(phi_bins)) 
+        idx = (numpy.isfinite(phi_bins))
 
         if not idx.any():
-            msg = "Average Error: No points inside ROI to average..." 
+            msg = "Average Error: No points inside ROI to average..."
             raise ValueError, msg
         #elif len(phi_bins[idx])!= self.nbins_phi:
@@ -656,8 +665,9 @@
         return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx])
     
+    
 def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11):
     """
     Returns the fraction of the pixel defined by
-    the four corners (q_00, q_01, q_10, q_11) that 
+    the four corners (q_00, q_01, q_10, q_11) that
     has q < qmax.::
     
@@ -714,8 +724,9 @@
     # this pixel and the constant-q ring. We only need
     # to know if we have to include it or exclude it.
-    elif (q_00 + q_01 + q_10 + q_11)/4.0 < qmax:
+    elif (q_00 + q_01 + q_10 + q_11) / 4.0 <  qmax:
         frac_max = 1.0
 
     return frac_max
+             
              
 def get_intercept(q, q_0, q_1):
@@ -727,7 +738,7 @@
     
     
-            A           B    
+            A           B
         +-----------+--------+    <--- pixel size
-        0                    1     
+        0                    1
         Q_0 -------- Q ----- Q_1   <--- equivalent Q range
         if Q_1 > Q_0, A is returned
@@ -738,20 +749,21 @@
     if q_1 > q_0:
         if (q > q_0 and q <= q_1):
-            return (q - q_0)/(q_1 - q_0)    
+            return (q - q_0) / (q_1 - q_0)
     else:
         if (q > q_1 and q <= q_0):
-            return (q - q_1)/(q_0 - q_1)
+            return (q - q_1) / (q_0 - q_1)
     return None
      
+     
 class _Sector:
     """
     Defines a sector region on a 2D data set.
     The sector is defined by r_min, r_max, phi_min, phi_max,
-    and the position of the center of the ring 
+    and the position of the center of the ring
     where phi_min and phi_max are defined by the right
     and left lines wrt central line
-    and phi_max could be less than phi_min. 
+    and phi_max could be less than phi_min.
    
-    Phi is defined between 0 and 2*pi in anti-clockwise 
+    Phi is defined between 0 and 2*pi in anti-clockwise
     starting from the x- axis on the left-hand side
     """
@@ -763,5 +775,4 @@
         self.nbins = nbins
         
-    
     def _agv(self, data2D, run='phi'):
         """
@@ -781,5 +792,5 @@
         q_data = data2D.q_data[numpy.isfinite(data2D.data)]
         err_data = data2D.err_data[numpy.isfinite(data2D.data)]
-        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
+        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)]
         qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
         dq_data = None
@@ -787,15 +798,15 @@
         # Get the dq for resolution averaging
         if data2D.dqx_data != None and data2D.dqy_data != None:
-            # The pinholes and det. pix contribution present 
+            # The pinholes and det. pix contribution present
             # in both direction of the 2D which must be subtracted when
             # converting to 1D: dq_overlap should calculated ideally at
-            # q = 0. 
+            # q = 0.
             # Extrapolate dqy(perp) at q = 0
             z_max = max(data2D.q_data)
             z_min = min(data2D.q_data)
             x_max = data2D.dqx_data[data2D.q_data[z_max]]
-            x_min = data2D.dqx_data[data2D.q_data[z_min]]    
+            x_min = data2D.dqx_data[data2D.q_data[z_min]]
             y_max = data2D.dqy_data[data2D.q_data[z_max]]
-            y_min = data2D.dqy_data[data2D.q_data[z_min]]          
+            y_min = data2D.dqy_data[data2D.q_data[z_min]]
             # Find qdx at q = 0
             dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min)
@@ -803,5 +814,5 @@
             if dq_overlap_x > min(data2D.dqx_data):
                 dq_overlap_x = min(data2D.dqx_data)
-            dq_overlap_x *= dq_overlap_x 
+            dq_overlap_x *= dq_overlap_x
             # Find qdx at q = 0
             dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min)
@@ -812,5 +823,5 @@
             dq_overlap_y *= dq_overlap_y
 
-            dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0)  
+            dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0)
             if dq_overlap < 0:
                 dq_overlap = y_min
@@ -819,5 +830,5 @@
             # def; dqx_data = dq_r dqy_data = dq_phi
             # Convert dq 2D to 1D here
-            dqx = dqx_data * dqx_data            
+            dqx = dqx_data * dqx_data
             dqy = dqy_data * dqy_data
             dq_data = numpy.add(dqx, dqy)
@@ -827,6 +838,6 @@
         x        = numpy.zeros(self.nbins)
         y        = numpy.zeros(self.nbins)
-        y_err    = numpy.zeros(self.nbins)   
-        x_err    = numpy.zeros(self.nbins)      
+        y_err    = numpy.zeros(self.nbins)
+        x_err    = numpy.zeros(self.nbins)
         y_counts = numpy.zeros(self.nbins)
                       
@@ -837,5 +848,5 @@
         q_data_max = numpy.max(q_data)
                       
-        for n in range(len(data)):     
+        for n in range(len(data)):
             frac = 0
             
@@ -848,14 +859,14 @@
             
             # phi-value of the pixel (j,i)
-            phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi 
+            phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi
             
             ## No need to calculate the frac when all data are within range
             if self.r_min <= q_value and q_value <= self.r_max:
-                frac = 1                       
+                frac = 1
             if frac == 0:
                 continue
             #In case of two ROIs (symmetric major and minor regions)(for 'q2')
-            if run.lower()=='q2':
-                ## For minor sector wing 
+            if run.lower() == 'q2':
+                ## For minor sector wing
                 # Calculate the minor wing phis
                 phi_min_minor = flip_phi(phi_min - Pi)
@@ -869,9 +880,9 @@
                               phi_value < phi_max_minor)
 
-            #For all cases(i.e.,for 'q', 'q2', and 'phi') 
-            #Find pixels within ROI  
-            if phi_min > phi_max: 
+            #For all cases(i.e.,for 'q', 'q2', and 'phi')
+            #Find pixels within ROI
+            if phi_min > phi_max:
                 is_in = is_in or (phi_value > phi_min or \
-                                   phi_value < phi_max)         
+                                   phi_value < phi_max)
             else:
                 is_in = is_in or (phi_value >= phi_min  and \
@@ -879,9 +890,9 @@
             
             if not is_in:
-                frac = 0                                                        
+                frac = 0
             if frac == 0:
                 continue
             # Check which type of averaging we need
-            if run.lower() == 'phi': 
+            if run.lower() == 'phi':
                 temp_x = (self.nbins) * (phi_value - self.phi_min)
                 temp_y = (self.phi_max - self.phi_min)
@@ -892,9 +903,9 @@
                 i_bin = int(math.floor(temp_x / temp_y))
 
-            # Take care of the edge case at phi = 2pi. 
-            if i_bin == self.nbins:  
-                i_bin =  self.nbins - 1
+            # Take care of the edge case at phi = 2pi.
+            if i_bin == self.nbins:
+                i_bin = self.nbins - 1
                 
-            ## Get the total y          
+            ## Get the total y
             y[i_bin] += frac * data_n
             x[i_bin] += frac * q_value
@@ -907,7 +918,7 @@
                 
             if dq_data != None:
-                # To be consistent with dq calculation in 1d reduction, 
-                # we need just the averages (not quadratures) because 
-                # it should not depend on the number of the q points 
+                # To be consistent with dq calculation in 1d reduction,
+                # we need just the averages (not quadratures) because
+                # it should not depend on the number of the q points
                 # in the qr bins.
                 x_err[i_bin] += frac * dq_data[n]
@@ -923,9 +934,9 @@
             # The type of averaging: phi,q2, or q
             # Calculate x[i]should be at the center of the bin
-            if run.lower()=='phi':               
+            if run.lower() == 'phi':
                 x[i] = (self.phi_max - self.phi_min) / self.nbins * \
                     (1.0 * i + 0.5) + self.phi_min
             else:
-                # We take the center of ring area, not radius.  
+                # We take the center of ring area, not radius.
                 # This is more accurate than taking the radial center of ring.
                 #delta_r = (self.r_max - self.r_min) / self.nbins
@@ -934,5 +945,5 @@
                 #x[i] = math.sqrt((r_inner * r_inner + r_outer * r_outer) / 2)
                 x[i] = x[i] / y_counts[i]
-        y_err[y_err==0] = numpy.average(y_err)        
+        y_err[y_err == 0] = numpy.average(y_err)
         idx = (numpy.isfinite(y) & numpy.isfinite(y_err))
         if x_err != None:
@@ -941,5 +952,5 @@
             d_x = None
         if not idx.any():
-            msg = "Average Error: No points inside sector of ROI to average..." 
+            msg = "Average Error: No points inside sector of ROI to average..."
             raise ValueError, msg
         #elif len(y[idx])!= self.nbins:
@@ -948,4 +959,5 @@
         return Data1D(x=x[idx], y=y[idx], dy=y_err[idx], dx=d_x)
                 
+                
 class SectorPhi(_Sector):
     """
@@ -963,6 +975,6 @@
         :return: Data1D object
         """
-
         return self._agv(data2D, 'phi')
+    
     
 class SectorQ(_Sector):
@@ -972,5 +984,5 @@
     
     A sector is defined by r_min, r_max, phi_min, phi_max.
-    r_min, r_max, phi_min, phi_max >0.  
+    r_min, r_max, phi_min, phi_max >0.
     The number of bin in Q also has to be defined.
     """
@@ -984,4 +996,5 @@
         """
         return self._agv(data2D, 'q2')
+
 
 class Ringcut(object):
@@ -993,8 +1006,8 @@
     The data returned is the region inside the ring
     
-    Phi_min and phi_max should be defined between 0 and 2*pi 
+    Phi_min and phi_max should be defined between 0 and 2*pi
     in anti-clockwise starting from the x- axis on the left-hand side
     """
-    def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0 ):
+    def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0):
         # Minimum radius
         self.r_min = r_min
@@ -1006,5 +1019,4 @@
         self.center_y = center_y
 
-        
     def __call__(self, data2D):
         """
@@ -1020,9 +1032,8 @@
 
         # Get data
-        qx_data = data2D.qx_data 
+        qx_data = data2D.qx_data
         qy_data = data2D.qy_data
         mask = data2D.mask
         q_data = numpy.sqrt(qx_data * qx_data + qy_data * qy_data)
-        #q_data_max = numpy.max(q_data)
 
         # check whether or not the data point is inside ROI
@@ -1051,5 +1062,5 @@
          
        :param data2D: Data2D object
-       :return: mask, 1d array (len = len(data)) 
+       :return: mask, 1d array (len = len(data))
            with Trues where the data points are inside ROI, otherwise False
         """
@@ -1060,14 +1071,14 @@
     def _find(self, data2D):
         """
-        Find a rectangular 2D region of interest. 
-        
-        :param data2D: Data2D object
-        
-        :return: out, 1d array (length = len(data)) 
+        Find a rectangular 2D region of interest.
+        
+        :param data2D: Data2D object
+        
+        :return: out, 1d array (length = len(data))
            with Trues where the data points are inside ROI, otherwise Falses
         """
         if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
             raise RuntimeError, "Boxcut take only plottable_2D objects"
-        # Get qx_ and qy_data 
+        # Get qx_ and qy_data
         qx_data = data2D.qx_data
         qy_data = data2D.qy_data
@@ -1080,12 +1091,13 @@
         return (outx & outy)
 
+
 class Sectorcut(object):
     """
     Defines a sector (major + minor) region on a 2D data set.
     The sector is defined by phi_min, phi_max,
-    where phi_min and phi_max are defined by the right 
-    and left lines wrt central line. 
+    where phi_min and phi_max are defined by the right
+    and left lines wrt central line.
    
-    Phi_min and phi_max are given in units of radian 
+    Phi_min and phi_max are given in units of radian
     and (phi_max-phi_min) should not be larger than pi
     """
@@ -1100,5 +1112,5 @@
         :param data2D: Data2D object
         
-        :return: mask, 1d array (len = len(data)) 
+        :return: mask, 1d array (len = len(data))
         
         with Trues where the data points are inside ROI, otherwise False
@@ -1110,18 +1122,18 @@
     def _find(self, data2D):
         """
-        Find a rectangular 2D region of interest. 
-        
-        :param data2D: Data2D object
-        
-        :return: out, 1d array (length = len(data)) 
+        Find a rectangular 2D region of interest.
+        
+        :param data2D: Data2D object
+        
+        :return: out, 1d array (length = len(data))
         
         with Trues where the data points are inside ROI, otherwise Falses
         """
         if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
-            raise RuntimeError, "Sectorcut take only plottable_2D objects" 
+            raise RuntimeError, "Sectorcut take only plottable_2D objects"
         Pi = math.pi
         # Get data 
         qx_data = data2D.qx_data
-        qy_data = data2D.qy_data        
+        qy_data = data2D.qy_data
         phi_data = numpy.zeros(len(qx_data))
 
@@ -1131,5 +1143,5 @@
         # Get the min and max into the region: -pi <= phi < Pi
         phi_min_major = flip_phi(self.phi_min + Pi) - Pi
-        phi_max_major = flip_phi(self.phi_max + Pi) - Pi  
+        phi_max_major = flip_phi(self.phi_max + Pi) - Pi
         # check for major sector
         if phi_min_major > phi_max_major:
@@ -1146,30 +1158,9 @@
         if phi_min_minor > phi_max_minor:
             out_minor = (phi_min_minor <= phi_data) + \
-                            (phi_max_minor >= phi_data) 
+                            (phi_max_minor >= phi_data)
         else:
             out_minor = (phi_min_minor <= phi_data) & \
-                            (phi_max_minor >= phi_data) 
+                            (phi_max_minor >= phi_data)
         out = out_major + out_minor
         
         return out
-
-if __name__ == "__main__": 
-
-    from loader import Loader
-    
-
-    d = Loader().load('test/MAR07232_rest.ASC')
-    #d = Loader().load('test/MP_New.sans')
-
-    
-    r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)
-    o = r(d)
-    
-    s = Ring(r_min=.000001, r_max=.01) 
-    p = s(d)
-    
-    for i in range(len(o.x)):
-        print o.x[i], o.y[i], o.dy[i]
-    
- 
-    
