/[escript]/trunk/modellib/py_src/tsunami.py
ViewVC logotype

Diff of /trunk/modellib/py_src/tsunami.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 247 by gross, Tue Nov 29 04:57:38 2005 UTC revision 251 by gross, Tue Nov 29 05:54:06 2005 UTC
# Line 1  Line 1 
1  # $Id:$  # $Id$
2    
3    
4  import os  import os
# Line 48  class GridData: Line 48  class GridData:
48            i_dx=1            i_dx=1
49            i_dy=1            i_dy=1
50            for i in range(x_array.shape[0]):            for i in range(x_array.shape[0]):
51               x0=x_array[i,0]-ox               x_long=x_array[i,0]-ox
52               x1=x_array[i,1]-oy               x_lat=x_array[i,1]-oy
53               j0=min(max(int(x0/dx),0),data_array.shape[1]-1-i_dy)               j0=min(max(int(x_long/dx),0),data_array.shape[1]-1-i_dy)
54               j1=min(max(int(x1/dy),0),data_array.shape[0]-1-i_dx)               j1=min(max(int(x_lat/dy),0),data_array.shape[0]-1-i_dx)
55               f01=(x0-j0*dx)/dx               f01=(x_long-j0*dx)/dx
56               f00=1.-f01               f00=1.-f01
57               f11=(x1-j1*dy)/dy               f11=(x_lat-j1*dy)/dy
58               f10=1.-f11               f10=1.-f11
59               H00=data_array[j1,j0]               H00=data_array[j1,j0]
60               H01=data_array[j1,j0+i_dx]               H01=data_array[j1,j0+i_dx]
# Line 76  class PointOnEarthSurface: Line 76  class PointOnEarthSurface:
76     """     """
77     coordinates of a  point on the surface of the Earth     coordinates of a  point on the surface of the Earth
78     """     """
79     def __init__(self,phi=0,theta=0):     def __init__(self,long=0,lat=0):
80         self.phi=phi         self.long=long
81         self.theta=theta         self.lat=lat
82     def __str__(self):     def __str__(self):
83         return "(%s,%s)"%(self.phi,self.theta)         return "(%s,%s)"%(self.long,self.lat)
84    
85     def __sub__(self,other):     def __sub__(self,other):
86         return self.dist(other)         return self.dist(other)
87    
88     def split(self,p,t):     def split(self,p,t):
89        phi=self.phi+t*(p.phi-self.phi)        return PointOnEarthSurface(long=self.long+t*(p.long-self.long),lat=self.lat+t*(p.lat-self.lat))
       theta=self.theta+t*(p.theta-self.theta)  
       return PointOnEarthSurface(phi,theta)  
