Index: sansmodels/src/sans/models/test/prototypes/test_sphere_err2.py
===================================================================
--- sansmodels/src/sans/models/test/prototypes/test_sphere_err2.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/test/prototypes/test_sphere_err2.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,156 @@
+try:
+    from sans.models.prototypes.SimSphereF import SimSphereF
+except:
+    print "This test uses the prototypes module."
+from sans.models.SphereModel import SphereModel
+import time, math, random
+
+def reset(seed=0):
+    sim = SimSphereF()
+    sim.setParam('radius', 60)
+    sim.setParam('seed', math.fmod(time.time()*100+seed,10000))
+    return sim
+
+# number of times to repeat each point
+N_q = 50
+# Relative error on the simulation
+sig_val = 0.10
+
+
+counter = 0
+
+sim = reset()
+sph = SphereModel()
+
+
+sph.setParam('radius', 60)
+sph.setParam('scale', 1)
+
+vol = 4/3*math.pi*60*60*60
+
+f = open('sim_err_4.txt', 'w')
+f.write("<q>  <npts>  <sigma>\n")
+
+
+
+
+for i_npts in range(15):
+
+    # Find which q has an error of sig_val
+    #print err    
+    q_sigma = 0
+    i_min = 0
+    i_max = 0
+    i_last_max = 0
+    val_min = 0
+    val_max = 0
+    last_error = 0
+    last_val = 0
+    error_vec = []
+    
+    found_it = False
+    i_vec = 0
+        
+    
+    t_0 = time.time()
+    npts = 50000*i_npts+1000
+    #npts = 1000 * math.pow(2.0,(i_npts+1))
+    
+    sim = reset()
+    sim.setParam('npoints', npts)
+    
+    # Vectors of errors
+    err = []
+    for i in range(80):
+        q = 0.40*(i) /50
+        
+        
+        # For each Q, repeat the simulation N_q times
+        sum_q_ana = 0
+        sum_q_sim = 0
+        sum = 0
+        vec = []
+        for i_q in range(N_q):
+            counter += 1
+            sim = reset(counter)
+            sim.setParam('npoints', npts)
+                    
+            simval = sim.run([q, 0])
+            sum_q_ana += sph.run([q, 0])
+            sum_q_sim += simval
+            
+            vec.append(simval)
+            
+        mean = sum_q_sim/N_q
+        for val in range(N_q):
+            diff = vec[val] - mean
+            sum += diff * diff
+        sum /= (N_q-1)
+
+        ana_value = sum_q_ana/N_q
+        sim_value = sum_q_sim/N_q
+        
+        
+        
+        try:    
+            err_i = [q, (sim_value-ana_value)/ana_value, ana_value, math.sqrt(sum)/ana_value]
+        except:
+            err_i = [q,0,0,0]
+        
+        err.append(err_i)
+                    
+         
+    
+        
+        # update extrema
+        if val_min > err_i[2]:
+            val_min = err_i[2]
+            
+        if val_max < err_i[2]:
+            val_max = err_i[2]
+        
+        # Find a minimum, where the error changes sign
+        if last_val==val_min and err_i[2] > val_min:
+            
+            
+            #print "  min", err_i[0]
+            if i_min > 0:
+                i_plateau = i_last_max
+      
+                # check error for plateau
+                #if math.fabs(err[i_plateau][1])>sig_val:
+                if math.fabs(err[i_max][3])>sig_val:
+                    error_vec.append(err[i_plateau])
+                    t_f = time.time()
+                    pred = sph.run([err[i_plateau][0],0])
+                    print '[',t_f-t_0,']',math.pow((vol/npts),-.6666), 'n=',npts, 'q=',err[i_plateau][0], \
+                        'err=',math.fabs(err[i_plateau][1]), i_min, 'sig=',err[i_plateau][3]
+                        
+                    f.write("%10g %10g %10g\n" % (err[i_plateau][0], npts, err[i_plateau][3]))
+                    found_it = True
+                    break
+                
+            i_min = i_vec
+            # Get ready for next maximum
+            val_max = val_min
+            
+        # check if we are at max
+        elif last_val==val_max and err_i[2]<val_max:
+            #print "   max", err[i_vec-1][0], err[i_vec-1][1]
+            i_last_max = i_max
+            i_max = i_vec-1
+            # get ready for next minimum
+            val_min = val_max
+
+        last_error = err_i[1]
+        last_val = err_i[2]
+        i_vec += 1
+
+        
+    if found_it == False:
+        print "Could not complete", npts
+    #print "Npts = %g; %s" %(npts, str(len(error_vec)))
+        
+
+f.close()
+
