/[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 4097 by caltinay, Fri Dec 7 01:18:35 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 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 508  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 518  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 562  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=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 = 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 583  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.5  
         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.4097  
changed lines
  Added in v.4108

  ViewVC Help
Powered by ViewVC 1.1.26