90    
91     def midPoint(self,other):     def midPoint(self,other):
92        return PointOnEarthSurface((self.phi+other.phi)/2,(self.theta+other.theta)/2)        return PointOnEarthSurface(long=(self.long+other.long)/2,lat=(self.lat+other.lat)/2)
93    
94     def dist(self,other):     def dist(self,other):
95        return math.sqrt((self.phi-other.phi)**2+(self.theta-other.theta)**2)        return math.sqrt((self.long-other.long)**2+(self.lat-other.lat)**2)
96    
97  class RegionOnEarthSurface:  class RegionOnEarthSurface:
98       """       """
# Line 106  class RegionOnEarthSurface: Line 104  class RegionOnEarthSurface:
104       NORTH=1       NORTH=1
105       EAST=2       EAST=2
106       SOUTH=3       SOUTH=3
107       def __init__(self,south_east=PointOnEarthSurface(0.,0.),north_west=PointOnEarthSurface(10.,10.),resolution=1.):       def __init__(self,west_south=PointOnEarthSurface(),east_north=PointOnEarthSurface(),resolution=1.):
108            if resolution<=0:            if resolution<=0:
109                raise ValueError,"resolution must be positive"                raise ValueError,"resolution must be positive"
110            if  south_east.phi>=north_west.phi:            if  west_south.long>=east_north.long:
111                raise ValueError,"east south corner must be further east than west north corner"                raise ValueError,"south west corner must be further east than west north corner"
112            if  south_east.theta>=north_west.theta:            if  west_south.lat>=east_north.lat:
113                raise ValueError,"east south corner must be further down than west north corner"                raise ValueError,"east south corner must be further down than west north corner"
114            if  north_west.theta-south_east.theta < resolution/2:            if  east_north.lat-west_south.lat < resolution/2:
115                raise ValueError,"latitude length of region must be 2*bigger then resolution"                raise ValueError,"latitude length of region must be 2*bigger then resolution"
116            if  north_west.phi-south_east.phi < resolution/2:            if  east_north.long-west_south.long < resolution/2:
117                raise ValueError,"longitude length of region must be 2*bigger then resolution"                raise ValueError,"longitude length of region must be 2*bigger then resolution"
118    
119            self.south_east=south_east            self.west_south=west_south
120            self.north_west=north_west            self.east_north=east_north
121            self.resolution=resolution            self.resolution=resolution
122    
123       def __str__(self):       def __str__(self):
124         return "RegionOnEarthSurface between %s and %s"%(str(self.south_east),str(self.north_west))         return "RegionOnEarthSurface between %s and %s"%(str(self.west_south),str(self.east_north))
125    
126       def isOnFace(self,p):       def isOnFace(self,p):
127           return self.isOnThisFace(p,self.SOUTH) or self.isOnThisFace(p,self.NORTH) or self.isOnThisFace(p,self.WEST) or self.isOnThisFace(p,self.EAST)           return self.isOnThisFace(p,self.SOUTH) or self.isOnThisFace(p,self.NORTH) or self.isOnThisFace(p,self.WEST) or self.isOnThisFace(p,self.EAST)
128       def isOnThisFace(self,p,face):       def isOnThisFace(self,p,face):
129          if face==self.WEST:          if face==self.WEST:
130             return self.south_east.phi==p.phi             return self.west_south.long==p.long
131          if face==self.EAST:          if face==self.EAST:
132           return p.phi==self.north_west.phi           return p.long==self.east_north.long
133          if face==self.SOUTH:          if face==self.SOUTH:
134           return self.south_east.theta==p.theta           return self.west_south.lat==p.lat
135          if face==self.NORTH:          if face==self.NORTH:
136             return p.theta==self.north_west.theta             return p.lat==self.east_north.lat
137    
138       def isBoxVertex(self,p):       def isBoxVertex(self,p):
139           return ( self.isOnThisFace(p,self.WEST) and  self.isOnThisFace(p,self.SOUTH) ) or \           return ( self.isOnThisFace(p,self.WEST) and  self.isOnThisFace(p,self.SOUTH) ) or \
# Line 146  class RegionOnEarthSurface: Line 144  class RegionOnEarthSurface:
144    
145       def getFace(self,p):       def getFace(self,p):
146          # order is critical          # order is critical
147          if self.south_east.phi==p.phi: return self.WEST          if self.west_south.long==p.long: return self.WEST
148          if p.theta==self.north_west.theta: return self.NORTH          if p.lat==self.east_north.lat: return self.NORTH
149          if p.phi==self.north_west.phi: return self.EAST          if p.long==self.east_north.long: return self.EAST
150          if self.south_east.theta==p.theta: return self.SOUTH          if self.west_south.lat==p.lat: return self.SOUTH
151    
152       def comparePointsOnFace(self,p0,p1):       def comparePointsOnFace(self,p0,p1):
153           f0=self.getFace(p0)           f0=self.getFace(p0)
# Line 161  class RegionOnEarthSurface: Line 159  class RegionOnEarthSurface:
159              return 1              return 1
160           else:           else:
161              if f0 == self.WEST:              if f0 == self.WEST:
162                  if p0.theta<p1.theta:                  if p0.lat<p1.lat:
163                     return -1                     return -1
164                  elif p0.theta>p1.theta:                  elif p0.lat>p1.lat:
165                     return 1                     return 1
166                  else:                  else:
167                     return 0                     return 0
168              elif f0 == self.EAST:              elif f0 == self.EAST:
169                  if p0.theta<p1.theta:                  if p0.lat<p1.lat:
170                     return 1                     return 1
171                  elif p0.theta>p1.theta:                  elif p0.lat>p1.lat:
172                     return -1                     return -1
173                  else:                  else:
174                     return 0                     return 0
175              elif f0 == self.NORTH:              elif f0 == self.NORTH:
176                  if p0.phi<p1.phi:                  if p0.long<p1.long:
177                     return -1                     return -1
178                  elif p0.phi>p1.phi:                  elif p0.long>p1.long:
179                     return 1                     return 1
180                  else:                  else:
181                     return 0                     return 0
182              else:              else:
183                  if p0.phi<p1.phi:                  if p0.long<p1.long:
184                     return 1                     return 1
185                  elif p0.phi>p1.phi:                  elif p0.long>p1.long:
186                     return -1                     return -1
187                  else:                  else:
188                     return 0                     return 0
189    
190    
191       def isInRegion(self,p):       def isInRegion(self,p):
192           return  self.south_east.phi<=p.phi \           return  self.west_south.long<=p.long \
193               and p.phi<=self.north_west.phi \               and p.long<=self.east_north.long \
194               and self.south_east.theta<=p.theta \               and self.west_south.lat<=p.lat \
195               and p.theta<=self.north_west.theta               and p.lat<=self.east_north.lat
196    
197       def cutSegment(self,p0,p1):       def cutSegment(self,p0,p1):
198              t=None              t=None
199              p=None              p=None
200              tmp=self.interseptOfSegment(p0,p1,d=0,v=self.south_east.phi)              tmp=self.interseptOfSegment(p0,p1,d=0,v=self.west_south.long)
201              if not tmp==None:              if not tmp==None:
202                  p_tmp=p0.split(p1,tmp)                  p_tmp=p0.split(p1,tmp)
203                  if self.isInRegion(p_tmp) and t<tmp:                  if self.isInRegion(p_tmp) and t<tmp:
204                     t=tmp                     t=tmp
205                     p=p_tmp                     p=p_tmp
206    
207              tmp=self.interseptOfSegment(p0,p1,d=0,v=self.north_west.phi)              tmp=self.interseptOfSegment(p0,p1,d=0,v=self.east_north.long)
208              if not tmp==None:              if not tmp==None:
209                  p_tmp=p0.split(p1,tmp)                  p_tmp=p0.split(p1,tmp)
210                  if self.isInRegion(p_tmp) and t<tmp:                  if self.isInRegion(p_tmp) and t<tmp:
211                     t=tmp                     t=tmp
212                     p=p_tmp                     p=p_tmp
213    
214              tmp=self.interseptOfSegment(p0,p1,d=1,v=self.south_east.theta)              tmp=self.interseptOfSegment(p0,p1,d=1,v=self.west_south.lat)
215              if not tmp==None:              if not tmp==None:
216                  p_tmp=p0.split(p1,tmp)                  p_tmp=p0.split(p1,tmp)
217                  if self.isInRegion(p_tmp) and t<tmp:                  if self.isInRegion(p_tmp) and t<tmp:
218                     t=tmp                     t=tmp
219                     p=p_tmp                     p=p_tmp
220    
221              tmp=self.interseptOfSegment(p0,p1,d=1,v=self.north_west.theta)              tmp=self.interseptOfSegment(p0,p1,d=1,v=self.east_north.lat)
222              if not tmp==None:              if not tmp==None:
223                  p_tmp=p0.split(p1,tmp)                  p_tmp=p0.split(p1,tmp)
224                  if self.isInRegion(p_tmp) and t<tmp:                  if self.isInRegion(p_tmp) and t<tmp:
# Line 233  class RegionOnEarthSurface: Line 231  class RegionOnEarthSurface:
231           find t isn [0,1] such that p0+t*(p0-p1) cuts x[d]=v. if it does nit exist None is returned.           find t isn [0,1] such that p0+t*(p0-p1) cuts x[d]=v. if it does nit exist None is returned.
232           """           """
233           if d==0:           if d==0:
234             a=p0.phi             a=p0.long
235             b=p1.phi             b=p1.long
236           else:           else:
237             a=p0.theta             a=p0.lat
238             b=p1.theta             b=p1.lat
239           if b==a:           if b==a:
240             if a==v:             if a==v:
241                t=0                t=0
# Line 321  class Polyline: Line 319  class Polyline:
319           for i in range(len(self.list_of_coordinates)-1):           for i in range(len(self.list_of_coordinates)-1):
320               p0=self.list_of_coordinates[i]               p0=self.list_of_coordinates[i]
321               p1=self.list_of_coordinates[i+1]               p1=self.list_of_coordinates[i+1]
322               integ+=(p1.theta-p0.theta)*(p1.phi+p0.phi)-(p1.phi-p0.phi)*(p1.theta+p0.theta)               integ+=(p1.lat-p0.lat)*(p1.long+p0.long)-(p1.long-p0.long)*(p1.lat+p0.lat)
323           if integ>=0.:           if integ>=0.:
324              return 1.              return 1.
325           else:           else:
# Line 343  class Coastline: Line 341  class Coastline:
341           for pl in self.polylines:           for pl in self.polylines:
342               out+="\n"+str(pl)               out+="\n"+str(pl)
343           return out           return out
344       def makeTriangulation(self,south_east_is_water=True,south_west_is_water=True,north_west_is_water=True,north_east_is_water=True):       def makeTriangulation(self,west_south_is_water=True,east_south_is_water=True,west_north_is_water=True,east_north_is_water=True):
345           self.clean()           self.clean()
346           vertices=[]           vertices=[]
347           segments=[]           segments=[]
# Line 371  class Coastline: Line 369  class Coastline:
369                         vertices+=v_tmp                         vertices+=v_tmp
370                         segments+=s_tmp                         segments+=s_tmp
371                         # find a point in the loop:                         # find a point in the loop:
372                         d_phi=v_tmp[1].phi-v_tmp[0].phi                         d_long=v_tmp[1].long-v_tmp[0].long
373                         d_theta=v_tmp[1].theta-v_tmp[0].theta                         d_lat=v_tmp[1].lat-v_tmp[0].lat
374                         l=math.sqrt(d_phi**2+d_theta**2)                         l=math.sqrt(d_long**2+d_lat**2)
375                         mid=v_tmp[0].midPoint(v_tmp[1])                         mid=v_tmp[0].midPoint(v_tmp[1])
376                         n_phi=-d_theta/l                         n_long=-d_lat/l
377                         n_theta=d_phi/l                         n_lat=d_long/l
378                         s=self.region.resolution                         s=self.region.resolution
379                         for seg in s_tmp:                         for seg in s_tmp:
380                           p0=vertices[seg[0]]                           p0=vertices[seg[0]]
381                           p1=vertices[seg[1]]                           p1=vertices[seg[1]]
382                           a_phi=p1.phi-p0.phi                           a_long=p1.long-p0.long
383                           a_theta=p1.theta-p0.theta                           a_lat=p1.lat-p0.lat
384                           d=a_theta*n_phi-a_phi*n_theta                           d=a_lat*n_long-a_long*n_lat
385                           if abs(d)>0.:                           if abs(d)>0.:
386                              t=((mid.theta-p0.theta)*n_phi-(mid.phi-p0.phi)*n_theta)/d                              t=((mid.lat-p0.lat)*n_long-(mid.long-p0.long)*n_lat)/d
387                              if 0<=t and t<=1:                              if 0<=t and t<=1:
388                                  s_tmp=((mid.theta-p0.theta)*a_phi-(mid.phi-p0.phi)*a_theta)/d                                  s_tmp=((mid.lat-p0.lat)*a_long-(mid.long-p0.long)*a_lat)/d
389                                  if s_tmp>EPS: s=min(s,s_tmp)                                  if s_tmp>EPS: s=min(s,s_tmp)
390                         h=PointOnEarthSurface(mid.phi+s/2*n_phi,mid.theta+s/2*n_theta)                         h=PointOnEarthSurface(long=mid.long+s/2*n_long,lat=mid.lat+s/2*n_lat)
391                         holes.append(h)                         holes.append(h)
392                  else:                  else:
393                     if len(short_pl)>1:                     if len(short_pl)>1:
# Line 401  class Coastline: Line 399  class Coastline:
399                            vertices.append(short_pl[i])                            vertices.append(short_pl[i])
400           # put into the bounding box:           # put into the bounding box:
401           new_vertices=[]           new_vertices=[]
402           if south_east_is_water:           if west_south_is_water:
403              new_vertices.append(PointOnEarthSurface(self.region.south_east.phi,self.region.south_east.theta))              new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long,lat=self.region.west_south.lat))
404           if north_east_is_water:           if east_south_is_water:
405              new_vertices.append(PointOnEarthSurface(self.region.south_east.phi,self.region.north_west.theta))              new_vertices.append(PointOnEarthSurface(long=self.region.east_north.long,lat=self.region.west_south.lat))
406           if north_west_is_water:           if west_north_is_water:
407              new_vertices.append(PointOnEarthSurface(self.region.north_west.phi,self.region.north_west.theta))              new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long,lat=self.region.east_north.lat))
408           if south_west_is_water:           if east_north_is_water:
409              new_vertices.append(PointOnEarthSurface(self.region.north_west.phi,self.region.south_east.theta))              new_vertices.append(PointOnEarthSurface(long=self.region.east_north.long,lat=self.region.east_north.lat))
410    
411           # add new vertices if they don't exist yet           # add new vertices if they don't exist yet
412           for q in new_vertices:           for q in new_vertices:
# Line 421  class Coastline: Line 419  class Coastline:
419                  vertices_on_face.append(q)                  vertices_on_face.append(q)
420           vertices_on_face.sort(self.region.comparePointsOnFace)           vertices_on_face.sort(self.region.comparePointsOnFace)
421           index=0           index=0
422           walking_on_water=south_east_is_water           walking_on_water=west_south_is_water
423           l=len(vertices_on_face)           l=len(vertices_on_face)
424           while index<l:           while index<l:
425               p1=vertices_on_face[(index+1)%l]               p1=vertices_on_face[(index+1)%l]
# Line 553  class EarthTriangulation: Line 551  class EarthTriangulation:
551             poly_file=self.fn+".poly"             poly_file=self.fn+".poly"
552             f=open(poly_file,"w")             f=open(poly_file,"w")
553             f.writelines("%d %d %d %d\n"%(len(vertices),2,0,0))             f.writelines("%d %d %d %d\n"%(len(vertices),2,0,0))
554             for i in range(len(vertices)): f.writelines("%d %e %e\n"%(i,vertices[i].phi,vertices[i].theta))             for i in range(len(vertices)): f.writelines("%d %e %e\n"%(i,vertices[i].long,vertices[i].lat))
555             f.writelines("%d %d\n"%(len(segments),0))             f.writelines("%d %d\n"%(len(segments),0))
556             for i in range(len(segments)): f.writelines("%d %d %d\n"%(i,segments[i][0],segments[i][1]))             for i in range(len(segments)): f.writelines("%d %d %d\n"%(i,segments[i][0],segments[i][1]))
557             f.writelines("%d\n"%(len(holes)))             f.writelines("%d\n"%(len(holes)))
558             for i in range(len(holes)):  f.writelines("%d %e %e\n"%(i,holes[i].phi,holes[i].theta))             for i in range(len(holes)):  f.writelines("%d %e %e\n"%(i,holes[i].long,holes[i].lat))
559             f.close()             f.close()
560             # start mesh generator:             # start mesh generator:
561             os.system(self.GENERATOR%(resolution**2,poly_file))             os.system(self.GENERATOR%(resolution**2,poly_file))
# Line 736  class OceanRegion(Model): Line 734  class OceanRegion(Model):
734                line=line.strip()                line=line.strip()
735                if line[:7]=="#region":                if line[:7]=="#region":
736                   data=line[line.index("[")+1:line.index("]")].split(",")                   data=line[line.index("[")+1:line.index("]")].split(",")
737                   reg=RegionOnEarthSurface(PointOnEarthSurface(self.south,self.west),PointOnEarthSurface(self.north,self.east),self.resolution)                   reg=RegionOnEarthSurface(PointOnEarthSurface(lat=self.south,long=self.west),PointOnEarthSurface(lat=self.north,long=self.east),self.resolution)
738                   self.trace(str(reg))                   self.trace(str(reg))
739                   self.coastline=Coastline(region=reg,name=data_name)                   self.coastline=Coastline(region=reg,name=data_name)
740                elif line.find("Shore Bin")>-1:                elif line.find("Shore Bin")>-1:
# Line 745  class OceanRegion(Model): Line 743  class OceanRegion(Model):
743                     segs=[]                     segs=[]
744                if not (line=="" or line[0]=="#" or line[0]==">") :                if not (line=="" or line[0]=="#" or line[0]==">") :
745                    x=line.split()                    x=line.split()
746                    segs.append(PointOnEarthSurface(float(x[0]),float(x[1])))                    segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))
747                line=f.readline()                line=f.readline()
748             self.coastline.append(Polyline(segs,name))             self.coastline.append(Polyline(segs,name))
749             d=self.bathymetry_data.interpolate([[self.east,self.south],[self.west,self.south],[self.east,self.north],[self.west,self.north]])             d=self.bathymetry_data.interpolate([[self.east,self.south],[self.west,self.south],[self.east,self.north],[self.west,self.north]])
750             self.domain=self.coastline.makeTriangulation(south_east_is_water=d[0]<0,south_west_is_water=d[1]<0,north_east_is_water=d[2]<0,north_west_is_water=d[3]<0).getFinleyDomain()             self.domain=self.coastline.makeTriangulation(east_south_is_water=d[0]<=0,west_south_is_water=d[1]<=0,east_north_is_water=d[2]<=0,west_north_is_water=d[3]<=0).getFinleyDomain()
751             self.bathymetry=maximum(-self.bathymetry_data.interpolate(Function(self.domain).getX()),0.)             self.bathymetry=maximum(-self.bathymetry_data.interpolate(Function(self.domain).getX()),0.)
752    
753    
# Line 776  class TsunamiSource(Model): Line 774  class TsunamiSource(Model):
774             """             """
775             beta=math.sqrt(-math.log(self.GAMMA))/self.decay_zone             beta=math.sqrt(-math.log(self.GAMMA))/self.decay_zone
776             x=self.domain.getX()             x=self.domain.getX()
777             x0=x[0]             x_long=x[0]
778             x1=x[1]             x_lat=x[1]
779             mid0=(self.start_long+self.end_long)/2             mid_long=(self.start_long+self.end_long)/2
780             mid1=(self.start_lat+self.end_lat)/2             mid_lat=(self.start_lat+self.end_lat)/2
781             dist=math.sqrt((mid0-self.end_long)**2+(mid1-self.end_lat)**2)             dist=math.sqrt((mid_long-self.end_long)**2+(mid_lat-self.end_lat)**2)
782             a=(self.start_long-mid0)/dist             a=(self.start_long-mid_long)/dist
783             b=(self.end_long-mid1)/dist             b=(self.start_lat-mid_lat)/dist
784             self.trace("source length = %s"%(dist*2.))             self.trace("source length = %s"%(dist*2.))
785               x_long_hat= a*(x_long-mid_long)+b*(x_lat-mid_lat)
786             x0_hat= a*(x0-mid0)+b*(x1-mid1)             x_lat_hat=-b*(x_long-mid_long)+a*(x_lat-mid_lat)
787             x1_hat=-b*(x0-mid0)+a*(x1-mid1)             # x_lat-direction
788             # x1-direction             s=abs(x_lat_hat)-self.width
            s=abs(x1_hat)-self.width  
