Index: sansmodels/test/prototypes/test_2dsim.py
===================================================================
--- sansmodels/test/prototypes/test_2dsim.py	(revision 18e250cb66e3283c50a57dddaf692b17263dd72c)
+++ 	(revision )
@@ -1,82 +1,0 @@
-from sans_extension.prototypes.c_models import CCanvas
-from sans.models.SphereModel import SphereModel
-from sans.models.CoreShellModel import CoreShellModel
-
-import math
-
-# Sphere ###############################
-sphere = SphereModel()
-sphere.setParam('radius',40.0)
-sphere.setParam('contrast',1.0)
-sphere.setParam('scale',1.0)
-sphere.setParam('background',0.0)
-
-can = CCanvas()
-
-print "ID", can.add()
-can.setParam(0,1,40.0)
- 
-
-for q in [1,0.1,0.05,0.02,0.01,0.001,0.0001]:
-    ana = sphere.run([q,1.0])
-    sim = can.evaluate([q,1.0])
-
-    print "Sphere: ANA = %g   SIM = %g   (%g)" %(ana, sim, sim/ana)
-
-
-
-
-# Core-shell ###########################
-# This is a core shell model where the weight-average
-# of the SLDs is zero. The I(q) should go to zero as q->0.
-
-import time
-t_0 = time.time()
-
-sphere = CoreShellModel()
-# Core radius
-radius = 10.0
-thickness = 5.0
-sphere.setParam('radius', radius)
-# Shell thickness
-sphere.setParam('thickness', thickness)
-sphere.setParam('core_sld', 1.0)
-
-core_vol = 4.0/3.0*math.pi*radius*radius*radius
-outer_radius = radius+thickness
-shell_vol = 4.0/3.0*math.pi*outer_radius*outer_radius*outer_radius - core_vol
-shell_sld = -1.0*core_vol/shell_vol
-print "Shell SLD", shell_sld
-
-sphere.setParam('shell_sld', shell_sld)
-sphere.setParam('solvent_sld',0.0)
-sphere.setParam('background',0.0)
-sphere.setParam('scale',1.0)
-
-can = CCanvas()
-
-print "ID", can.add()
-can.setParam(0,1, outer_radius)
-can.setParam(0,2, shell_sld)
-
-
-print "ID", can.add()
-can.setParam(1,1, radius)
-can.setParam(1,2, 1.0)
-
-print "VOLS inner = %g, outer = %g, shell = %g" % (core_vol, shell_vol+core_vol, shell_vol)
-
-f = open('core-shell.txt','w')
-f.write("<q_cs>   <ans_cs>   <sim_cs>\n")
-for q in [1,0.1,0.05,0.04,0.03,0.02,0.01,0.005,0.001,0.0001]:
-#for i in range(400):
-    #q = 0.0001+i*0.005
-    ana = sphere.run([q,1.0])
-    sim = can.evaluate([q,1.0])
-    f.write("%g  %g  %g\n" %(q, ana, sim))
-    print "Core-shell: q=%g  ANA = %g   SIM = %g   (%g)" %(q, ana, sim, sim/ana)
-
-f.close()
-
-
-print "Time: ", time.time()-t_0
Index: sansmodels/test/prototypes/test_cyl_err2.py
===================================================================
--- sansmodels/test/prototypes/test_cyl_err2.py	(revision 18e250cb66e3283c50a57dddaf692b17263dd72c)
+++ 	(revision )
@@ -1,164 +1,0 @@
-try:
-    from sans.models.prototypes.SimCylinderF import SimCylinderF
-except:
-    print "This test uses the prototypes module."
-from sans.models.CylinderModel import CylinderModel
-import time, math, random
-
-def reset(seed=0):
-    sim = SimCylinderF()
-    sim.setParam('radius', 20)
-    sim.setParam('length', 200)
-    sim.setParam('phi', 1)
-    sim.setParam('theta', 1)
-    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.20
-
-
-counter = 0
-
-sim = reset()
-sph = CylinderModel()
-
-
-sph.setParam('radius', 20)
-sph.setParam('length', 200)
-sph.setParam('scale', 1)
-sph.setParam('contrast',1)
-sph.setParam('length', 200)
-sph.setParam('cyl_phi', 1)
-sph.setParam('cyl_theta', 1)
-
-vol = math.pi*20*20*200
-
-f = open('sim_err_cyl.txt', 'w')
-f.write("<q>  <npts>  <sigma>\n")
-
-
-
-
-for i_npts in range(1):
-
-    # 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
-    #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)
-                    
-        print 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()
-
Index: sansmodels/test/prototypes/test_cyl_iq.py
===================================================================
--- sansmodels/test/prototypes/test_cyl_iq.py	(revision 18e250cb66e3283c50a57dddaf692b17263dd72c)
+++ 	(revision )
@@ -1,77 +1,0 @@
-try:
-    from sans.models.prototypes.TestSphere2 import TestSphere2
-    from sans.models.prototypes.SimCylinderF import SimCylinderF
-    from sans.models.prototypes.SimCylinder import SimCylinder
-except:
-    print "This test uses the prototypes module."
-from sans.models.SphereModel import SphereModel
-from sans.models.CylinderModel import CylinderModel
-import math, time
-
-
-def test_lores():
-
-    lores = open('lores.iq', 'r')
-    f = open('lores_check.txt', 'w')
-    content = lores.read()
-    lines = content.split('\n')
-    inum = 0
-    norma = 0.
-    for line in lines:
-        toks = line.split()
-        if len(toks)==2:
-            inum += 1
-            q = float(toks[0])
-            ana = sph.run(q)
-            if inum==1:
-                continue
-            elif inum==2:
-                norma = ana/float(toks[1])
-            lor = float(toks[1])*norma
-            f.write("%g %g %g\n" % (q, ana, lor))
-    f.close()
-            
-
-#sim = TestSphere2()
-#sph = SphereModel()
-sim = SimCylinderF()
-sph = CylinderModel()
-
-sim.setParam('radius', 20)
-sim.setParam('length', 200)
-sim.setParam('phi', 1)
-sim.setParam('theta', 1)
-sim.setParam('npoints',2500)
-#print sph.params
-
-sph.setParam('radius', 20)
-sph.setParam('scale', 1)
-sph.setParam('contrast',1)
-sph.setParam('length', 200)
-sph.setParam('cyl_phi', 1)
-sph.setParam('cyl_theta', 1)
-
-#print "ANA", sph.run(0.1)
-#print "SIM", sim.run(0.1)
-
-f = open('sim_iq.txt', 'w')
-f.write("<q>  <ana>  <sim>\n")
-
-t_0 = time.time()
-
-for i in range(60):
-    q = 0.01 * (i+1) /3.0
-    #ana_value = sph.run([q,0])
-    #sim_value = sim.run([q,0])
-    ana_value = sph.run(q)
-    sim_value = sim.run(q)
-    ratio = 0
-    if sim_value>0:
-        ratio = sim_value/ana_value
-    print q, ana_value, sim_value, ratio
-    f.write("%10g %10g %10g\n" % (q, ana_value, sim_value))
-
-t_f = time.time()
-print "Time = ", t_f-t_0
-f.close()
-
Index: sansmodels/test/prototypes/test_simulation.py
===================================================================
--- sansmodels/test/prototypes/test_simulation.py	(revision 18e250cb66e3283c50a57dddaf692b17263dd72c)
+++ 	(revision )
@@ -1,275 +1,0 @@
-try:
-    from sans.models.prototypes.SimCylinder import SimCylinder
-except:
-    print "This test uses the prototypes module."
-    
-from sans.models.CylinderModel import CylinderModel
-
-import os, sys
-
-def simulate(filename="simul.txt", npts=500):
-    output_file = open(filename, 'w')
-    #output_file = open("simul_phi=0_theta=1_57_l=400_r=20_9.txt", 'w')
-    
-    sim = SimCylinder()
-    sim.setParam('scale', 1.0)
-    #sim.setParam('length', 100.0)
-    #sim.setParam('radius', 40.0)
-    sim.setParam('length', 400.0)
-    sim.setParam('radius', 20.0)
-    sim.setParam('theta', 1.57)
-    sim.setParam('phi', 0.0)
-    sim.setParam('qmax', 0.2)
-    
-    output_file.write("<q_value>  <ana_value>  <sim_value>\n")
-    norma = 0
-    for i in range(npts):
-        q = 1.0/npts*(i+1)
-        
-        cyl_val = sim.run(q)
-        sim_val = sim.run([q, 0.0])
-        if i==0:
-            norma = cyl_val/sim_val
-        
-        #print q, cyl_val, sim_val, sim_val/cyl_val
-        output_file.write("%10g %10g %10g\n" % (q, cyl_val, sim_val*norma))
-        
-    output_file.close()
-     
-def create_model():
-    sim = SimCylinder()
-    sim.setParam('scale', 1.0)
-    sim.setParam('length', 100.0)
-    sim.setParam('radius', 40.0)
-    # simul and simul_2
-    #sim.setParam('theta', 1.57)
-    #sim.setParam('phi', 0.0)
-    sim.setParam('theta', 0.0)
-    sim.setParam('phi', 1.0)
-    sim.setParam('qmax', 0.2)
-    return sim
-    
-def comp(): 
-    """
-        Simple test that should give 1 for the ratio
-        of simulated and analytical
-    """
-    sim = create_model()
-    
-    ana_val = sim.run(.05)
-    sim_val = sim.run([0.05, 0])
-    print ana_val, sim_val, ana_val/sim_val
-    
-    
-#    ana_val = sim.run(.05)
-#    sim_val = sim.run([0.05, 1])
-#    print ana_val, sim_val, ana_val/sim_val
-    
-#    sim = create_model()
-#    sim.setParam('theta', 1.0)
-#    sim.setParam('phi', 1.0)   
-#    ana_val = sim.run(.05)
-#    sim_val = sim.run([0.05, 1])
-#    print ana_val, sim_val, ana_val/sim_val
-    
-#    sim = create_model()
-#    sim.setParam('theta', 1.0)
-#    sim.setParam('phi', 2.0)
-#    ana_val = sim.run(.05)
-#    sim_val = sim.run([0.05, 1])
-#    print ana_val, sim_val, ana_val/sim_val
-    
-    sim = create_model()
-    sim.setParam('theta', 1.0)
-    sim.setParam('phi', -2.0)
-    ana_val = sim.run(.05)
-    sim_val = sim.run([0.15, 1])
-    print ana_val, sim_val, ana_val/sim_val
-    
-def add_error_to_file(prefix = "", dir = "output"):
-    import os, math
-    
-    # Open standard deviation file
-    std_file = open("%s/std_simul.txt" % dir, 'r')
-    std_content = std_file.read()
-    std_lines = std_content.split('\n')
-    
-    std_values = {}
-    for line in std_lines:
-        if line == std_lines[0]:
-            continue
-        
-        toks = line.split()
-        if not len(toks) > 3:
-            continue
-
-        std_values[toks[0]] = float(toks[2])
-    
-    std_file.close()
-    
-    # Read the files
-    
-    file_list = os.listdir(dir)
-    
-    for file in file_list:
-        if file.count(prefix)>0:
-            print "Reading", file
-    
-            file_obj = open("%s/%s" % (dir, file), 'r')
-            
-            # open new file
-            file_new = open("%s/error-%s" % (dir, file), 'w')
-            file_content = file_obj.read()
-            lines = file_content.split('\n')
-            for line in lines:
-                if line == lines[0]:
-                    file_new.write(line+'\n')
-                    continue
-                    
-                toks = line.split()
-                
-                if not len(toks) == 3:
-                    file_new.write(line+'\n')
-                    continue
-            
-                error = 0
-                if toks[0].lower() in std_values:
-                    error = float(toks[2]) * std_values[toks[0].lower()]
-                    
-                file_new.write(line+" %10g\n" % error)
-                
-            file_obj.close()
-            file_new.close()
-    
-    
-    
-def std_estimate(prefix = "", dir="output"):
-    import os, math
-    
-    file_list = os.listdir(dir)
-    output_file = open("%s/std_simul.txt" % dir, 'w')
-    output_file.write("<q_value>   <std>  <frac>\n")
-    values = {}
-    
-    for file in file_list:
-        if file.count(prefix)>0:
-            print "Reading", file
-            file_obj = open("%s/%s" % (dir, file), 'r')
-            file_content = file_obj.read()
-            lines = file_content.split('\n')
-            for line in lines:
-                if line == lines[0]:
-                    continue
-                toks = line.split()
-                
-                if not len(toks) == 3:
-                    continue
-                
-                if toks[0] not in values:
-                    values[toks[0]] = []
-                    
-                values[toks[0]].append(float(toks[2]))
-                                       
-    q_list = values.keys()
-    q_list.sort()                 
-    for q in q_list:
-        num = len(values[q])
-        sum = 0
-        mean = 0
-        for val in range(num):
-            mean += values[q][val] 
-        mean /= num
-            
-        for val in range(num):
-            diff = values[q][val] - mean
-            sum += diff * diff
-        sum /= (num-1)
-        #print q, math.sqrt(sum)
-    
-        output_file.write("%10s %10g %10g %10g\n" % (q, math.sqrt(sum), math.sqrt(sum)/mean, math.sqrt(sum)/mean*100.0))
-        
-    output_file.close()
-    
-    
-    # sqrt( 1/(N-1) sum[ (x-mean)**2 ] )
-    
-def test_volume_pts():
-    import random, math, time
-    
-    print "Generating points"
-    pt_list = []
-    for i in range(50000):
-        x = random.random()
-        y = random.random()
-        z = random.random()
-        pt_list.append([x, y, z])
-        
-    # Compute list of closest distance for each point
-    print "Computing distances"
-    closest_list = []
-    t_init = time.time()
-    
-    for i in range(len(pt_list)):
-        if math.fmod(i,100)==0:
-            print i, time.time()-t_init
-            t_init = time.time()
-        pt = pt_list[i]
-        
-        min_dist = -1.0
-        for j in range(i+1, len(pt_list)):
-            pt_2 = pt_list[j]
-            #print pt_2
-            dx = pt[0]-pt_2[0]
-            dy = pt[1]-pt_2[1]
-            dz = pt[2]-pt_2[2]
-            dist = math.sqrt( dx*dx + dy*dy + dz*dz ) 
-            if min_dist<0 or dist < min_dist:
-                min_dist = dist
-             
-        closest_list.append(min_dist)        
-        
-    # Analyze list of distances
-    print "Analyzing"
-    sum = 0
-    for d in closest_list:
-        sum += d
-    average = sum/len(closest_list)
-    
-    sum = 0
-    for d in closest_list:
-        diff = d - average
-        sum += diff*diff
-    sum /= len(closest_list)
-    
-    print "%5g +/- %5g)" % (average, math.sqrt(sum))
-    
-                              
-            
-    
-if __name__ == '__main__':
-    #test_volume_pts()
-    #simulate("simul_phi=0_theta=1_57_l=400_r=20_test.txt")
-    #std_estimate("simul_phi=0_theta=1_57_l=400_r=20")
-    #add_error_to_file("l=400_r=20.txt")
-    
-    option = "error"
-    
-    #if len(sys.argv)>1:
-    if len(option)>0:
-        if option=="gen":
-            ifile = 0
-            for i in range(100):
-                ifile += 1
-                #print "simul_qmax=1_phi=0_theta=1_57_l=400_r=20_%g.txt" % ifile
-                #simulate("simul_phi=0_theta=1_57_l=400_r=20_%g.txt" % ifile)
-                print "simul_qmax=1_phi=0_theta=1_57_l=400_r=20_%g.txt" % ifile
-                simulate("output/pt5000/simul_phi=0_theta=1_57_l=400_r=20_%g.txt" % ifile)
-        elif option=="comp":
-            comp()
-        elif option=="std":
-            #std_estimate("simul_phi=0_theta=1_57_l=400_r=20")
-            std_estimate("simul_phi=0_theta=1_57_l=400_r=20", "output/pt5000")
-        elif option=="error":
-            #add_error_to_file("l=400_r=20.txt")
-            add_error_to_file("l=400_r=20_1.txt", "output/pt5000")
-            
Index: sansmodels/test/prototypes/test_sphere_err2.py
===================================================================
--- sansmodels/test/prototypes/test_sphere_err2.py	(revision 18e250cb66e3283c50a57dddaf692b17263dd72c)
+++ 	(revision )
@@ -1,156 +1,0 @@
-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()
-
Index: sansmodels/test/prototypes/test_sphere_iq.py
===================================================================
--- sansmodels/test/prototypes/test_sphere_iq.py	(revision 18e250cb66e3283c50a57dddaf692b17263dd72c)
+++ 	(revision )
@@ -1,43 +1,0 @@
-try:
-    from sans.models.prototypes.SimSphereF import SimSphereF
-except:
-    print "This test uses the prototypes module."
-from sans.models.SphereModel import SphereModel
-import math, time
-
-
-
-sim = SimSphereF()
-sph = SphereModel()
-
-sim.setParam('radius', 60)
-sim.setParam('npoints', 50000)
-
-sph.setParam('radius', 60)
-sph.setParam('scale', 1)
-
-
-f = open('sim_iq.txt', 'w')
-f.write("<q>  <ana>  <sim> <err>\n")
-
-t_0 = time.time()
-
-#for i in range(60):
-for i in range(100):
-    #q = 0.01 * (i+1) /3.0
-    q = 0.05 + 0.20*(i+1) /100
-    
-    ana_value = sph.run([q,0])
-    sim_value = sim.run([q,0])
-#    ana_value = sph.run(q)
-#    sim_value = sim.run(q)
-    ratio = 0
-    if sim_value>0:
-        ratio = sim_value/ana_value
-    print q, ana_value, sim_value, ratio
-    f.write("%10g %10g %10g %10g\n" % (q, ana_value, sim_value, (sim_value-ana_value)/ana_value))
-
-t_f = time.time()
-print "Time = ", t_f-t_0
-f.close()
-
