/[escript]/trunk/downunder/py_src/datasources.py
ViewVC logotype

Diff of /trunk/downunder/py_src/datasources.py

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

revision 4080 by jfenwick, Mon Nov 19 01:45:38 2012 UTC revision 4108 by caltinay, Thu Dec 13 06:38:11 2012 UTC
# Line 22  __license__="""Licensed under the Open S Line 22  __license__="""Licensed under the Open S
22  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
23  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
24    
25  __all__ = ['DataSource','ErMapperData','SyntheticData','SmoothAnomaly']  __all__ = ['simpleBackgroundMagneticField', 'DataSource','ErMapperData','SyntheticFeatureData','SmoothAnomaly']
26    
27  import logging  import logging
28  import numpy as np  import numpy as np
# Line 78  def LatLonToUTM(lon, lat, wkt_string=Non Line 78  def LatLonToUTM(lon, lat, wkt_string=Non
78      x,y=pyproj.transform(p_src, p_dest, lon, lat)      x,y=pyproj.transform(p_src, p_dest, lon, lat)
79      return x,y      return x,y
80    
81    def simpleBackgroundMagneticField(latitude, longitude=0.):
82            theta = (90-latitude)/180.*np.pi
83            B_0=U.Mu_0  * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi *  U.R_Earth**3)
84            B_theta= B_0 * sin(theta)
85            B_r= 2 * B_0 * cos(theta)
86            return B_r, B_theta, 0.
87    
88  class DataSource(object):  class DataSource(object):
89      """      """
# Line 85  class DataSource(object): Line 91  class DataSource(object):
91      This is an abstract base class that implements common functionality.      This is an abstract base class that implements common functionality.
92      Methods to be overwritten by subclasses are marked as such.      Methods to be overwritten by subclasses are marked as such.
93      This class assumes 2D data which is mapped to a slice of a 3D domain.      This class assumes 2D data which is mapped to a slice of a 3D domain.
94      For other setups override the `_createDomain` method.      For other setups override the methods as required.
95      """      """
96    
97      GRAVITY, MAGNETIC = list(range(2))      GRAVITY, MAGNETIC = list(range(2))
# Line 96  class DataSource(object): Line 102  class DataSource(object):
102          """          """
103          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
104          self.__subsampling_factor=1          self.__subsampling_factor=1
105            self.__background_magnetic_field=None
106    
107      def getDataExtents(self):      def getDataExtents(self):
108          """          """
# Line 410  class NetCdfData(DataSource): Line 417  class NetCdfData(DataSource):
417          #lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]          #lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]
418          lon_range=longitude.data.min(),longitude.data.max()          lon_range=longitude.data.min(),longitude.data.max()
419          lat_range=latitude.data.min(),latitude.data.max()          lat_range=latitude.data.min(),latitude.data.max()
420          lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)          if lon_range[1]<180:
421                lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
422          lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]          lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]
423          f.close()          f.close()
424    
# Line 507  class SmoothAnomaly(SourceFeature): Line 515  class SmoothAnomaly(SourceFeature):
515          return m          return m
516    
517  ##############################################################################  ##############################################################################
518  class SyntheticData(DataSource):  class SyntheticFeatureData(DataSource):
519      def __init__(self, datatype, DIM, NE, l, features):      def __init__(self, datatype, DIM, NE, l, features, B_b=None):
520          super(SyntheticData,self).__init__()          super(SyntheticFeatureData,self).__init__()
521          if not datatype in [self.GRAVITY,self.MAGNETIC]:          if not datatype in [self.GRAVITY,self.MAGNETIC]:
522              raise ValueError("Invalid value for datatype parameter")              raise ValueError("Invalid value for datatype parameter")
523          self.__datatype = datatype          self.__datatype = datatype
# Line 517  class SyntheticData(DataSource): Line 525  class SyntheticData(DataSource):
525          self.__origin = [0]*(DIM-1)          self.__origin = [0]*(DIM-1)
526          self.__nPts = [NE]*(DIM-1)          self.__nPts = [NE]*(DIM-1)
527          self.__delta = [float(l)/NE]*(DIM-1)          self.__delta = [float(l)/NE]*(DIM-1)
528            self.__B_b =None
529          self.DIM=DIM          self.DIM=DIM
530          self.NE=NE          self.NE=NE
531          self.l=l          self.l=l
532            # this is for Cartesian (FIXME ?)
533            if datatype  ==  self.MAGNETIC:
534                if self.DIM<3:
535                   self.__B_b =  np.array([-B_b[2],  -B_b[0]])
536                else:
537                   self.__B_b = ([-B_b[1],  -B_b[2],  -B_b[0]])
538    
539      def getDataExtents(self):      def getDataExtents(self):
540          return (list(self.__origin), list(self.__nPts), list(self.__delta))          return (list(self.__origin), list(self.__nPts), list(self.__delta))
# Line 561  class SyntheticData(DataSource): Line 576  class SyntheticData(DataSource):
576                  m=f.getMask(x)                  m=f.getMask(x)
577                  k_ref = k_ref * (1-m) + f.getValue(x) * m                  k_ref = k_ref * (1-m) + f.getValue(x) * m
578              self._k=k_ref              self._k=k_ref
579              B_b = self.getBackgroundMagneticField()              pde.setValue(A=kronecker(domain), X=k_ref*self.__B_b, q=m_psi_ref)
             pde.setValue(A=kronecker(domain), X=(1+k_ref)*B_b, q=m_psi_ref)  
580          pde.setSymmetryOn()          pde.setSymmetryOn()
581          psi_ref=pde.getSolution()          psi_ref=pde.getSolution()
582          del pde          del pde
583          if self.getDataType()==DataSource.GRAVITY:          if self.getDataType()==DataSource.GRAVITY:
584              data = -grad(psi_ref, ReducedFunction(domain))              data = -grad(psi_ref, ReducedFunction(domain))
585          else:          else:
586              data = (1+self._k)*B_b-grad(psi_ref, ReducedFunction(domain))              data = self._k*self.__B_b-grad(psi_ref, ReducedFunction(domain))
587    
588          sigma=1.          sigma=1.
589          # limit mask to non-padding in horizontal area          # limit mask to non-padding in horizontal area
# Line 582  class SyntheticData(DataSource): Line 596  class SyntheticData(DataSource):
596                  * whereNonPositive(x[self.DIM-1]-spacing[self.DIM-1])                  * whereNonPositive(x[self.DIM-1]-spacing[self.DIM-1])
597          return data,sigma          return data,sigma
598    
599      def getBackgroundMagneticField(self):  
         #FIXME:  
         latitude=-28.0  
         theta = (90-latitude)/180.*np.pi  
         B_0=U.Mu_0  * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi *  U.R_Earth**3)  
         B_theta= B_0 * sin(theta)  
         B_r= 2 * B_0 * cos(theta)  
         if self.DIM<3:  
             return np.array([0.,  -B_r])  
         else:  
             return np.array([-B_theta, 0.,  -B_r])  
600    

Legend:
Removed from v.4080  
changed lines
  Added in v.4108

  ViewVC Help
Powered by ViewVC 1.1.26