789             m=s.whereNegative()             m=s.whereNegative()
790             f1=(1.-m)*exp(-(s*beta)**2)+m             f1=(1.-m)*exp(-(s*beta)**2)+m
791             # x0-direction             # x_long-direction
792             s=abs(x0_hat)-dist             s=abs(x_long_hat)-dist
793             m=s.whereNegative()             m=s.whereNegative()
794             f0=(1.-m)*exp(-(s*beta)**2)+m             f0=(1.-m)*exp(-(s*beta)**2)+m
795             self.wave_height=f1*f0*self.amplitude             self.wave_height=f1*f0*self.amplitude
# Line 931  if __name__=="__main__": Line 928  if __name__=="__main__":
928     from esys.modellib.input import Sequencer     from esys.modellib.input import Sequencer
929    
930     oc=OceanRegionCollector()     oc=OceanRegionCollector()
931     oc.resolution=0.5     oc.resolution=2
932     oc.south=-45.5     oc.south=-45.5
933     oc.north=-5.4     oc.north=-5.4
934     oc.east=161.1     oc.east=161.1
# Line 954  if __name__=="__main__": Line 951  if __name__=="__main__":
951     src=TsunamiSource()     src=TsunamiSource()
952     src.domain=Link(oreg,"domain")     src.domain=Link(oreg,"domain")
953     src.start_lat=-10.     src.start_lat=-10.
    src.start_long=110.  
