Changeset 3190 for production/sydney_2006

Timestamp:

Jun 21, 2006, 9:31:57 AM (19 years ago)

Author:

sexton

Message:

MOST and ANUGA comparisons and updates

Location:

production/sydney_2006

Files:

• export_results.py (modified) (1 diff)
• find_max.py (modified) (1 diff)
• get_timeseries.py (modified) (10 diffs)
• plot_integral.py (modified) (1 diff)
• process_gauges.py (modified) (2 diffs)
• project.py (modified) (3 diffs)
• run_sydney_smf.py (modified) (7 diffs)
• run_sydney_smf_polyingo.py (modified) (13 diffs)
• run_timeseries.py (modified) (1 diff)
• smf_volume_conservation.py (modified) (2 diffs)

production/sydney_2006/export_results.py

@@ -4 +4 @@
 from pyvolution.ermapper_grids import read_ermapper_grid
 
-name = project.outputname2
+name = project.outputname
 
 #print 'Which variable do you want to export?'

production/sydney_2006/find_max.py

@@ -35 +35 @@
 c = 0
 direction = 'negative'
+xstar = x0
+print 'start of loop', deriv
 while c < 100000 and deriv > 0:
     deriv = dfuncdx(x,0.83)
-    if deriv < 0: xstar = x
+    if deriv < 0:
+        xstar = x
+        print 'hello', deriv, xstar/lam0, c
     if direction == 'positive': x += step
     if direction == 'negative': x -= step
     c += 1
+    
 
 print 'location of maximum of surface elevation function: xstar = %f' % (xstar/lam0)
 const = 1.0 #a3D = ? #sydney = 86? and kappad=1
 
-x = arange(-20,25,0.001)
-y = arange(-10,10,0.1)
-X,Y = meshgrid(x,y)
+#x = arange(-20,25,0.001)
+#y = arange(-10,10,0.1)
+#X,Y = meshgrid(x,y)
 
-test = 0.0
-for xi in x:
-    testi = func(xi,0.0,0.83)
-    if direction == 'positive':
-        if testi > test:
-            test = testi
-            xstar = xi
-    else:
-        if testi < test:
-            test = testi
-            xstar = xi
+#test = 0.0
+#for xi in x:
+#    testi = func(xi,0.0,0.83)
+#    if direction == 'positive':
+#        if testi > test:
+#            test = testi
+#            xstar = xi
+#    else:
+#        if testi < test:
+#            test = testi
+#            xstar = xi
     
-print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)
+#print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)

production/sydney_2006/get_timeseries.py

@@ -5 +5 @@
 import project
 from pyvolution.util import file_function
-from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
+#from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
 from pylab import *
 from matplotlib.ticker import MultipleLocator, FormatStrFormatter
 
-swwfile = project.outputname + '.sww' 
-
-
+swwfile = project.outputname + '.sww'
+#swwfile = project.outputdir + timestampdir + sep + project.outputname + '.sww'
+# place to store figures
+#graphloc = project.outputdir + timestampdir + sep
+graphloc = project.outputdir
 
 #Time interval to plot
-tmin = 13000
-tmax = 21000
+tmin = 2*60
+tmax = 10*60
 
 def get_gauges_from_file(filename):
@@ -34 +36 @@
         gaugelocation.append(location)
         
-    #Return gauges and raw data for subsequent storage
     return gauges, lines, gaugelocation
 
@@ -40 +41 @@
 gauges, lines, locations = get_gauges_from_file(project.gauge_filename)
 
-print 'number of gauges for Benfield:   ', len(gauges)
+print 'number of gauges:   ', len(gauges)
 
 #Read model output
@@ -50 +51 @@
                   use_cache = True)
 
-
-               
+print 'size f', size(f.quantities['stage'],axis=0), size(f.quantities['stage'],axis=1)
+
 from math import sqrt, atan, degrees
 from Numeric import ones
@@ -85 +86 @@
         uh = f(t, point_id = k)[2]
         vh = f(t, point_id = k)[3]
-        myloc = locations[k]
+        gaugeloc = locations[k]
         depth = w-z
-
+        
         m = sqrt(uh*uh + vh*vh)   #Absolute momentum
         vel = sqrt(uh*uh + vh*vh) / (w-z + 1.e-30) #Absolute velocity
@@ -135 +136 @@
              model_time, elevations, '-k')
     #name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
