Changeset 3190 for production/sydney_2006
- Timestamp:
- Jun 21, 2006, 9:31:57 AM (19 years ago)
- Location:
- production/sydney_2006
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
production/sydney_2006/export_results.py
r2875 r3190 4 4 from pyvolution.ermapper_grids import read_ermapper_grid 5 5 6 name = project.outputname 26 name = project.outputname 7 7 8 8 #print 'Which variable do you want to export?' -
production/sydney_2006/find_max.py
r2878 r3190 35 35 c = 0 36 36 direction = 'negative' 37 xstar = x0 38 print 'start of loop', deriv 37 39 while c < 100000 and deriv > 0: 38 40 deriv = dfuncdx(x,0.83) 39 if deriv < 0: xstar = x 41 if deriv < 0: 42 xstar = x 43 print 'hello', deriv, xstar/lam0, c 40 44 if direction == 'positive': x += step 41 45 if direction == 'negative': x -= step 42 46 c += 1 47 43 48 44 49 print 'location of maximum of surface elevation function: xstar = %f' % (xstar/lam0) 45 50 const = 1.0 #a3D = ? #sydney = 86? and kappad=1 46 51 47 x = arange(-20,25,0.001)48 y = arange(-10,10,0.1)49 X,Y = meshgrid(x,y)52 #x = arange(-20,25,0.001) 53 #y = arange(-10,10,0.1) 54 #X,Y = meshgrid(x,y) 50 55 51 test = 0.052 for xi in x:53 testi = func(xi,0.0,0.83)54 if direction == 'positive':55 if testi > test:56 test = testi57 xstar = xi58 else:59 if testi < test:60 test = testi61 xstar = xi56 #test = 0.0 57 #for xi in x: 58 # testi = func(xi,0.0,0.83) 59 # if direction == 'positive': 60 # if testi > test: 61 # test = testi 62 # xstar = xi 63 # else: 64 # if testi < test: 65 # test = testi 66 # xstar = xi 62 67 63 print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)68 #print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test) -
production/sydney_2006/get_timeseries.py
r2885 r3190 5 5 import project 6 6 from pyvolution.util import file_function 7 from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn7 #from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn 8 8 from pylab import * 9 9 from matplotlib.ticker import MultipleLocator, FormatStrFormatter 10 10 11 swwfile = project.outputname + '.sww' 12 13 11 swwfile = project.outputname + '.sww' 12 #swwfile = project.outputdir + timestampdir + sep + project.outputname + '.sww' 13 # place to store figures 14 #graphloc = project.outputdir + timestampdir + sep 15 graphloc = project.outputdir 14 16 15 17 #Time interval to plot 16 tmin = 1300017 tmax = 2100018 tmin = 2*60 19 tmax = 10*60 18 20 19 21 def get_gauges_from_file(filename): … … 34 36 gaugelocation.append(location) 35 37 36 #Return gauges and raw data for subsequent storage37 38 return gauges, lines, gaugelocation 38 39 … … 40 41 gauges, lines, locations = get_gauges_from_file(project.gauge_filename) 41 42 42 print 'number of gauges for Benfield: ', len(gauges)43 print 'number of gauges: ', len(gauges) 43 44 44 45 #Read model output … … 50 51 use_cache = True) 51 52 52 53 53 print 'size f', size(f.quantities['stage'],axis=0), size(f.quantities['stage'],axis=1) 54 54 55 from math import sqrt, atan, degrees 55 56 from Numeric import ones … … 85 86 uh = f(t, point_id = k)[2] 86 87 vh = f(t, point_id = k)[3] 87 myloc = locations[k]88 gaugeloc = locations[k] 88 89 depth = w-z 89 90 90 91 m = sqrt(uh*uh + vh*vh) #Absolute momentum 91 92 vel = sqrt(uh*uh + vh*vh) / (w-z + 1.e-30) #Absolute velocity … … 135 136 model_time, elevations, '-k') 136 137 #name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1]) 137 name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], myloc)138 name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], gaugeloc) 138 139 title(name) 139 140 … … 144 145 shadow=True, 145 146 loc='upper right') 146 #savefig('Gauge_%d_stage' %k) 147 savefig('Gauge_%s_stage' %myloc) 148 149 # raw_input('Next') 147 #('Gauge_%d_stage' %k) 148 savefig('%sGauge_%s_stage' %(graphloc, gaugeloc)) 149 #savefig('Gauge_%s_stage.eps' %gaugeloc) 150 150 151 151 #Momentum plot … … 159 159 ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]') 160 160 #savefig('Gauge_%d_momentum' %k) 161 savefig(' Gauge_%s_momentum' %myloc)161 savefig('%sGauge_%s_momentum' %(graphloc, gaugeloc)) 162 162 163 # raw_input('Next')164 165 163 #Bearing plot 166 164 ion() … … 181 179 ylabel(' atan(vh/uh) [degrees from North]') 182 180 #savefig('Gauge_%d_bearing' %k) 183 savefig('Gauge_%s_bearing' %myloc) 184 185 # raw_input('Next') 181 savefig('%sGauge_%s_bearing' %(graphloc, gaugeloc)) 186 182 187 183 #Speed plot … … 195 191 ylabel('sqrt( uh^2 + vh^2 ) / depth [m/s]') 196 192 #savefig('Gauge_%d_speed' %k) 197 savefig('Gauge_%s_speed' %myloc) 198 199 # raw_input('Next') 200 201 whichone = '_%s' %myloc 193 savefig('%sGauge_%s_speed' %(graphloc, gaugeloc)) 194 195 whichone = '_%s' %gaugeloc 202 196 thisfile = project.gaugetimeseries+whichone+'.csv' 203 197 fid = open(thisfile, 'w') -
production/sydney_2006/plot_integral.py
r2487 r3190 30 30 hold(False) 31 31 plot(t, integral, '.-b') 32 axis([0, max(t), min(integral)-0.1, min(integral) + 1]) 32 33 title('Checking for water conservation') 33 34 xlabel('Time (sec)') -
production/sydney_2006/process_gauges.py
r2885 r3190 5 5 import project 6 6 from pyvolution.util import file_function 7 from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn7 #from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn 8 8 from pylab import * 9 9 … … 11 11 plot = False 12 12 13 swwfile = project.outputname + '.sww'13 swwfile = project.outputname4 + '.sww' 14 14 15 15 def get_gauges_from_file(filename): -
production/sydney_2006/project.py
r2875 r3190 30 30 nmaxviz = 6283000 31 31 32 basename = 'slump_ 315res'33 basename2 = 'slump_1000res '32 basename = 'slump_duncan' 33 basename2 = 'slump_1000res_1min' 34 34 basename4 = 'slump_poly_ingo_test' 35 35 … … 62 62 63 63 gauge_filename = gaugedir + 'sydney_gauges.xya' 64 buildings_filename = gaugedir + 'all_bld_ind.csv' 64 65 #gauge_filename = gaugedir + 'sydney_gauges_test.xya' 65 66 gauge_outname = gaugedir + 'gauges_max_output.xya' … … 137 138 138 139 # clipping used for demo purposes to fit around poly3 and poly4 (Ingo's files) 139 j0 = [318000, 6249000] 140 j1 = [387000, 6249000] 141 j2 = [387000, 6270000] 142 j3 = [318000, 6270000] 143 144 demopoly = [j0, j1, j2, j3] 140 #j0 = [318000, 6249000] 141 #j1 = [387000, 6249000] 142 #j2 = [387000, 6270000] 143 #j3 = [318000, 6270000] 144 145 #demopoly = [j0, j1, j2, j3] 146 147 # second demo poly which isn't rectangular 148 j0 = [385000, 6280000] 149 j1 = [360000, 6272500] 150 j2 = [335000, 6272500] 151 j3 = [330000, 6265000] 152 j31 = [325000, 6260000] 153 j4 = [316000, 6260000] 154 j5 = [316000, 6247000] 155 j6 = [350000, 6247000] 156 j7 = [385000, 6238000] 157 158 demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7] 145 159 146 160 # to put chunk back in -
production/sydney_2006/run_sydney_smf.py
r2640 r3190 27 27 from pyvolution.quantity import Quantity 28 28 from Numeric import allclose 29 from utilities.polygon import inside_polygon 29 30 30 31 # Application specific imports … … 70 71 71 72 # original issue to Benfield 72 #interior_res = 500073 #interior_regions = [[project.harbour_polygon_2, interior_res],74 #[project.botanybay_polygon_2, interior_res]]73 interior_res = 5000 74 interior_regions = [[project.harbour_polygon_2, interior_res], 75 [project.botanybay_polygon_2, interior_res]] 75 76 76 77 # used for finer mesh 77 interior_res1 = 5000 78 interior_res2 = 315 79 interior_regions = [[project.newpoly1, interior_res1], 80 [project.south1, interior_res1], 81 [project.finepolymanly, interior_res2], 82 [project.finepolyquay, interior_res2]] 83 84 print 'number of interior regions', len(interior_regions) 78 #interior_res1 = 5000 79 #interior_res2 = 315 80 #interior_regions = [[project.newpoly1, interior_res1], 81 # [project.south1, interior_res1], 82 # [project.finepolymanly, interior_res2], 83 # [project.finepolyquay, interior_res2]] 84 85 # used for coastal polygons constructed by GIS guys 86 def get_polygon_from_file(filename): 87 """ Function to read in output from GIS determined polygon 88 """ 89 fid = open(filename) 90 lines = fid.readlines() 91 fid.close() 92 93 polygon = [] 94 for line in lines[1:]: 95 fields = line.split(',') 96 x = float(fields[1]) 97 y = float(fields[2]) 98 polygon.append([x, y]) 99 100 return polygon 101 102 num_polygons = 9 103 fileext = '.csv' 104 filename = project.polygonptsfile 105 106 #interior_res = 1000 107 #interior_regions = [] 108 #bounding_polygon = project.diffpolygonall#project.demopoly 109 #count = 0 110 #for p in range(1, num_polygons+1): 111 # thefilename = filename + str(p) + fileext 112 # print 'reading in polygon points', thefilename 113 # interior_polygon = get_polygon_from_file(thefilename) 114 # interior_regions.append([interior_polygon, interior_res]) 115 # n = len(interior_polygon) 116 # # check interior polygon falls in bounding polygon 117 # if len(inside_polygon(interior_polygon, bounding_polygon, 118 # closed = True, verbose = False)) <> len(interior_polygon): 119 # print 'WARNING: interior polygon %d is outside bounding polygon' %(p) 120 # count += 1 121 # # check for duplicate points in interior polygon 122 123 print 'number of interior polygons: ', len(interior_regions) 124 #if count == 0: print 'interior polygons OK' 85 125 86 126 #FIXME: Fix caching of this one once the mesh_interface is ready … … 89 129 # original + refined region 90 130 _ = cache(create_mesh_from_regions, 91 project.diffpolygonall, 131 #project.demopoly, 132 project.diffpolygonall2, 133 #{'boundary_tags': {'bottom': [0], 134 # 'right1': [1], 'right0': [2], 135 # 'right2': [3], 'top': [4], 'left1': [5], 136 # 'left2': [6], 'left3': [7]}, 92 137 {'boundary_tags': {'bottom': [0], 93 ' right1': [1], 'right0': [2],94 ' right2': [3], 'top': [4], 'left1': [5],138 'bottom1': [1], 'right': [2], 139 'top1': [3], 'top': [4], 'left1': [5], 95 140 'left2': [6], 'left3': [7]}, 141 #{'boundary_tags': {'bottom': [0], 'right': [1], 142 # 'top': [2], 'left': [3]}, 96 143 'maximum_triangle_area': 100000, 97 144 'filename': meshname, … … 126 173 radius=3330, 127 174 dphi=0.23, 128 x0=project.slump_origin [0],129 y0=project.slump_origin [1],175 x0=project.slump_origin2[0], 176 y0=project.slump_origin2[1], 130 177 alpha=0.0, 131 domain=domain )132 178 domain=domain, 179 verbose=True) 133 180 134 181 #------------------------------------------------------------------------------- … … 138 185 # apply slump after 30 mins, initialise to water level of tide = 0 139 186 domain.set_quantity('stage', 0.0) 140 domain.set_quantity('friction', 0.0 3)187 domain.set_quantity('friction', 0.04) 141 188 domain.set_quantity('elevation', 142 189 filename = project.combineddemname + '.pts', … … 158 205 # 'right2': Br, 'top': Br, 'left1': Br, 159 206 # 'left2': Br, 'left3': Br} ) 160 207 # for new tests 4 April 2006 208 #domain.set_boundary( {'bottom': Br, 'bottom1': Br, 'right': Br, 209 # 'top1': Br, 'top': Br, 'left1': Br, 210 # 'left2': Br, 'left3': Br} ) 211 161 212 # enforce Dirichlet BC - from 30/03/06 Benfield visit 162 domain.set_boundary( {'bottom': Bd, ' right1': Bd, 'right0': Bd,163 ' right2': Bd, 'top': Bd, 'left1': Bd,213 domain.set_boundary( {'bottom': Bd, 'bottom1': Bd, 'right': Bd, 214 'top1': Bd, 'top': Bd, 'left1': Bd, 164 215 'left2': Bd, 'left3': Bd} ) 165 216 166 217 # increasingly finer interior regions 167 #domain.set_boundary( {'bottom': B r, 'right': Br, 'left': Br, 'top': Br} )218 #domain.set_boundary( {'bottom': Bd, 'right': Bd, 'left': Bd, 'top': Bd} ) 168 219 169 220 … … 177 228 fid = open(thisfile, 'w') 178 229 179 for t in domain.evolve(yieldstep = 120, finaltime = 18000): 230 # save every 10 secs leading up to slump initiation 231 for t in domain.evolve(yieldstep = 10, finaltime = 60): # 6 steps 180 232 domain.write_time() 181 233 domain.write_boundary_statistics(tags = 'bottom') 182 183 # calculate integral 184 thisstagestep = domain.get_quantity('stage') 185 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 186 fid.write(s) 187 188 # add slump after 30 mins 189 if allclose(t, 30*60): 234 # calculate integral 235 thisstagestep = domain.get_quantity('stage') 236 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 237 fid.write(s) 238 # add slump after 30 secs 239 if allclose(t, 30): 190 240 slump = Quantity(domain) 191 241 slump.set_values(tsunami_source) 192 242 domain.set_quantity('stage', slump + thisstagestep) 193 243 #test_stage = domain.get_quantity('stage') 244 #test_elevation = domain.get_quantity('elevation') 245 #test_depth = test_stage - test_elevation 246 #test_max = max(test_depth.get_values()) 247 #print 'testing', test_max 248 249 import sys; sys.exit() 250 251 # save every two minutes leading up to interesting period 252 for t in domain.evolve(yieldstep = 120, finaltime = 660, # steps 253 skip_initial_step = True): 254 domain.write_time() 255 domain.write_boundary_statistics(tags = 'bottom') 256 # calculate integral 257 thisstagestep = domain.get_quantity('stage') 258 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 259 fid.write(s) 260 261 262 # save every thirty secs during interesting period 263 for t in domain.evolve(yieldstep = 60, finaltime = 5000, # steps 264 skip_initial_step = True): 265 domain.write_time() 266 domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage') 267 # calculate integral 268 thisstagestep = domain.get_quantity('stage') 269 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 270 fid.write(s) 271 272 273 import sys; sys.exit() 274 # save every two mins for next 5000 secs 275 for t in domain.evolve(yieldstep = 120, finaltime = 10000, # about 42 steps 276 skip_initial_step = True): 277 domain.write_time() 278 domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage') 279 # calculate integral 280 thisstagestep = domain.get_quantity('stage') 281 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 282 fid.write(s) 283 284 # save every half hour to end of simulation 285 for t in domain.evolve(yieldstep = 1800, finaltime = 10*60*60, # 14 steps 286 skip_initial_step = True): 287 domain.write_time() 288 domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage' 289 # calculate integral 290 thisstagestep = domain.get_quantity('stage') 291 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 292 fid.write(s) 194 293 195 294 print 'That took %.2f seconds' %(time.time()-t0) -
production/sydney_2006/run_sydney_smf_polyingo.py
r2732 r3190 20 20 21 21 # Related major packages 22 from pyvolution.shallow_water import Domain, Reflective_boundary, Time_boundary 22 from pyvolution.shallow_water import Domain, Reflective_boundary, Time_boundary, Dirichlet_boundary 23 23 from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts 24 24 from pyvolution.combine_pts import combine_rectangular_points_files … … 26 26 from pyvolution.quantity import Quantity 27 27 from geospatial_data import * 28 from utilities.polygon import inside_polygon 29 from Numeric import allclose 28 30 29 31 # Application specific imports … … 81 83 meshname = project.meshname4+'.msh' 82 84 83 interior_res1 = 2500085 interior_res1 = 1000 84 86 interior_res2 = 1315 85 print 'interior resolution', interior_res1, interior_res286 87 87 88 # Read in files containing polygon points (generated by GIS) … … 104 105 return polygon 105 106 106 num_polygons = 3107 num_polygons = 4 107 108 fileext = '.csv' 108 109 filename = project.polygonptsfile 109 110 110 111 interior_regions = [] 111 112 bounding_polygon = project.diffpolygonall 113 count = 0 112 114 for p in range(3, num_polygons+1): 113 115 thefilename = filename + str(p) + fileext … … 115 117 interior_polygon = get_polygon_from_file(thefilename) 116 118 interior_regions.append([interior_polygon, interior_res1]) 117 119 n = len(interior_polygon) 120 # check interior polygon falls in bounding polygon 121 if len(inside_polygon(interior_polygon, bounding_polygon, 122 closed = True, verbose = False)) <> len(interior_polygon): 123 print 'WARNING: interior polygon %d is outside bounding polygon' %(p) 124 count += 1 125 # check for duplicate points in interior polygon 126 118 127 print 'number of interior polygons: ', len(interior_regions) 119 120 #print 'Number of interior polygons: ', len(interior_regions) 128 if count == 0: print 'interior polygons OK' 121 129 122 130 # original … … 132 140 # [project.finepolyquay, interior_res2]] 133 141 142 134 143 #FIXME: Fix caching of this one once the mesh_interface is ready 135 144 from caching import cache 145 136 146 137 147 #_ = cache(create_mesh_from_regions, … … 141 151 # 'right2': [3], 'top': [4], 'left1': [5], 142 152 # 'left2': [6], 'left3': [7]}, 143 # 'maximum_triangle_area': 1 00000,153 # 'maximum_triangle_area': 150000, 144 154 # 'filename': meshname, 145 155 # 'interior_regions': interior_regions}, 146 156 # verbose = True) 147 157 158 # for demo construction 148 159 _ = cache(create_mesh_from_regions, 149 project.d iffpolygonall_test,160 project.demopoly, 150 161 {'boundary_tags': {'bottom': [0], 'right': [1], 151 162 'top': [2], 'left': [3]}, … … 171 182 domain.set_datadir(project.outputdir) 172 183 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 173 174 175 184 176 185 … … 190 199 domain=domain) 191 200 192 #print 'testing', tsunami_source.get_integral()193 194 201 #------------------------------------------------------------------------------- 195 202 # Setup initial conditions … … 199 206 domain.set_quantity('stage', 0) #tsunami_source) 200 207 domain.set_quantity('friction', 0.0) # supplied by Benfield, initial value 0.03 201 domain.set_quantity('elevation', alpha = 10.0,#G, alpha = 550.0, 202 #filename = project.combineddemname + '.pts', 208 domain.set_quantity('elevation', 203 209 filename = project.coarsedemname + '.pts', 204 use_cache = False,210 use_cache = True, 205 211 verbose = True) 206 212 … … 231 237 232 238 Br = Reflective_boundary(domain) 233 234 #domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br, 235 # 'right2': Br, 'top': Br, 'left1': Br, 236 # 'left2': Br, 'left3': Br} ) 237 domain.set_boundary( {'bottom': Br, 'right': Br, 'left': Br, 'top': Br} ) 239 Bd = Dirichlet_boundary([0, 0, 0]) 240 241 # for demo 242 domain.set_boundary( {'bottom': Bd, 'right': Bd, 'left': Bd, 'top': Bd} ) 238 243 #domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br, 239 244 # 'right2': Br, 'top': Br, 'left1': Br, … … 250 255 fid = open(thisfile, 'w') 251 256 252 # original 253 for t in domain.evolve(yieldstep = 1, finaltime = 10): 257 for t in domain.evolve(yieldstep = 10, finaltime = 600): 254 258 domain.write_time() 255 259 domain.write_boundary_statistics(tags = 'bottom') … … 257 261 s = '%.2f, %.2f\n' %(t, stagestep.get_integral()) 258 262 fid.write(s) 263 if allclose(t, 30): 264 slump = Quantity(domain) 265 slump.set_values(tsunami_source) 266 domain.set_quantity('stage', slump + stagestep) 267 268 # save every two minutes leading up to interesting period 269 for t in domain.evolve(yieldstep = 120, finaltime = 660, # steps 270 skip_initial_step = True): 271 domain.write_time() 272 domain.write_boundary_statistics(tags = 'bottom') 273 # calculate integral 274 thisstagestep = domain.get_quantity('stage') 275 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 276 fid.write(s) 277 278 # save every thirty secs during interesting period 279 for t in domain.evolve(yieldstep = 30, finaltime = 3500, # steps 280 skip_initial_step = True): 281 domain.write_time() 282 domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage') 283 # calculate integral 284 thisstagestep = domain.get_quantity('stage') 285 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 286 fid.write(s) 287 288 # save every two mins for next 5000 secs 289 for t in domain.evolve(yieldstep = 120, finaltime = 10000, # about 42 steps 290 skip_initial_step = True): 291 domain.write_time() 292 domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage') 293 # calculate integral 294 thisstagestep = domain.get_quantity('stage') 295 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 296 fid.write(s) 297 298 # save every half hour to end of simulation 299 for t in domain.evolve(yieldstep = 1800, finaltime = 10*60*60, # 14 steps 300 skip_initial_step = True): 301 domain.write_time() 302 domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage' 303 # calculate integral 304 thisstagestep = domain.get_quantity('stage') 305 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 306 fid.write(s) 259 307 260 308 print 'That took %.2f seconds' %(time.time()-t0) -
production/sydney_2006/run_timeseries.py
r2795 r3190 14 14 15 15 # False to generate latex output, True otherwise 16 title_on = False16 title_on = True 17 17 18 18 # generate figures - see sww2timeseries for documentation -
production/sydney_2006/smf_volume_conservation.py
r2673 r3190 7 7 from Numeric import ones 8 8 from pylab import * 9 import scipy10 from scipy.special import erf9 #import scipy 10 #from scipy.special import erf 11 11 12 12 ion() … … 39 39 eta2 = func(X,0,1) 40 40 41 #import sys; sys.exit() 41 42 #eta = func(X,Y,kappad) 42 43 #im1 = imshow(eta, cmap=cm.jet, interpolation='nearest')
Note: See TracChangeset
for help on using the changeset viewer.