954     src.end_lat=-12.     src.end_lat=-12.
955       src.start_long=110.
956     src.end_long=120.     src.end_long=120.
957     src.width=0.1     src.width=0.1
958     src.decay_zone=0.01     src.decay_zone=0.01
# Line 968  if __name__=="__main__": Line 965  if __name__=="__main__":
965     ts.bathymetry=Link(oreg,"bathymetry")     ts.bathymetry=Link(oreg,"bathymetry")
966    
967     sq=Sequencer()     sq=Sequencer()
968     sq.t_end=1000.     sq.t_end=100000.
969    
970     sm=SurfMovie()     sm=SurfMovie()
971     sm.bathymetry=Link(b,"bathymetry")     sm.bathymetry=Link(b,"bathymetry")
972     sm.wave_height=Link(ts,"wave_height")     sm.wave_height=Link(ts,"wave_height")
973     sm.coastline=Link(oreg,"coastline")     sm.coastline=Link(oreg,"coastline")
974     sm.t=Link(sq,"t")     sm.t=Link(sq,"t")
975     sm.dt=100.     sm.dt=5000.
976     sm.filename="movie.mpg"     sm.filename="movie.mpg"
977        
978     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
Legend:
Removed from v.247  
changed lines
  Added in v.251
  ViewVC Help
Powered by ViewVC 1.1.26