-    name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], myloc)
+    name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], gaugeloc)
     title(name)
 
@@ -144 +145 @@
            shadow=True,
            loc='upper right')
-    #savefig('Gauge_%d_stage' %k)
-    savefig('Gauge_%s_stage' %myloc)
-
-    # raw_input('Next')
+    #('Gauge_%d_stage' %k)
+    savefig('%sGauge_%s_stage' %(graphloc, gaugeloc))
+    #savefig('Gauge_%s_stage.eps' %gaugeloc)
     
     #Momentum plot
@@ -159 +159 @@
     ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')    
     #savefig('Gauge_%d_momentum' %k)
-    savefig('Gauge_%s_momentum' %myloc)
+    savefig('%sGauge_%s_momentum' %(graphloc, gaugeloc))
     
-    # raw_input('Next')
-
     #Bearing plot
     ion()
@@ -181 +179 @@
     ylabel(' atan(vh/uh) [degrees from North]')    
     #savefig('Gauge_%d_bearing' %k)
-    savefig('Gauge_%s_bearing' %myloc)
-    
-    # raw_input('Next')
+    savefig('%sGauge_%s_bearing' %(graphloc, gaugeloc))
 
     #Speed plot
@@ -195 +191 @@
     ylabel('sqrt( uh^2 + vh^2 ) / depth [m/s]')    
     #savefig('Gauge_%d_speed' %k)
-    savefig('Gauge_%s_speed' %myloc)
-    
-    # raw_input('Next')
-
-    whichone = '_%s' %myloc
+    savefig('%sGauge_%s_speed' %(graphloc, gaugeloc))
+
+    whichone = '_%s' %gaugeloc
     thisfile = project.gaugetimeseries+whichone+'.csv'
     fid = open(thisfile, 'w')

production/sydney_2006/plot_integral.py

@@ -30 +30 @@
 hold(False)
 plot(t, integral, '.-b')
+axis([0, max(t), min(integral)-0.1, min(integral) + 1])
 title('Checking for water conservation')
 xlabel('Time (sec)')

production/sydney_2006/process_gauges.py

@@ -5 +5 @@
 import project
 from pyvolution.util import file_function
-from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
+#from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
 from pylab import *
 
@@ -11 +11 @@
 plot = False
     
-swwfile = project.outputname + '.sww'
+swwfile = project.outputname4 + '.sww'
 
 def get_gauges_from_file(filename):

production/sydney_2006/project.py

@@ -30 +30 @@
 nmaxviz = 6283000
 
-basename = 'slump_315res'
-basename2 = 'slump_1000res'
+basename = 'slump_duncan'
+basename2 = 'slump_1000res_1min'
 basename4 = 'slump_poly_ingo_test'
 
@@ -62 +62 @@
 
 gauge_filename = gaugedir + 'sydney_gauges.xya'
+buildings_filename = gaugedir + 'all_bld_ind.csv'
 #gauge_filename = gaugedir + 'sydney_gauges_test.xya'
 gauge_outname = gaugedir + 'gauges_max_output.xya'
@@ -137 +138 @@
 
 # clipping used for demo purposes to fit around poly3 and poly4 (Ingo's files)
-j0 = [318000, 6249000]
-j1 = [387000, 6249000]
-j2 = [387000, 6270000]
-j3 = [318000, 6270000]
-
-demopoly = [j0, j1, j2, j3]
+#j0 = [318000, 6249000]
+#j1 = [387000, 6249000]
+#j2 = [387000, 6270000]
+#j3 = [318000, 6270000]
+
+#demopoly = [j0, j1, j2, j3]
+
+# second demo poly which isn't rectangular
+j0 = [385000, 6280000]
+j1 = [360000, 6272500]
+j2 = [335000, 6272500]
+j3 = [330000, 6265000]
+j31 = [325000, 6260000]
+j4 = [316000, 6260000]
+j5 = [316000, 6247000]
+j6 = [350000, 6247000]
+j7 = [385000, 6238000]
+
+demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7]
 
 # to put chunk back in

production/sydney_2006/run_sydney_smf.py

@@ -27 +27 @@
 from pyvolution.quantity import Quantity
 from Numeric import allclose
+from utilities.polygon import inside_polygon
 
 # Application specific imports
@@ -70 +71 @@
 
 # original issue to Benfield
-#interior_res = 5000
-#interior_regions = [[project.harbour_polygon_2, interior_res],
-#                    [project.botanybay_polygon_2, interior_res]]
+interior_res = 5000
+interior_regions = [[project.harbour_polygon_2, interior_res],
+                    [project.botanybay_polygon_2, interior_res]]
 
 # used for finer mesh
-interior_res1 = 5000
-interior_res2 = 315
-interior_regions = [[project.newpoly1, interior_res1],
-                    [project.south1, interior_res1],
-                    [project.finepolymanly, interior_res2],
-                    [project.finepolyquay, interior_res2]]
-
-print 'number of interior regions', len(interior_regions)
+#interior_res1 = 5000
+#interior_res2 = 315
+#interior_regions = [[project.newpoly1, interior_res1],
+#                    [project.south1, interior_res1],
+#                    [project.finepolymanly, interior_res2],
+#                    [project.finepolyquay, interior_res2]]
+
+# used for coastal polygons constructed by GIS guys
+def get_polygon_from_file(filename):
+    """ Function to read in output from GIS determined polygon
+    """
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+    
+    polygon = []
+    for line in lines[1:]:
+       fields = line.split(',')
+       x = float(fields[1])
+       y = float(fields[2])
+       polygon.append([x, y])
+        
+    return polygon
+
+num_polygons = 9
+fileext = '.csv'
+filename = project.polygonptsfile
+
+#interior_res = 1000
+#interior_regions = []
+#bounding_polygon = project.diffpolygonall#project.demopoly
+#count = 0
+#for p in range(1, num_polygons+1):
+#    thefilename = filename + str(p) + fileext
+#    print 'reading in polygon points', thefilename
+#    interior_polygon = get_polygon_from_file(thefilename)
+#    interior_regions.append([interior_polygon, interior_res])
+#    n = len(interior_polygon)
+#    # check interior polygon falls in bounding polygon
+#    if len(inside_polygon(interior_polygon, bounding_polygon,
+#                   closed = True, verbose = False)) <> len(interior_polygon):
+#        print 'WARNING: interior polygon %d is outside bounding polygon' %(p)
+#        count += 1
+#    # check for duplicate points in interior polygon
+    
+print 'number of interior polygons: ', len(interior_regions)
+#if count == 0: print 'interior polygons OK'
 
 #FIXME: Fix caching of this one once the mesh_interface is ready
@@ -89 +129 @@
 # original + refined region
 _ = cache(create_mesh_from_regions,
-          project.diffpolygonall,
+          #project.demopoly,
+          project.diffpolygonall2,
+          #{'boundary_tags': {'bottom': [0], 
+          #                   'right1': [1], 'right0': [2],
+          #                   'right2': [3], 'top': [4], 'left1': [5],
+          #                   'left2': [6], 'left3': [7]},
           {'boundary_tags': {'bottom': [0], 
-                             'right1': [1], 'right0': [2],
-                             'right2': [3], 'top': [4], 'left1': [5],
+                             'bottom1': [1], 'right': [2],
+                             'top1': [3], 'top': [4], 'left1': [5],
                              'left2': [6], 'left3': [7]},
+          #{'boundary_tags': {'bottom': [0], 'right': [1],
+          #                   'top': [2], 'left': [3]},
            'maximum_triangle_area': 100000,
            'filename': meshname,           
@@ -126 +173 @@
                                radius=3330,
                                dphi=0.23,
-                               x0=project.slump_origin[0], 
-                               y0=project.slump_origin[1], 
+                               x0=project.slump_origin2[0], 
+                               y0=project.slump_origin2[1], 
                                alpha=0.0, 
-                               domain=domain)
-
+                               domain=domain,
+                               verbose=True)
 
 #-------------------------------------------------------------------------------                                 
@@ -138 +185 @@
 # apply slump after 30 mins, initialise to water level of tide = 0
 domain.set_quantity('stage', 0.0)
-domain.set_quantity('friction', 0.03) 
+domain.set_quantity('friction', 0.04) 
 domain.set_quantity('elevation',
                     filename = project.combineddemname + '.pts',
@@ -158 +205 @@
 #                      'right2': Br, 'top': Br, 'left1': Br,
 #                      'left2': Br, 'left3': Br} )
-
+# for new tests 4 April 2006
+#domain.set_boundary( {'bottom': Br, 'bottom1': Br, 'right': Br,
+#                      'top1': Br, 'top': Br, 'left1': Br,
+#                      'left2': Br, 'left3': Br} )
+                             
 # enforce Dirichlet BC - from 30/03/06 Benfield visit
-domain.set_boundary( {'bottom': Bd, 'right1': Bd, 'right0': Bd,
-                      'right2': Bd, 'top': Bd, 'left1': Bd,
+domain.set_boundary( {'bottom': Bd, 'bottom1': Bd, 'right': Bd,
+                      'top1': Bd, 'top': Bd, 'left1': Bd,
                       'left2': Bd, 'left3': Bd} )
 
 # increasingly finer interior regions
-#domain.set_boundary( {'bottom': Br, 'right': Br, 'left': Br, 'top': Br} )
+#domain.set_boundary( {'bottom': Bd, 'right': Bd, 'left': Bd, 'top': Bd} )
 
 
@@ -177 +228 @@
 fid = open(thisfile, 'w')
 
-for t in domain.evolve(yieldstep = 120, finaltime = 18000):
+# save every 10 secs leading up to slump initiation
+for t in domain.evolve(yieldstep = 10, finaltime = 60): # 6 steps
     domain.write_time()
     domain.write_boundary_statistics(tags = 'bottom')
-
-    # calculate integral
-    thisstagestep = domain.get_quantity('stage') 
-    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
-    fid.write(s)
-
-    # add slump after 30 mins
-    if allclose(t, 30*60):
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+    # add slump after 30 secs
+    if allclose(t, 30):
         slump = Quantity(domain)
         slump.set_values(tsunami_source)
         domain.set_quantity('stage', slump + thisstagestep)
-    
+        #test_stage = domain.get_quantity('stage')
+        #test_elevation = domain.get_quantity('elevation')
+        #test_depth = test_stage - test_elevation
+        #test_max = max(test_depth.get_values())
+        #print 'testing', test_max
+
+import sys; sys.exit()
+
+# save every two minutes leading up to interesting period
+for t in domain.evolve(yieldstep = 120, finaltime = 660, # steps
+                       skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+
+
+# save every thirty secs during interesting period
+for t in domain.evolve(yieldstep = 60, finaltime = 5000, # steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+
+
+import sys; sys.exit()
+# save every two mins for next 5000 secs
+for t in domain.evolve(yieldstep = 120, finaltime = 10000, # about 42 steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+    
+# save every half hour to end of simulation    
+for t in domain.evolve(yieldstep = 1800, finaltime = 10*60*60, # 14 steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage'
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
     
 print 'That took %.2f seconds' %(time.time()-t0)

production/sydney_2006/run_sydney_smf_polyingo.py

@@ -20 +20 @@
 
 # Related major packages
-from pyvolution.shallow_water import Domain, Reflective_boundary, Time_boundary
+from pyvolution.shallow_water import Domain, Reflective_boundary, Time_boundary, Dirichlet_boundary
 from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
 from pyvolution.combine_pts import combine_rectangular_points_files 
@@ -26 +26 @@
 from pyvolution.quantity import Quantity
 from geospatial_data import *
+from utilities.polygon import inside_polygon
+from Numeric import allclose
 
 # Application specific imports
@@ -81 +83 @@
 meshname = project.meshname4+'.msh'
 
-interior_res1 = 25000
+interior_res1 = 1000
 interior_res2 = 1315
-print 'interior resolution', interior_res1, interior_res2
 
 # Read in files containing polygon points (generated by GIS)
@@ -104 +105 @@
     return polygon
 
-num_polygons = 3
+num_polygons = 4
 fileext = '.csv'
 filename = project.polygonptsfile
 
 interior_regions = []
-
+bounding_polygon = project.diffpolygonall
+count = 0
 for p in range(3, num_polygons+1):
     thefilename = filename + str(p) + fileext
@@ -115 +117 @@
     interior_polygon = get_polygon_from_file(thefilename)
     interior_regions.append([interior_polygon, interior_res1])
-
+    n = len(interior_polygon)
+    # check interior polygon falls in bounding polygon
+    if len(inside_polygon(interior_polygon, bounding_polygon,
+                   closed = True, verbose = False)) <> len(interior_polygon):
+        print 'WARNING: interior polygon %d is outside bounding polygon' %(p)
+        count += 1
+    # check for duplicate points in interior polygon
+    
 print 'number of interior polygons: ', len(interior_regions)
-
-#print 'Number of interior polygons: ', len(interior_regions)
+if count == 0: print 'interior polygons OK'
 
 # original
@@ -132 +140 @@
 #                    [project.finepolyquay, interior_res2]]
 
+
 #FIXME: Fix caching of this one once the mesh_interface is ready
 from caching import cache
+
 
 #_ = cache(create_mesh_from_regions,
@@ -141 +151 @@
 #                             'right2': [3], 'top': [4], 'left1': [5],
 #                             'left2': [6], 'left3': [7]},
-#           'maximum_triangle_area': 100000,
+#           'maximum_triangle_area': 150000,
 #           'filename': meshname,           
 #           'interior_regions': interior_regions},
 #          verbose = True)
 
+# for demo construction
 _ = cache(create_mesh_from_regions,
-          project.diffpolygonall_test,
+          project.demopoly,
           {'boundary_tags': {'bottom': [0], 'right': [1],
                              'top': [2], 'left': [3]},
@@ -171 +182 @@
 domain.set_datadir(project.outputdir)
 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-
-
 
 
@@ -190 +199 @@
                                domain=domain)
 
-#print 'testing', tsunami_source.get_integral()
-
 #-------------------------------------------------------------------------------                                 
 # Setup initial conditions
@@ -199 +206 @@
 domain.set_quantity('stage', 0) #tsunami_source)
 domain.set_quantity('friction', 0.0) # supplied by Benfield, initial value 0.03
-domain.set_quantity('elevation', alpha = 10.0,#G, alpha = 550.0, 
-                    #filename = project.combineddemname + '.pts',
+domain.set_quantity('elevation', 
                     filename = project.coarsedemname + '.pts',
-                    use_cache = False,
+                    use_cache = True,
                     verbose = True)
 
@@ -231 +237 @@
 
 Br = Reflective_boundary(domain)
-
-#domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
-#                      'right2': Br, 'top': Br, 'left1': Br,
-#                      'left2': Br, 'left3': Br} )
-domain.set_boundary( {'bottom': Br, 'right': Br, 'left': Br, 'top': Br} )
+Bd = Dirichlet_boundary([0, 0, 0])
+
+# for demo
+domain.set_boundary( {'bottom': Bd, 'right': Bd, 'left': Bd, 'top': Bd} )
 #domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
 #                      'right2': Br, 'top': Br, 'left1': Br,
@@ -250 +255 @@
 fid = open(thisfile, 'w')
     
-# original
-for t in domain.evolve(yieldstep = 1, finaltime = 10):
+for t in domain.evolve(yieldstep = 10, finaltime = 600):
     domain.write_time()
     domain.write_boundary_statistics(tags = 'bottom')
@@ -257 +261 @@
     s = '%.2f, %.2f\n' %(t, stagestep.get_integral())
     fid.write(s)
+    if allclose(t, 30):
+        slump = Quantity(domain)
+        slump.set_values(tsunami_source)
+        domain.set_quantity('stage', slump + stagestep)
+
+# save every two minutes leading up to interesting period
+for t in domain.evolve(yieldstep = 120, finaltime = 660, # steps
+                       skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+    
+# save every thirty secs during interesting period
+for t in domain.evolve(yieldstep = 30, finaltime = 3500, # steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+
+# save every two mins for next 5000 secs
+for t in domain.evolve(yieldstep = 120, finaltime = 10000, # about 42 steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+    
+# save every half hour to end of simulation    
+for t in domain.evolve(yieldstep = 1800, finaltime = 10*60*60, # 14 steps
+                       skip_initial_step = True):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage'
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
 
 print 'That took %.2f seconds' %(time.time()-t0)

production/sydney_2006/run_timeseries.py

@@ -14 +14 @@
 
 # False to generate latex output, True otherwise
-title_on = False
+title_on = True
 
 # generate figures - see sww2timeseries for documentation

production/sydney_2006/smf_volume_conservation.py

@@ -7 +7 @@
 from Numeric import ones
 from pylab import *
-import scipy
-from scipy.special import erf
+#import scipy
+#from scipy.special import erf
 
 ion()
@@ -39 +39 @@
 eta2 = func(X,0,1)
 
+#import sys; sys.exit()
 #eta = func(X,Y,kappad)
 #im1 = imshow(eta, cmap=cm.jet, interpolation='nearest')