Package SatStress :: Module SatStress
[frames] | no frames]

Source Code for Module SatStress.SatStress

   1  """ 
   2  A framework for calculating the surface stresses at a particular place and time 
   3  on a satellite resulting from one or more tidal potentials. 
   4   
   5  1 Input and Output 
   6  ================== 
   7  Because C{SatStress} is a "library" module, it doesn't do a lot of input and 
   8  output - it's mostly about doing calculations.  It does need to read in the 
   9  specification of a L{Satellite} object though, and it can write the same kind 
  10  of specification out.  To do this, it uses name-value files, and a function 
  11  called L{nvf2dict}, which creates a Python dictionary (or "associative array"). 
  12   
  13  A name-value file is just a file containing a bunch of name-value pairs, like:: 
  14   
  15    ORBIT_ECCENTRICITY = 0.0094   # e must be < 0.25 
  16   
  17  It can also contain comments to enhance human readability (anything following a 
  18  '#' on a line is ignored, as with the note in the line above). 
  19   
  20  2 Satellites 
  21  ============ 
  22  Obviously if we want to calculate the stresses on the surface of a satellite, 
  23  we need to define the satellite, this is what the L{Satellite} object does. 
  24   
  25  2.1 Specifying a Satellite 
  26  -------------------------- 
  27    In order to specify a satellite, we need: 
  28   
  29      - an ID of some kind for the planet/satellite pair of interest 
  30      - the charactaristics of the satellite's orbital environment 
  31      - the satellite's internal structure and material properties 
  32      - the forcings to which the satellite is subjected 
  33   
  34    From a few basic inputs, we can calculate many derived characteristics, such 
  35    as the satellite's orbital period or the surface gravity. 
  36   
  37    The internal structure and material properties are specified by a series of 
  38    concentric spherical shells (layers), each one being homogeneous throughout 
  39    its extent.  Given the densities and thicknesses of the these layers, we can 
  40    calculate the satellite's overall size, mass, density, etc. 
  41   
  42    Specifying a tidal forcing may be simple or complex.  For instance, the 
  43    L{Diurnal} forcing depends only on the orbital eccentricity (and other 
  44    orbital parameters already supplied), and the L{NSR} forcing requires only 
  45    the addition of the non-synchronous rotation period of the shell.  Specifying 
  46    an arbitrary true polar wander trajectory would be much more complex. 
  47   
  48    For the moment, becuase we are only including simple forcings, their 
  49    specifying parameters are read in from the satellite definition file.  If 
  50    more, and more complex forcings are eventually added to the model, their 
  51    specification will probably be split into a separate input file. 
  52   
  53  2.2 Internal Structure and Love Numbers 
  54  --------------------------------------- 
  55    SatStress treats the solid portions of the satellite as U{viscoelastic 
  56    Maxwell solids <http://en.wikipedia.org/wiki/Maxwell_material>}, that respond 
  57    differently to forcings having different frequencies (S{omega}).  Given the a 
  58    specification of the internal structure and material properties of a 
  59    satellite as a series of layers, and information about the tidal forcings the 
  60    body is subject to, it's possible to calculate appropriate Love numbers, 
  61    which describe how the body responds to a change in the gravitational 
  62    potential. 
  63   
  64    Currently the calculation of Love numbers is done by an external program 
  65    written in Fortran by John Wahr and others, with roots reaching deep into the 
  66    Dark Ages of computing.  As that code (or another Love number code) is more 
  67    closely integrated with the model, the internal structure of the satellite 
  68    will become more flexible, but for the moment, we are limited to assuming 
  69    a 4-layer structure: 
  70   
  71      - B{C{ICE_UPPER}}: The upper portion of the shell (cooler, stiffer) 
  72      - B{C{ICE_LOWER}}: The lower portion of the shell (warmer, softer) 
  73      - B{C{OCEAN}}: An inviscid fluid decoupling the shell from the core. 
  74      - B{C{CORE}}: The silicate interior of the body. 
  75   
  76  3 Stresses 
  77  ========== 
  78   
  79  C{SatStress} can calculate the following stress fields: 
  80   
  81    1. B{L{Diurnal}}: stresses arising from an eccentric orbit, having a 
  82       forcing frequency equal to the orbital frequency. 
  83   
  84    2. B{L{NSR}}: stresses arising due to the faster-than-synchronous rotation 
  85       of a floating shell that is decoupled from the satellite's interior by a 
  86       fluid layer (an ocean). 
  87   
  88  The expressions defining these stress fields are derived in "Modeling Stresses 
  89  on Satellites due to Non-Synchronous Rotation and Orbital Eccentricity Using 
  90  Gravitational Potential Theory" (U{preprint, 15MB PDF 
  91  <http://satstress.googlecode.com/files/Wahretal2008.pdf>}) by Wahr et al. 
  92  (submitted to I{Icarus}, in March, 2008). 
  93   
  94  3.1 Stress Fields Live in L{StressDef} Objects 
  95  ---------------------------------------------- 
  96    Each of the above stress fields is defined by a similarly named L{StressDef} 
  97    object.  These objects contain the formulae necessary to calculate the 
  98    surface stress.  The expressions for the stresses depend on many parameters 
  99    which are defined within the L{Satellite} object, and so to create a 
 100    L{StressDef} object, you need to provide a L{Satellite} object. 
 101   
 102    There are many formulae which are identical for both the L{NSR} and 
 103    L{Diurnal} stress fields, and so instead of duplicating them in both classes, 
 104    they reside in the L{StressDef} I{base class}, from which all L{StressDef} 
 105    objects inherit many properties. 
 106   
 107    The main requirement for each L{StressDef} object is that it must define the 
 108    three components of the stress tensor S{tau}: 
 109   
 110      - C{Ttt} (S{tau}_S{theta}S{theta}) the north-south (latitudinal) component 
 111   
 112      - C{Tpt} (S{tau}_S{phi}S{theta} = S{tau}_S{theta}S{phi}) the shear 
 113        component 
 114   
 115      - C{Tpp} (S{tau}_S{phi}S{phi}) the east-west (longitudinal) component 
 116   
 117  3.2 Stress Calculations are Performed by L{StressCalc} Objects 
 118  -------------------------------------------------------------- 
 119    Once you've I{instantiated} a L{StressDef} object, or several of them (one 
 120    for each stress you want to include), you can compose them together into a 
 121    L{StressCalc} object, which will actually do calculations at given points on 
 122    the surface, and given times, and return a 2x2 matrix containing the 
 123    resulting stress tensor (each component of which is the sum of all of the 
 124    corresponding components of the stress fields that were used to instantiated 
 125    the L{StressCalc} object). 
 126   
 127    This is (hopefully) easier than it sounds.  With the following few lines, you 
 128    can construct a satellite, do a single calculation on its surface, and see 
 129    what it looks like: 
 130   
 131    >>> from SatStress.SatStress import * 
 132    >>> the_sat = Satellite(open("input/Europa.satellite")) 
 133    >>> the_stresses = StressCalc([Diurnal(the_sat), NSR(the_sat)]) 
 134    >>> Tau = the_stresses.tensor(theta=pi/4.0, phi=pi/3.0, t=10000) 
 135    >>> print(Tau) 
 136   
 137    The C{test} program included in the SatStress distribution shows a slightly 
 138    more complex example, which should be enough to get you started using the 
 139    package. 
 140   
 141  3.3 Extending the Model   
 142  ----------------------- 
 143    Other stress fields can (and hopefully will!), be added easily, so long as 
 144    they use the same mathematical definition of the membrane stress tensor 
 145    (S{tau}), as a function of co-latitude (S{theta}) (measured south from the 
 146    north pole), east-positive longitude (S{phi}), measured from the meridian 
 147    on the satellite which passes through the point on the satellite directly 
 148    beneath the parent planet (assuming a synchronously rotating satellite), 
 149    and time (B{M{t}}), defined as seconds elapsed since pericenter. 
 150   
 151    This module could also potentially be extended to also calculate the 
 152    surface strain (S{epsilon}) and displacement (B{M{s}}) fields, or to 
 153    calculate the stresses at any point within the satellite. 
 154   
 155    @group Exceptions (error handling classes): *Error 
 156  """ 
 157   
 158  # regular expressions 
 159  import re 
 160   
 161  # required for command line parsing 
 162  import sys 
 163   
 164  # Scientific functions... like complex exponentials 
 165  import scipy 
 166  import numpy 
 167   
 168  # Physical constants, like the Newton's Gravitational Constant Big G 
 169  import physcon as pc 
 170   
 171  ############################## 
 172  #      HELPER FUNCTIONS      # 
 173  ############################## 
174 -def nvf2dict(nvf, comment='#'):
175 """ 176 Reads from a file object listing name value pairs, creating and returning a 177 corresponding Python dictionary. 178 179 The file should contain a series of name value pairs, one per line 180 separated by the '=' character, with names on the left and values on the 181 right. Blank lines are ignored, as are lines beginning with the comment 182 character (assumed to be the pound or hash character '#', unless otherwise 183 specified). End of line comments are also allowed. String values should 184 not be quoted in the file. Names are case sensitive. 185 186 Returns a Python dictionary that uses the names as keys and the values as 187 values, and so all Python limitations on what can be used as a dictionary 188 key apply to the name fields. 189 190 Leading and trailing whitespace is stripped from all names and values, and 191 all values are returned as strings. 192 193 @param nvf: an open file object from which to read the name value pairs 194 @param comment: character which begins comments 195 @type nvf: file 196 @type comment: str 197 @return: a dictionary containing the name value pairs read in from C{nvf}. 198 @rtype: dict 199 @raise NameValueFileParseError: if a non-comment input line does not 200 contain an '=' character, or if a non-comment line has nothing but 201 whitespace preceeding or following the '=' character. 202 @raise NameValueFileDuplicateNameError: if more than one instance of the 203 same name is found in the input file C{nvf}. 204 205 """ 206 207 params = {} 208 for line in nvf: 209 210 line = line.strip() 211 212 # strip comments 213 line = re.sub("\s*%s.*" % (comment), "", line) 214 215 # if the line is all whitespace, ignore it 216 if re.match("^\s*$", line): 217 continue 218 219 # if we don't have two non-space strings separated by '=', Error! 220 if not (re.match("[^\s]*.*=.*[^\s]", line)): 221 raise NameValueFileParseError(nvf, line) 222 # split the NAME=VALUE pair into two variables if possible. 223 try: 224 (name,value) = re.split("=",line,1) 225 except ValueError: 226 raise NameValueFileParseError(nvf, line) 227 228 # remove leading and trailing whitespace 229 name = name.strip() 230 value = value.strip() 231 232 # If the file we're parsing has the same name specified more than once, 233 # then we're going to end up clobbering the earlier value here, which 234 # is probably not what we want to do. 235 if params.has_key(name): 236 raise NameValueFileDuplicateNameError(nvf, name) 237 else: 238 params[name] = value 239 240 return params
241 # end nvf2dict() 242 243 244 ############################## 245 # CLASSES # 246 ##############################
247 -class Satellite(object):
248 """An object describing the physical structure and context of a satellite. 249 250 Defines a satellite's material properties, internal structure, orbital 251 context, and the tidal forcings to which it is subjected. 252 253 @ivar sourcefile: the file object which was read in order to create 254 the L{Satellite} instance. 255 @type sourcefile: file 256 @ivar satParams: dictionary containing the name value pairs read in from 257 the input file. 258 @type satParams: dict 259 @ivar system_id: string identifying the planet/satellite system, 260 corresponds to C{SYSTEM_ID} in the input file. 261 @type system_id: str 262 @ivar orbit_semimajor_axis: semimajor axis of the satellite's orbit [m], 263 corresponds to C{ORBIT_SEMIMAJOR_AXIS} in the input file. 264 @type orbit_semimajor_axis: float 265 @ivar orbit_eccentricity: the satellite's orbital eccentricity, 266 corresponds to C{ORBIT_ECCENTRICITY} in the input file. 267 @type orbit_eccentricity: float 268 @ivar planet_mass: the mass of the planet the satellite orbits [kg], 269 corresponds to C{PLANET_MASS} in the input file. 270 @type planet_mass: float 271 @ivar nsr_period: the time it takes for the decoupled ice shell to 272 complete one full rotation [s], corresponds to C{NSR_PERIOD} in the input 273 file. May also be set to infinity (inf, infinity, INF, INFINITY). 274 @type nsr_period: float 275 @ivar num_layers: the number of layers making up the satellite, as 276 indicated by the number of keys within the satParams dictionary contain 277 the string 'LAYER_ID'. Currently this must equal 4. 278 @type num_layers: int 279 @ivar layers: a list of L{SatLayer} objects, describing the layers making 280 up the satellite. The layers are ordered from the center of the satellite 281 outward, with layers[0] corresponding to the core. 282 @type layers: list 283 """ 284
285 - def __init__(self, satFile): #
286 """ 287 Construct a Satellite object from a satFile 288 289 Required input file parameters: 290 =============================== 291 The Satellite is initialized from a name value file (as described under 292 L{nvf2dict}). The file must define the following parameters, all of which 293 are specified in SI (MKS) units. 294 295 - B{C{SYSTEM_ID}}: A string identifying the planetary system, e.g. 296 JupiterEuropa. 297 298 - B{C{PLANET_MASS}}: The mass of the planet the satellite orbits 299 [kg]. 300 301 - B{C{ORBIT_ECCENTRICITY}}: The eccentricity of the satellite's 302 orbit. Must not exceed 0.25. 303 304 - B{C{ORBIT_SEMIMAJOR_AXIS}}: The semimajor axis of the satellite's 305 orbit [m]. 306 307 - B{C{NSR_PERIOD}}: The time it takes for the satellite's icy shell 308 to undergo one full rotation [s]. If you don't want to have any 309 NSR stresses, just put INFINITY here. 310 311 An example satFile is included with the SatStress package:: 312 313 SatStress-X.Y.Z/input/Europa.satellite 314 315 Where X.Y.Z refers to the version numbers. 316 317 @param satFile: Open file object containing name value pairs specifying 318 the satellite's internal structure and orbital context, and the tidal 319 forcings to which the satellite is subjected. 320 @type satFile: file 321 @return: a Satellite object corresponding to the proffered input file. 322 @rtype: L{Satellite} 323 324 @raise NameValueFileError: if parsing of the input file fails. 325 @raise MissingSatelliteParamError: if a required input field is not 326 found within the input file. 327 @raise NonNumberSatelliteParamError: if a required input which is of 328 a numeric type is found within the file, but its value is not 329 convertable to a float. 330 @raise LoveLayerNumberError: if the number of layers specified is not 331 exactly 4. 332 @raise GravitationallyUnstableSatelliteError: if the layers are found not 333 to decrease in density from the core to the surface. 334 @raise ExcessiveSatelliteMassError: if the mass of the satellite's parent 335 planet is not at least 10 times larger than the mass of the satellite. 336 @raise LargeEccentricityError: if the satellite's orbital eccentricity is 337 greater than 0.25 338 @raise NegativeNSRPeriodError: if the NSR period of the satellite is less 339 than zero. 340 """ 341 342 # Make a note of the input file defining this satellite: 343 self.sourcefile = satFile 344 345 # Slurp the proffered parameters into a dictionary: 346 try: 347 self.satParams = nvf2dict(self.sourcefile, comment='#') 348 except NameValueFileError: 349 raise 350 351 # Try to assign our data members to the values we read in from 352 # the file. If we didn't get all the values we need, raise an 353 # exception saying so. 354 355 try: 356 self.system_id = self.satParams['SYSTEM_ID'] 357 except KeyError: 358 raise MissingSatelliteParamError(self, 'SYSTEM_ID') 359 360 try: 361 self.orbit_semimajor_axis = float(self.satParams['ORBIT_SEMIMAJOR_AXIS']) 362 except KeyError: 363 raise MissingSatelliteParamError(self, 'ORBIT_SEMIMAJOR_AXIS') 364 except ValueError: 365 raise NonNumberSatelliteParamError(self, 'ORBIT_SEMIMAJOR_AXIS') 366 367 try: 368 self.orbit_eccentricity = float(self.satParams['ORBIT_ECCENTRICITY']) 369 except KeyError: 370 raise MissingSatelliteParamError(self, 'ORBIT_ECCENTRICITY') 371 except ValueError: 372 raise NonNumberSatelliteParamError(self, 'ORBIT_SEMIMAJOR_AXIS') 373 374 try: 375 self.planet_mass = float(self.satParams['PLANET_MASS']) 376 except KeyError: 377 raise MissingSatelliteParamError(self, 'PLANET_MASS') 378 except ValueError: 379 raise NonNumberSatelliteParamError(self, 'PLANET_MASS') 380 381 try: 382 self.nsr_period = float(self.satParams['NSR_PERIOD']) 383 except KeyError: 384 raise MissingSatelliteParamError(self, 'NSR_PERIOD') 385 except ValueError: 386 raise NonNumberSatelliteParamError(self, 'NSR_PERIOD') 387 388 # Figure out how many layers are listed in the input file, and give the 389 # user an error if it's not exactly 4 (the current Love number code can 390 # only deal with 4): 391 self.num_layers = 0 392 for key in self.satParams.keys(): 393 if re.match('LAYER_ID', key): self.num_layers += 1 394 395 if self.num_layers != 4: 396 raise LoveLayerNumberError(self) 397 398 # construct the layer objects: 399 self.layers = [] 400 n = 0 401 while (n < self.num_layers): 402 self.layers.append(SatLayer(self,layer_n=n)) 403 404 # Verify that layers increase in density toward the center of the 405 # satellite. If they don't the satellite will be gravitationally 406 # unstable, and the Love number code will blow up: 407 if (n > 0): 408 if (self.layers[n].density > self.layers[n-1].density): 409 raise GravitationallyUnstableSatelliteError(self, n) 410 411 n = n+1 412 413 if self.planet_mass < 10*self.mass(): 414 raise ExcessiveSatelliteMassError(self) 415 416 if self.orbit_eccentricity > 0.25: 417 raise LargeEccentricityError(self) 418 419 if self.nsr_period < 0: 420 raise NegativeNSRPeriodError(self)
421 422 # end __init__ 423 424 # Methods defining derived quantities:
425 - def mass(self):
426 """Calculate the mass of the satellite. (the sum of the layer masses)""" 427 mass = 0.0 428 radius = 0.0 429 for layer in self.layers: 430 mass += ((4.0/3.0)*scipy.pi)*((radius+layer.thickness)**3 - (radius**3))*layer.density 431 radius += layer.thickness 432 return(mass)
433
434 - def radius(self):
435 """Calculate the radius of the satellite (the sum of the layer thicknesses).""" 436 radius = 0.0 437 for layer in self.layers: 438 radius += layer.thickness 439 return(radius)
440
441 - def density(self):
442 """Calculate the mean density of the satellite in [kg m^-3].""" 443 return(self.mass() / ((4.0/3.0)*scipy.pi*(self.radius()**3)))
444
445 - def surface_gravity(self):
446 """Calculate the satellite's surface gravitational acceleration in [m s^-2].""" 447 return(pc.G*self.mass() / (self.radius()**2))
448
449 - def orbit_period(self):
450 """Calculate the satellite's Keplerian orbital period in seconds.""" 451 return(2.0*scipy.pi*scipy.sqrt( (self.orbit_semimajor_axis**3) / (pc.G*self.planet_mass) ))
452
453 - def mean_motion(self):
454 """Calculate the orbital mean motion of the satellite [rad s^-1].""" 455 return(2.0*scipy.pi / (self.orbit_period()))
456 457 # end derived quantities 458
459 - def __str__(self): #
460 """Output a satellite definition file equivalent to the object.""" 461 462 myStr = """# 463 # Satellite system definition file for use with the Python SatStress package. 464 # All quantities are in SI (meters, kilograms, seconds) units. For more 465 # information, see: 466 # 467 # http://code.google.com/p/satstress 468 469 ############################################################ 470 # Basic Satellite parameters: 471 ############################################################ 472 473 # Satellite system name 474 SYSTEM_ID = %s 475 476 # System orbital parameters: 477 PLANET_MASS = %g 478 ORBIT_ECCENTRICITY = %g 479 ORBIT_SEMIMAJOR_AXIS = %g 480 481 # Additional parameters required to calculate Love numbers: 482 NSR_PERIOD = %g # seconds (== %g yr) 483 484 ############################################################ 485 # Derived properties of the Satellite: 486 ############################################################ 487 488 # RADIUS = %g 489 # MASS = %g 490 # DENSITY = %g 491 # SURFACE_GRAVITY = %g 492 # ORBIT_PERIOD = %g 493 # MEAN_MOTION = %g 494 """ % (self.system_id,\ 495 self.planet_mass,\ 496 self.orbit_eccentricity,\ 497 self.orbit_semimajor_axis,\ 498 self.nsr_period,\ 499 self.nsr_period/(31556926),\ 500 self.radius(),\ 501 self.mass(),\ 502 self.density(),\ 503 self.surface_gravity(),\ 504 self.orbit_period(),\ 505 self.mean_motion()) 506 507 myStr += """ 508 ############################################################ 509 # Layering structure of the Satellite: 510 ############################################################ 511 """ 512 n = 0 513 while (n < self.num_layers): 514 myStr += self.layers[n].ordered_str(n) 515 n += 1 516 517 return(myStr) 518 # end __str__ 519 # end class Satellite 520
521 -class SatLayer(object): #
522 """ 523 An object describing a uniform material layer within a satellite. 524 525 Note that a layer by itself has no knowledge of where within the satellite 526 it resides. That information is contained in the ordering of the list of 527 layers within the satellite object. 528 529 @ivar layer_id: a string identifying the layer, e.g. C{CORE}, or C{ICE_LOWER} 530 @type layer_id: str 531 @ivar density: the density of the layer, at zero pressure [kg m^-3], 532 @type density: float 533 @ivar lame_mu: the layer's Lame parameter, S{mu} (the shear modulus) [Pa]. 534 @type lame_mu: float 535 @ivar lame_lambda: the layer's Lame parameter, S{lambda} [Pa]. 536 @type lame_lambda: float 537 @ivar thickness: the radial thickness of the layer [m]. 538 @type thickness: float 539 @ivar viscosity: the viscosity of the layer [Pa s]. 540 @type viscosity: float 541 @ivar tensile_str: the tensile failure strength of the layer [Pa]. 542 @type tensile_str: float 543 """ 544 545
546 - def __init__(self, sat, layer_n=0): #
547 """ 548 Construct an object representing a layer within a L{Satellite}. 549 550 Gets values from the L{Satellite.satParams} dictionary for the 551 layer that corresponds to the value of C{layer_n}. 552 553 Each layer is defined by seven parameter values, and each layer has a 554 unique numeric identifier, appended to the end of all the names of its 555 parameters. Layer zero is the core, with the number increasing as the 556 satellite is built up toward the surface. In the below list the "N" at 557 the end of the parameter names should be replaced with the number of 558 the layer (C{layer_n}). Currently, because of the constraints of the 559 Love number code that we are using, you must specify 4 layers (C{CORE}, 560 C{OCEAN}, C{ICE_LOWER}, C{ICE_UPPER}). 561 562 - B{C{LAYER_ID_N}}: A string identifying the layer, e.g. C{OCEAN} or 563 C{ICE_LOWER} 564 565 - B{C{DENSITY_N}}: The density of the layer at zero pressure [m 566 kg^-3] 567 568 - B{C{LAME_MU_N}}: The real-valued Lame parameter S{mu} (shear 569 modulus) [Pa] 570 571 - B{C{LAME_LAMBDA_N}}: The real-valued Lame parameter S{lambda} [Pa] 572 573 - B{C{THICKNESS_N}}: The thickness of the layer [m] 574 575 - B{C{VISCOSITY_N}}: The viscosity of the layer [Pa s] 576 577 - B{C{TENSILE_STRENGTH_N}}: The tensile strength of the layer [Pa] 578 579 Not all layers necessarily require all parameters in order for the 580 calculation to succeed, but it is required that they be provided. 581 Parameters that will currently be ignored: 582 583 - B{C{TENSILE_STRENGTH}} of all layers except for the surface (which is 584 only used when creating fractures). 585 586 - B{C{VISCOSITY}} of the ocean and the core. 587 588 - B{C{LAME_MU}} of the ocean, assumed to be zero. 589 590 @param sat: the L{Satellite} object to which the layer belongs. 591 @type sat: L{Satellite} 592 @param layer_n: layer_n indicates which layer in the satellite is being 593 defined, with n=0 indicating the center of the satellite, and 594 increasing outward. This is needed in order to select the appropriate 595 values from the L{Satellite.satParams} dictionary. 596 @type layer_n: int 597 598 @raise MissingSatelliteParamError: if any of the seven input parameters 599 listed above is not found in the L{Satellite.satParams} dictionary 600 belonging to the C{sat} object passed in. 601 @raise NonNumberSatelliteParamError: if any numeric instance variable 602 list above has a value that cannot be converted to a float. 603 @raise LowLayerDensityError: if a layer is specififed with a 604 density of less than 100 [kg m^-3] 605 @raise LowLayerThicknessError: if a layer is specified with a thickness 606 of less than 100 [m]. 607 @raise NegativeLayerParamError: if either of the Lame parameters, the 608 viscosity, or the tensile strength of a layer is less than zero. 609 610 @return: a L{SatLayer} object having the properties specified in the 611 @rtype: L{SatLayer} 612 613 """ 614 615 # Assign all of the SatLayer's data attributes, converting those 616 # that should be numbers into floats, and making sure that if we 617 # didn't get a required parameter, an exception is raised. 618 try: 619 self.layer_id = sat.satParams['LAYER_ID_' + str(layer_n)] 620 except KeyError: 621 raise MissingSatelliteParamError(sat, 'LAYER_ID_' + str(layer_n)) 622 623 try: 624 self.density = float(sat.satParams['DENSITY_' + str(layer_n)]) 625 except KeyError: 626 raise MissingSatelliteParamError(sat, 'DENSITY_' + str(layer_n)) 627 except ValueError: 628 raise NonNumberSatelliteParamError(sat, 'DENSITY_' + str(layer_n)) 629 630 try: 631 self.lame_mu = float(sat.satParams['LAME_MU_' + str(layer_n)]) 632 except KeyError: 633 raise MissingSatelliteParamError(sat, 'LAME_MU_' + str(layer_n)) 634 except ValueError: 635 raise NonNumberSatelliteParamError(sat, 'LAME_MU_' + str(layer_n)) 636 637 try: 638 self.lame_lambda = float(sat.satParams['LAME_LAMBDA_' + str(layer_n)]) 639 except KeyError: 640 raise MissingSatelliteParamError(sat, 'LAME_LAMBDA_' + str(layer_n)) 641 except ValueError: 642 raise NonNumberSatelliteParamError(sat, 'LAME_LAMBDA' + str(layer_n)) 643 644 try: 645 self.thickness = float(sat.satParams['THICKNESS_' + str(layer_n)]) 646 except KeyError: 647 raise MissingSatelliteParamError(sat, 'THICKNESS_' + str(layer_n)) 648 except ValueError: 649 raise NonNumberSatelliteParamError(sat, 'THICKNESS_' + str(layer_n)) 650 651 try: 652 self.viscosity = float(sat.satParams['VISCOSITY_' + str(layer_n)]) 653 except KeyError: 654 raise MissingSatelliteParamError(sat, 'VISCOSITY_' + str(layer_n)) 655 except ValueError: 656 raise NonNumberSatelliteParamError(sat, 'VISCOSITY_' + str(layer_n)) 657 658 try: 659 self.tensile_str = float(sat.satParams['TENSILE_STR_' + str(layer_n)]) 660 except KeyError: 661 raise MissingSatelliteParamError(sat, 'TENSILE_STR_' + str(layer_n)) 662 except ValueError: 663 raise NonNumberSatelliteParamError(sat, 'TENSILE_STR_' + str(layer_n)) 664 665 # Test the SatLayer data attributes to make sure they have nominally 666 # reasonable values: 667 668 if self.density <= 100.0: 669 raise LowLayerDensityError(sat, layer_n) 670 671 if self.thickness <= 100.0: 672 raise LowLayerThicknessError(sat, layer_n) 673 674 if self.lame_mu < 0.0: 675 raise NegativeLayerParamError(sat, 'LAME_MU_' + str(layer_n)) 676 677 if self.lame_lambda < 0.0: 678 raise NegativeLayerParamError(sat, 'LAME_LAMBDA_' + str(layer_n)) 679 680 if self.viscosity < 0.0: 681 raise NegativeLayerParamError(sat, 'VISCOSITY_' + str(layer_n)) 682 683 if self.tensile_str < 0.0: 684 raise NegativeLayerParamError(sat, 'TENSILE_STRENGTH_' + str(layer_n)) 685 686 # end __init__() 687
688 - def __str__(self):
689 """ 690 Output a human and machine readable text description of the layer. 691 692 Note that because the layer object does not know explicitly where it is 693 within the stratified L{Satellite} (that information is contained in 694 the ordering of Satellite.layers list). 695 696 Derived quantities are also output for clarity, but are preceeded by 697 the comment character '#' to avoid accidental use in re-constituting a 698 satellite. 699 700 """ 701 702 return(self.ordered_str())
703 # end __str__() 704
705 - def ordered_str(self, layer_n=None):
706 """ 707 Return a string representation of the L{SatLayer}, including information 708 about where in the L{Satellite} the layer lies. 709 710 When being output as part of a L{Satellite} object, the L{SatLayer} needs 711 to communicate which layer it is. 712 713 @param layer_n: Layer location within satellite. Zero is the core. Defaults 714 to None, in which case no ordering information is conveyed in output. 715 @type layer_n: int 716 717 """ 718 if layer_n is not None: 719 order_str = "_"+str(layer_n) 720 else: 721 order_str = "" 722 723 myStr = """ 724 LAYER_ID%s = %s 725 DENSITY%s = %g 726 LAME_MU%s = %g 727 LAME_LAMBDA%s = %g 728 THICKNESS%s = %g 729 VISCOSITY%s = %g 730 TENSILE_STR%s = %g 731 # MAXWELL_TIME%s = %g 732 # BULK_MODULUS%s = %g 733 # YOUNGS_MODULUS%s = %g 734 # POISSONS_RATIO%s = %g 735 # P_WAVE_VELOCITY%s = %g 736 """ % (order_str, str(self.layer_id),\ 737 order_str, self.density,\ 738 order_str, self.lame_mu,\ 739 order_str, self.lame_lambda,\ 740 order_str, self.thickness,\ 741 order_str, self.viscosity,\ 742 order_str, self.tensile_str,\ 743 order_str, self.maxwell_time(),\ 744 order_str, self.bulk_modulus(),\ 745 order_str, self.youngs_modulus(),\ 746 order_str, self.poissons_ratio(),\ 747 order_str, self.p_wave_velocity()) 748 749 return(myStr)
750 751 # Methods for calculating quantities that depend on data attributes
752 - def maxwell_time(self):
753 """Calculate the Maxwell relaxation time of the layer [s] 754 (viscosity/lame_mu).""" 755 if(self.viscosity==0.0 or self.viscosity=="NONE" or self.lame_mu==0.0): 756 return(0.0) 757 else: 758 return(self.viscosity / self.lame_mu)
759
760 - def bulk_modulus(self):
761 """Calculate the bulk modulus (S{kappa}) of the layer [Pa].""" 762 return(self.lame_lambda + (2.0/3.0) * self.lame_mu)
763
764 - def youngs_modulus(self):
765 """Calculate the Young's modulus (E) of the layer [Pa].""" 766 return(self.lame_mu * (3.0*self.lame_lambda + 2.0*self.lame_mu) / (self.lame_lambda + self.lame_mu))
767
768 - def poissons_ratio(self):
769 """Calculate poisson's ratio (S{nu}) of the layer [Pa].""" 770 return(self.lame_lambda / (2.0 * (self.lame_lambda + self.lame_mu) ))
771
772 - def p_wave_velocity(self):
773 """Calculate the velocity of a compression wave in the layer [m s^-1]""" 774 return(scipy.sqrt(self.bulk_modulus()/self.density))
775 # end derived quantities 776 777 # end class SatLayer 778
779 -class LoveNum(object): #
780 """A container class for the complex Love numbers: h2, k2, and l2. 781 782 @ivar h2: the degree 2 complex, frequency dependent Love number h. 783 @type h2: complex 784 @ivar k2: the degree 2 complex, frequency dependent Love number k. 785 @type k2: complex 786 @ivar l2: the degree 2 complex, frequency dependent Love number l. 787 @type l2: complex 788 """ 789
790 - def __init__(self, h2_real, h2_imag, k2_real, k2_imag, l2_real, l2_imag):
791 """Using the real and imaginary parts, create complex values.""" 792 self.h2 = h2_real + h2_imag*1.0j 793 self.k2 = k2_real + k2_imag*1.0j 794 self.l2 = l2_real + l2_imag*1.0j
795
796 - def __str__(self):
797 """Return a human readable string representation of the Love numbers""" 798 return("h2 = %s\nk2 = %s\nl2 = %s" % (self.h2, self.k2, self.l2))
799 800 # end class LoveNum 801
802 -class StressDef(object): #
803 """A base class from which particular tidal stress field objects descend. 804 805 Different tidal forcings are specified as sub-classes of this superclass (one 806 for each separate forcing). 807 808 In the expressions of the stress fields, the time M{t} is specified in 809 seconds, with zero occuring at periapse, in order to be compatible with the 810 future inclusion of stressing mechanisms which may have explicit time 811 dependence instead of being a function of the satellite's orbital position 812 (e.g. a true polar wander trajectory). 813 814 Location is specified within a polar coordinate system having its origin at the 815 satellite's center of mass, using the following variables: 816 817 - co-latitude (S{theta}): The arc separating a point on the surface of 818 the satellite from the north pole (0 < S{theta} < S{pi}). 819 820 - longitude (S{phi}): The arc separating the meridian of a point and the 821 meridian which passes under the average location of the primary 822 (planet) in the sky over the course of an orbit (0 < S{phi} < 2S{pi}). 823 B{East is taken as positive.} 824 825 Each subclass must define its own version of the three components of the 826 membrane stress tensor, C{Ttt}, C{Tpp}, and C{Tpt} (the north-south, 827 east-west, and shear stress components) as methods. 828 829 @cvar satellite: the satellite which the stress is being applied to. 830 @type satellite: L{Satellite} 831 @cvar omega: the forcing frequency associated with the stress. 832 @type omega: float 833 @cvar love: the Love numbers which result from the given forcing frequency 834 and the specified satellite structure. 835 @type love: L{LoveNum} 836 837 """ 838 839 # These variables (objects) are used in the superclass and in the subclasses, 840 # and so need to be defined here. (Is this true?) 841 omega = 0.0 842 satellite = None 843 love = LoveNum(0,0,0,0,0,0) 844 845 # Common StressDef Methods:
846 - def __str__(self):
847 """ 848 Output information about this stress field, including frequency 849 dependent parameters. 850 851 """ 852 myStr = """ 853 %s_SURFACE_DELTA = %g 854 %s_OMEGA = %g 855 %s_FORCING_PERIOD = %g 856 %s_MU_TWIDDLE = %g + %gj 857 %s_LAMBDA_TWIDDLE = %g + %gj 858 %s_LOVE_H2 = %s 859 %s_LOVE_K2 = %s 860 %s_LOVE_L2 = %s 861 """ % (self.__name__, self.Delta(),\ 862 self.__name__, self.omega,\ 863 self.__name__, self.forcing_period(),\ 864 self.__name__, self.mu_twiddle().real, self.mu_twiddle().imag,\ 865 self.__name__, self.lambda_twiddle().real, self.lambda_twiddle().imag,\ 866 self.__name__, self.love.h2,\ 867 self.__name__, self.love.k2,\ 868 self.__name__, self.love.l2) 869 870 return(myStr)
871
872 - def forcing_period(self):
873 """ 874 Calculate the forcing period based on the forcing frequency (omega). 875 """ 876 return(2*scipy.pi/self.omega)
877
878 - def calcLove(self): #
879 """ 880 Calculate the Love numbers for the satellite and the given forcing. 881 882 If an infinite forcing period is given, return zero valued Love 883 numbers. 884 885 This is a wrapper function, which can be used to call different Love 886 number codes in the future. 887 888 @raise InvalidLoveNumberError: if the magnitude of the imaginary part 889 of any Love number is larger than its real part, if the real part is 890 ever less than zero, or if the real coefficient of the imaginary part 891 is ever positive. 892 893 """ 894 # I am separating the Love number calculation out into a separate 895 # method because I anticipate having more than one option here, 896 # possibly allowing the user to explicitly set their Love numbers in 897 # the input file, and possibly wrapping the Fortran code in Python, 898 # or translating James Roberts' code into pure python, so I don't have 899 # to worry about this thing compiling everywhere... 900 901 # If the period is infinite (omega = 0) then we aren't going to get any 902 # stresses from the forcing at all, and our Love number code will probably 903 # blow up anyway, so for the moment we'll just set the Love numbers to be 904 # null, and deal with them appropriately in the stress calculation elsewhere. 905 906 if self.omega == 0.0: 907 self.calcLoveInfinitePeriod() 908 else: 909 self.calcLoveWahr4LayerExternal() 910 911 # Do a little sanity checking on the Love numbers to make sure that we have 912 # realistic values: 913 # - Imaginary parts should always have smaller magnitudes than real parts 914 # - Real parts should be greater than zero 915 # - Imaginary parts should be less than zero 916 if (abs(self.love.h2.real) < abs(self.love.h2.imag)) or\ 917 (abs(self.love.l2.real) < abs(self.love.l2.imag)) or\ 918 (self.love.h2.real < 0) or\ 919 (self.love.l2.real < 0) or\ 920 ((self.love.h2.imag*-1j).real > 0) or\ 921 ((self.love.l2.imag*-1j).real > 0): 922 raise InvalidLoveNumberError(self) 923 924 # end calcLove() 925
926 - def calcLoveInfinitePeriod(self): #
927 """ 928 Return a set of zero Love numbers constructed statically. 929 930 This method is included so we don't have to worry about whether the 931 Love number code can deal with being given an infinite period. All 932 stresses will relax to zero with an infinite period (since the shear 933 modulus S{mu} goes to zero), so it doesn't really matter what we set 934 the Love numbers to here. 935 """ 936 # Yes, I know this actually corresponds to a rigid body. 937 self.love = LoveNum(0.0, 0.0, 0.0, 0.0, 0.0, 0.0) 938
939 - def calcLoveWahr4LayerExternal(self): #
940 """Use John Wahr's Love number code to calculate h, k, and l. 941 942 At the moment, the code is fairly limited in the kind of input it can 943 take. The specified satellite must: 944 945 - use a Maxwell rheology 946 947 - have a liquid water ocean underlying the ice shell 948 949 - have a 4-layer structure (ice_upper, ice_lower, ocean, core) 950 951 Eventually the Love number code will be more closely integrated with 952 this package, allowing more flexibility in the interior structure of 953 the satellite. 954 955 A temporary directory named lovetmp-XXXXXXX (where the X's are a random 956 hexadecimal number) is created in the current working directory, within 957 which the Love number code is run. The directory is deleted 958 immediately following the calculation. 959 960 @raise LoveExcessiveDeltaError: if L{StressDef.Delta}() > 10^9 for 961 either of the ice layers. 962 """ 963 # random number generator for creating random temporary filenames: 964 import random 965 # methods for interfacing with the operating system: 966 import os 967 968 # make sure that Delta isn't so big that the Love number code will fail. 969 # For John's Love number code, things work up to about Delta = 10^9 970 # Note that we need to check for both of the ice layers. 971 # look at the ice layers to see if they're too gooey for this Love number 972 # code to deal with. Negative array indices mean count from the END of the 973 # list, so, the surface of the satellite is -1, the next-to-surface layer 974 # is -2, etc. 975 for layer_n in (-1, -2): 976 if self.Delta(layer_n) > 1e9: 977 raise LoveExcessiveDeltaError(self, layer_n) 978 979 orig_dir = os.getcwd() 980 981 # write an input file that John's Love number code can read 982 # Because we're doing a bunch of stuff outside of this package's 983 # confines here... and we really have no idea what the OS is like 984 # we do this inside a try/except structure... and do our best to 985 # clean up after ourselves. 986 try: 987 lovetmp = "lovetmp-%s" % hex(int(random.random()*2**32))[2:-1] 988 os.mkdir(lovetmp) 989 990 # We need to be reasonably sure that the Love number code is 991 # actually in the current PATH. If the package has been installed 992 # in the normal way, (using setup.py) then we're fine (the Love 993 # number code ought to have been installed elsewhere on the 994 # system), but if we're in the code repository, or if we're just 995 # doing the tests pre-installation, then the code is in a funny 996 # spot - so let's add that spot to our PATH explicitly here. 997 # Unfortunately, this will only work on Unix machines... 998 dirname = os.path.dirname(os.path.abspath(__file__)) 999 os.environ['PATH'] += ":%s/Love/JohnWahr" % (dirname,) 1000 1001 os.chdir(lovetmp) 1002 love_infile = open("in.love", 'w') 1003 1004 love_infile.write("""%g Mean Density of Satellite (g/cm^3) 1005 1 Rheology (1=Maxwell, 0=elastic) 1006 %g Forcing period (earth days, 86400 seconds) 1007 %g Viscosity of upper ice layer (Pa sec) 1008 %g Viscosity of lower ice layer (Pa sec) 1009 %g Young's modulus for the rocky core 1010 %g Poisson's ratio for the rocky core 1011 %g Density of the rocky core (g/cm^3) 1012 %g Young's modulus for ice 1013 %g Poisson's ratio for ice 1014 %g Density of ice 1015 1 Decoupling fluid layer (e.g. global ocean)? (1=yes, 0=no) 1016 %g Thickness of fluid layer (km) 1017 %g Density of fluid layer (g/cm^3) 1018 %g P-wave velocity in fluid layer (km/sec) 1019 %g Total radius of satellite (km) 1020 %g Thickness of upper (cold) ice layer (km) 1021 %g Thickness of lower (warm) ice layer (km) 1022 165 Total number of calculation nodes 1023 3 Number of density discontinuities 1024 154 Node number of innermost boundary 1025 161 Node number of 2nd innermost boundary 1026 163 Node number of 3rd innermost boundary""" % (self.satellite.density()/1000.0,\ 1027 self.forcing_period()/86400,\ 1028 self.satellite.layers[3].viscosity,\ 1029 self.satellite.layers[2].viscosity,\ 1030 self.satellite.layers[0].youngs_modulus(),\ 1031 self.satellite.layers[0].poissons_ratio(),\ 1032 self.satellite.layers[0].density/1000.0,\ 1033 self.satellite.layers[3].youngs_modulus(),\ 1034 self.satellite.layers[3].poissons_ratio(),\ 1035 self.satellite.layers[3].density/1000.0,\ 1036 self.satellite.layers[1].thickness/1000.0,\ 1037 self.satellite.layers[1].density/1000.0,\ 1038 self.satellite.layers[1].p_wave_velocity(),\ 1039 self.satellite.radius()/1000.0,\ 1040 self.satellite.layers[3].thickness/1000.0,\ 1041 self.satellite.layers[2].thickness/1000.0) ) 1042 1043 love_infile.close() 1044 1045 # Actually call the code (it looks at in.love and writes to out.love 1046 # in our current working directory) 1047 os.popen('calcLoveWahr4Layer') 1048 1049 # read in the output file that it creates 1050 # take the last line and strip() off the newline at the end 1051 # then split() it by whitespace 1052 luv = open("out.love", 'r').readlines()[-1].strip().split() 1053 # extract those fields that correspond to the love numbers 1054 (h2_real, h2_imag, k2_real, k2_imag, l2_real, l2_imag) = luv[3:5]+luv[6:8]+luv[9:11] 1055 # generate a LoveNum Object made of our new love numbers, and return it. 1056 self.love = LoveNum(float(h2_real), float(h2_imag), float(k2_real), float(k2_imag), float(l2_real), float(l2_imag)) 1057 1058 except: 1059 # If the above attemped to run the love2 program did cause some 1060 # kind of error, we should just raise the exception up the chain. 1061 raise 1062 1063 finally: 1064 # Get back to our original working directory, cleanup after the 1065 # love2 program. We need to do these things even if our Love number 1066 # calculation failed somehow... 1067 os.chdir(orig_dir) 1068 try: 1069 os.remove(os.path.join(lovetmp, "in.love")) 1070 os.remove(os.path.join(lovetmp, "out.love")) 1071 os.remove(os.path.join(lovetmp, "out.model")) 1072 os.rmdir(lovetmp) 1073 except OSError: 1074 # If the Love number calculation DID fail, then the 1075 # above attempt to get rid of the files associated with it 1076 # will probably raise exceptions (directory not found, etc) 1077 # and that's okay - we just do nothing. 1078 pass 1079 1080 # end calcLoveJohnWahr4Layer() 1081
1082 - def Delta(self, layer_n=-1):
1083 """ 1084 Calculate S{Delta}, a measure of how viscous the layer's response is. 1085 1086 @param layer_n: indicates which satellite layer Delta should be 1087 calculated for, defaulting to the surface (recall that layer 0 is the 1088 core) 1089 @type layer_n: int 1090 @return: S{Delta}= S{mu}/(S{omega}*S{eta}) 1091 @rtype: float 1092 1093 """ 1094 return(self.satellite.layers[layer_n].lame_mu/(self.omega*self.satellite.layers[layer_n].viscosity))
1095
1096 - def Z(self):
1097 """ 1098 Calculate the value of Z, a constant that sits in front of many terms 1099 in the potential defined by Wahr et al. (2008). 1100 1101 @return: Z, a common constant in many of the Wahr et al. potential terms. 1102 @rtype: float 1103 """ 1104 numerator = 3.0*pc.G*self.satellite.planet_mass*self.satellite.radius()**2 1105 denominator = 2.0*self.satellite.orbit_semimajor_axis**3 1106 return (numerator/denominator)
1107
1108 - def mu_twiddle(self, layer_n=-1):
1109 """ 1110 Calculate the frequency-dependent Lame parameter S{mu} for a Maxwell 1111 rheology. 1112 1113 @param layer_n: number of layer for which we want to calculate S{mu}, defaults 1114 to the surface (with the core being layer zero). 1115 @return: the frequency-dependent Lame parameter S{mu} for a Maxwell rheology 1116 @rtype: complex 1117 """ 1118 1119 return(self.satellite.layers[layer_n].lame_mu*(1.0/(1-1.0j*self.Delta(layer_n))))
1120
1121 - def lambda_twiddle(self, layer_n=-1):
1122 """ 1123 Calculate the frequency-dependent Lame parameter S{lambda} for a 1124 Maxwell rheology. 1125 1126 @param layer_n: number of layer for which we want to calculate S{mu}, defaults 1127 to the surface (with the core being layer zero). 1128 @return: the frequency-dependent Lame parameter S{lambda} for a Maxwell 1129 rheology. 1130 @rtype: complex 1131 """ 1132 1133 lame_lambda = self.satellite.layers[layer_n].lame_lambda 1134 lame_mu = self.satellite.layers[layer_n].lame_mu 1135 1136 numerator = (1.0 - 1j*self.Delta(layer_n)*((2.0*lame_mu+3.0*lame_lambda)/(3.0*lame_lambda))) 1137 denominator = (1.0 - 1j*self.Delta(layer_n)) 1138 return(lame_lambda * (numerator/denominator))
1139
1140 - def alpha(self):
1141 """ 1142 Calculate the coefficient alpha twiddle for the surface layer (see 1143 Wahr et al. 2008). 1144 1145 @return: Calculate the coefficient alpha twiddle for the surface layer 1146 (see Wahr et al. 2008). 1147 @rtype: complex 1148 1149 """ 1150 numerator = 3.0*self.lambda_twiddle()+2.0*self.mu_twiddle() 1151 denominator = self.lambda_twiddle()+2.0*self.mu_twiddle() 1152 return(numerator/denominator)
1153
1154 - def Gamma(self):
1155 """ 1156 Calculate the coefficient capital Gamma twiddle for the surface layer 1157 (see Wahr et al. 2008). 1158 1159 @return: the coefficient capital Gamma twiddle for the surface layer 1160 (see Wahr et al. 2008). 1161 @rtype: complex 1162 """ 1163 return(self.mu_twiddle()*self.love.l2)
1164
1165 - def b1(self):
1166 """ 1167 Calculate the coefficient beta one twiddle for the surface layer (see Wahr et al. 2008). 1168 1169 @return: the coefficient beta one twiddle for the surface layer (see Wahr et al. 2008). 1170 @rtype: complex 1171 """ 1172 return(self.mu_twiddle()*(self.alpha()*(self.love.h2-3.0*self.love.l2)+3.0*self.love.l2))
1173
1174 - def g1(self):
1175 """ 1176 Calculate the coefficient gamma one twiddle for the surface layer (see Wahr et al. (2008)). 1177 1178 @return: the coefficient gamma one twiddle for the surface layer (see Wahr et al. (2008)). 1179 @rtype: complex 1180 """ 1181 return(self.mu_twiddle()*(self.alpha()*(self.love.h2-3.0*self.love.l2)-self.love.l2))
1182
1183 - def b2(self):
1184 """ 1185 Calculate the coefficient beta two twiddle for the surface layer (see Wahr et al. (2008)). 1186 1187 @return: the coefficient beta two twiddle for the surface layer (see Wahr et al. (2008)). 1188 @rtype: complex 1189 """ 1190 return(self.mu_twiddle()*(self.alpha()*(self.love.h2-3.0*self.love.l2)-3.0*self.love.l2))
1191
1192 - def g2(self):
1193 """ 1194 Calculate the coefficient gamma two twiddle for the surface layer (see Wahr et al. (2008)). 1195 @return: the coefficient gamma two twiddle for the surface layer (see Wahr et al. (2008)). 1196 @rtype: complex 1197 """ 1198 return(self.mu_twiddle()*(self.alpha()*(self.love.h2-3.0*self.love.l2)+self.love.l2))
1199 1200 # end Common StressDef Methods 1201
1202 - def Ttt(self, theta, phi, t):
1203 """ 1204 Calculates the S{tau}_S{theta}S{theta} (north-south) component of the 1205 stress tensor. 1206 1207 In the base class, this is a purely virtual method - it must be defined 1208 by the subclasses that describe particular tidal stresses. 1209 1210 @param theta: the co-latitude of the point at which to calculate the stress [rad]. 1211 @type theta: float 1212 @param phi: the east-positive longitude of the point at which to 1213 calculate the stress [rad]. 1214 @type phi: float 1215 @param t: the time, in seconds elapsed since pericenter, at which to perform the 1216 stress calculation [s]. 1217 @type t: float 1218 @return: the S{tau}_S{theta}S{theta} component of the 2x2 membrane stress tensor. 1219 @rtype: float 1220 1221 """ 1222 pass
1223
1224 - def Tpp(self, theta, phi, t):
1225 """ 1226 Calculates the S{tau}_S{phi}S{phi} (east-west) component of the stress tensor. 1227 1228 In the base class, this is a purely virtual method - it must be defined 1229 by the subclasses that describe particular tidal stresses. 1230 1231 @param theta: the co-latitude of the point at which to calculate the stress [rad]. 1232 @type theta: float 1233 @param phi: the east-positive longitude of the point at which to 1234 calculate the stress [rad]. 1235 @type phi: float 1236 @param t: the time, in seconds elapsed since pericenter, at which to perform the 1237 stress calculation [s]. 1238 @type t: float 1239 @return: the S{tau}_S{phi}S{phi} component of the 2x2 membrane stress tensor. 1240 @rtype: float 1241 1242 """ 1243 pass
1244
1245 - def Tpt(self, theta, phi, t):
1246 """ 1247 Calculates the S{tau}_S{phi}S{theta} (off-diagonal) component of the 1248 stress tensor. 1249 1250 In the base class, this is a purely virtual method - it must be defined 1251 by the subclasses that describe particular tidal stresses. 1252 1253 @param theta: the co-latitude of the point at which to calculate the stress [rad]. 1254 @type theta: float 1255 @param phi: the east-positive longitude of the point at which to 1256 calculate the stress [rad]. 1257 @type phi: float 1258 @param t: the time in seconds elapsed since pericenter, at which to perform the 1259 stress calculation [s]. 1260 @type t: float 1261 @return: the S{tau}_S{phi}S{theta} component of the 2x2 membrane stress tensor. 1262 @rtype: float 1263 """ 1264 pass
1265 1266 # end class StressDef 1267
1268 -class NSR(StressDef): #
1269 """An object defining the stress field which arises from the 1270 non-synchronous rotation (NSR) of a satellite's icy shell. 1271 1272 NSR is a subclass of L{StressDef}. See the derivation and detailed 1273 discussion of this stress field in in Wahr et al. (2008). 1274 """ 1275
1276 - def __init__(self, satellite):
1277 """Initialize the definition of the stresses due to NSR of the ice shell. 1278 1279 The forcing frequency S{omega} is the frequency with which a point on 1280 the surface passes through a single hemisphere, because the NSR stress 1281 field is degree 2 (that is, it's 2x the expected S{omega} from a full 1282 rotation) 1283 1284 Because the core is not subject to the NSR forcing (it remains tidally 1285 locked and synchronously rotating), all stresses within it are presumed 1286 to relax away, allowing it to deform into a tri-axial ellipsoid, with 1287 its long axis pointing toward the parent planet. In order to allow for 1288 this relaxation the shear modulus (S{mu}) of the core is set to an 1289 artificially low value for the purpose of the Love number calculation. 1290 This increases the magnitude of the radial deformation (and the Love 1291 number h2) significantly. See Wahr et al. (2008) for complete 1292 discussion. 1293 1294 @param satellite: the satellite to which the stress is being applied. 1295 @type satellite: L{Satellite} 1296 @return: an object defining the NSR stresses for a particular satellite. 1297 @rtype: L{NSR} 1298 1299 """ 1300 1301 self.__name__ = 'NSR' 1302 self.satellite = satellite 1303 # Set the core's shear modulus to be low, so we get a core that responds to 1304 # the NSR forcing nearly as a fluid: 1305 self.satellite.layers[0].lame_mu = self.satellite.layers[0].lame_mu/1000 1306 1307 # For NSR, the forcing period is only half the NSR period (or the forcing 1308 # frequency is twice the NSR frequency) because the shell goes through an 1309 # an entire oscillation in half a rotation. That's why this has 4*pi instead 1310 # of only 2*pi. 1311 self.omega = 4.0*scipy.pi/self.satellite.nsr_period 1312 self.calcLove() 1313 1314 # Restore the satellite's core to normal... 1315 self.satellite.layers[0].lame_mu = self.satellite.layers[0].lame_mu*1000
1316
1317 - def Ttt(self, theta, phi, t):
1318 """ 1319 Calculates the S{tau}_S{theta}S{theta} (north-south) component of the 1320 stress tensor. 1321 """ 1322 1323 if self.omega == 0: 1324 return(0.0) 1325 TttN = (self.b1()-self.g1()*scipy.cos(2.0*theta))*scipy.exp(1j*(2.0*phi+self.omega*t)) 1326 TttN = TttN.real 1327 TttN = TttN * (self.Z()/(2.0*self.satellite.surface_gravity()*self.satellite.radius())) 1328 return(TttN)
1329
1330 - def Tpp(self, theta, phi, t):
1331 """ 1332 Calculates the S{tau}_S{phi}S{phi} (east-west) component of the stress tensor. 1333 """ 1334 if self.omega == 0: 1335 return(0.0) 1336 TppN = (self.b2()-self.g2()*scipy.cos(2.0*theta))*scipy.exp(1j*(2.0*phi+self.omega*t)) 1337 TppN = TppN.real 1338 TppN = TppN * (self.Z()/(2.0*self.satellite.surface_gravity()*self.satellite.radius())) 1339 return(TppN)
1340
1341 - def Tpt(self, theta, phi, t):
1342 """ 1343 Calculates the S{tau}_S{phi}S{theta} (off-diagonal) component of the 1344 stress tensor. 1345 """ 1346 1347 if self.omega == 0: 1348 return(0.0) 1349 TptN = self.Gamma()*1j*scipy.exp(1j*(2.0*phi+self.omega*t))*scipy.cos(theta) 1350 TptN = TptN.real 1351 TptN = TptN * (2.0*self.Z()/(self.satellite.surface_gravity()*self.satellite.radius())) 1352 return(TptN)
1353 1354 # end class NSR 1355
1356 -class Diurnal(StressDef):
1357 """ 1358 An object defining the stress field that arises on a satellite due to an 1359 eccentric orbit. 1360 1361 Diurnal is a subclass of L{StressDef}. See the derivation and detailed 1362 discussion of this stress field in in Wahr et al. (2008). 1363 """ 1364
1365 - def __init__(self, satellite):
1366 """ 1367 Sets the object's satellite and omega attributes; calculates Love numbers. 1368 1369 @param satellite: the satellite to which the stress is being applied. 1370 @type satellite: L{Satellite} 1371 @return: an object defining the NSR stresses for a particular satellite. 1372 @rtype: L{NSR} 1373 """ 1374 1375 self.__name__ = 'Diurnal' 1376 self.satellite = satellite 1377 self.omega = self.satellite.mean_motion() 1378 self.calcLove()
1379
1380 - def Ttt(self, theta, phi, t):
1381 """ 1382 Calculates the S{tau}_S{theta}S{theta} (north-south) component of the 1383 stress tensor. 1384 """ 1385 1386 Ttt1 = 3.0*(self.b1()-self.g1()*scipy.cos(2.0*theta))*scipy.exp(1j*self.omega*t)*scipy.cos(2.0*phi) 1387 Ttt2 = -1.0*(self.b1()+3.0*self.g1()*scipy.cos(2.0*theta))*scipy.exp(1j*self.omega*t) 1388 Ttt3 = -4.0*(self.b1()-self.g1()*scipy.cos(2.0*theta))*1j*scipy.exp(1j*self.omega*t)*scipy.sin(2.0*phi) 1389 TttD = Ttt1 + Ttt2 + Ttt3 1390 TttD = TttD.real*self.satellite.orbit_eccentricity*self.Z()/(2.0*self.satellite.surface_gravity()*self.satellite.radius()) 1391 return(TttD)
1392
1393 - def Tpp(self, theta, phi, t):
1394 """ 1395 Calculates the S{tau}_S{phi}S{phi} (east-west) component of the stress tensor. 1396 """ 1397 Tpp1 = 3.0*(self.b2()-self.g2()*scipy.cos(2.0*theta))*scipy.exp(1j*self.omega*t)*scipy.cos(2.0*phi) 1398 Tpp2 = -1.0*(self.b2()+3.0*self.g2()*scipy.cos(2.0*theta))*scipy.exp(1j*self.omega*t) 1399 Tpp3 = -4.0*(self.b2()-self.g2()*scipy.cos(2.0*theta))*1j*scipy.exp(1j*self.omega*t)*scipy.sin(2.0*phi) 1400 TppD = Tpp1 + Tpp2 + Tpp3 1401 TppD = TppD.real*self.satellite.orbit_eccentricity*self.Z()/(2.0*self.satellite.surface_gravity()*self.satellite.radius()) 1402 return(TppD)
1403
1404 - def Tpt(self, theta, phi, t):
1405 """ 1406 Calculates the S{tau}_S{phi}S{theta} (off-diagonal) component of the 1407 stress tensor. 1408 """ 1409 Tpt1 = -4.0*self.Gamma()*1j*scipy.exp(1j*(self.omega*t))*scipy.cos(theta)*scipy.cos(2.0*phi) 1410 Tpt2 = -3.0*self.Gamma()*scipy.exp(1j*self.omega*t)*scipy.cos(theta)*scipy.sin(2.0*phi) 1411 TptD = Tpt1 + Tpt2 1412 TptD = TptD.real*2.0*self.satellite.orbit_eccentricity*self.Z()/(self.satellite.surface_gravity()*self.satellite.radius()) 1413 return(TptD)
1414 1415 # end class Diurnal 1416
1417 -class StressCalc(object): #
1418 """ 1419 An object which calculates the stresses on the surface of a L{Satellite} 1420 that result from one or more stress fields. 1421 1422 @ivar stresses: a list of L{StressDef} objects, corresponding to the stresses which are to 1423 be included in the calculations done by the L{StressCalc} object. 1424 @type stresses: list 1425 """ 1426
1427 - def __init__(self, stressdefs):
1428 """ 1429 Defines the list of stresses which are to be calculated at a given point. 1430 1431 @param stressdefs: a list of L{StressDef} objects, corresponding to the 1432 stresses which are to be included in the calculation. 1433 @type stressdefs: list 1434 @return: 1435 @rtype: L{StressCalc} 1436 """ 1437 1438 self.stresses = stressdefs
1439
1440 - def tensor(self, theta, phi, t): #
1441 """ 1442 Calculates surface stresses and returns them as a 2x2 stress tensor. 1443 @param theta: the co-latitude of the point at which to calculate the stress [rad]. 1444 @type theta: float 1445 @param phi: the east-positive longitude of the point at which to 1446 calculate the stress [rad]. 1447 @type phi: float 1448 @param t: the time in seconds elapsed since pericenter, at which to perform the 1449 stress calculation [s]. 1450 @type t: float 1451 @return: symmetric 2x2 surface (membrane) stress tensor S{tau} 1452 @rtype: Numpy.array 1453 """ 1454 1455 # This naming is just to make things a little more readable. 1456 tau = numpy.array([ [0.0, 0.0],\ 1457 [0.0, 0.0] ]) 1458 1459 # Build up each component of the stress tensor, adding the contributions 1460 # of each forcing listed in self.stresses one at a time: 1461 for stress in self.stresses: 1462 Ttt = stress.Ttt(theta, phi, t) 1463 Tpt = Ttp = stress.Tpt(theta, phi, t) 1464 Tpp = stress.Tpp(theta, phi, t) 1465 1466 tau += numpy.array([ [Ttt, Tpt],\ 1467 [Ttp, Tpp] ]) 1468 1469 return(tau) 1470 # end tensor 1471 1472 # end class StressCalc 1473 1474 # Exception classes (for error handling): 1475 # 1476 # We're being pretty strict here: if we get anything dubious, we just raise the 1477 # exception and allow it to kill the program. 1478 1479
1480 -class Error(Exception):
1481 """Base class for errors within the SatStress module.""" 1482 pass
1483
1484 -class NameValueFileError(Error): #
1485 """Base class for errors related to NAME=VALUE style input files.""" 1486
1487 -class NameValueFileParseError(NameValueFileError): #
1488 """Indicates a poorly formatted NAME=VALUE files."""
1489 - def __init__(self, nvf, line):
1490 """Stores the file and line that generated the parse error. 1491 1492 The file object (nvf) and the contents of the poorly formed line 1493 (badline) are stored within the exception, so we can print an error 1494 message with useful debugging information to the user.""" 1495 1496 self.nvf = nvf 1497 self.line = line
1498
1499 - def __str__(self):
1500 return(""" 1501 Poorly formatted NAME=VALUE pair found in the following file: 1502 1503 %s 1504 1505 Each non-comment line must consist of a non-whitespace string, followed by 1506 an '=', followed by another non-whitespace string. For example: 1507 1508 MY_PARAMETER = 12345.678e+09 1509 1510 Here's the line I got stuck on: 1511 1512 %s""" % (self.nvf.name, self.line))
1513 # end NameValueFileParseError 1514
1515 -class NameValueFileDuplicateNameError(NameValueFileError): #
1516 """Indicates multiple copies of the same name in an input file."""
1517 - def __init__(self, nvf, name):
1518 """Stores the file and the NAME that was found to be multiply defined. 1519 1520 NAME is the key, which has been found to be multiply defined in 1521 the input file, nvf.""" 1522 self.nvf = nvf 1523 self.name = name
1524
1525 - def __str__(self):
1526 return(""" 1527 The name %s was defined multiple times in the input file: 1528 1529 %s 1530 1531 Each NAME in a NAME=VALUE file must be defined only once, otherwise it's 1532 unclear which VALUE ought to be stored. 1533 """ % (self.name, self.nvf.name))
1534 # end NameValueFileDuplicateNameError 1535 1536 # end NameValueFileError classes 1537
1538 -class SatelliteParamError(Error): #
1539 """Indicates a problem with the Satellite initialization.""" 1540 pass 1541
1542 -class MissingSatelliteParamError(SatelliteParamError):
1543 """Indicates a required parameter was not found in the input file."""
1544 - def __init__(self, sat, missingname):
1545 self.sat = sat 1546 self.missingname = missingname
1547 - def __str__(self):
1548 return(""" 1549 The required parameter %s was not found in the satellite definition file: 1550 1551 %s 1552 """ % (self.missingname, self.sat.sourcefile.name))
1553
1554 -class InvalidSatelliteParamError(SatelliteParamError):
1555 """Raised when a required parameter is not found in the input file."""
1556 - def __init__(self, sat):
1557 """ 1558 Default initialization of an InvalidSatelliteParamError 1559 1560 Simply sets self.sat = sat (a Satellite object). Most errors can be well 1561 described using only the parameters stored in the Satellite object. 1562 """ 1563 self.sat = sat
1564
1565 -class LargeEccentricityError(InvalidSatelliteParamError):
1566 """Raised when satellite orbital eccentricity is > 0.25"""
1567 - def __str__(self):
1568 return(""" 1569 The specified orbital eccentricity (e = %g) found in the file 1570 1571 %s 1572 1573 is too big. The model implemented by this module, and described in Wahr et al. 1574 (2008) requires that orbital eccentricity be relatively small (< 0.25), for 1575 larger values of eccentricity, some of the mathematical approximations used 1576 break down. 1577 """ % (self.sat.orbit_eccentricity, self.sat.sourcefile.name))
1578
1579 -class NegativeNSRPeriodError(InvalidSatelliteParamError):
1580 """Raised if the satellite's NSR period is less than zero"""
1581 - def __str__(self):
1582 return(""" 1583 The input file: 1584 1585 %s 1586 1587 specifies an NSR_PERIOD < 0.0: 1588 1589 NSR_PERIOD = %g 1590 1591 but NSR_PERIOD can't be negative... 1592 """ % (self.sat.sourcefile.name, self.sat.nsr_period))
1593
1594 -class ExcessiveSatelliteMassError(InvalidSatelliteParamError):
1595 """ 1596 Raised if the satellite's parent planet is less than 10x as massive as the 1597 satellite.""" 1598
1599 - def __str__(self):
1600 return(""" 1601 The mass of the planet you specified: 1602 1603 PLANET_MASS = %g kg 1604 1605 is not that much larger than the mass of the satellite: 1606 1607 SATELLITE_MASS = %g kg 1608 1609 That seems unlikely to be correct. Maybe you have a units problem, or 1610 accidentally set the density of one of the layers in the satellite to be 1611 unrealistically high? 1612 """ % (self.sat.planet_mass, self.sat.mass()))
1613
1614 -class LoveLayerNumberError(InvalidSatelliteParamError):
1615 """ 1616 Raised if the number of layers specified in the satellite definition file 1617 is incompatible with the Love number code. 1618 """
1619 - def __str__(self):
1620 return(""" 1621 You specified a satellite with %d layers in the file: 1622 1623 %s 1624 1625 Unfortunately, at the moment the Love number code that the model uses (based on 1626 Dahlen, 1976) can only deal with a 4 layer satellite (core, ocean, and a 1627 2-layer ice shell). You can, however, set both of the ice layers to have 1628 identical properties, and vary the thicknesses of the layers significantly to 1629 minimize their effects if you want. 1630 """ % (self.sat.num_layers, self.sat.sourcefile.name))
1631
1632 -class InvalidLoveNumberError(InvalidSatelliteParamError):
1633 """Raised if the Love numbers are found to be suspicious."""
1634 - def __init__(self, stress, love):
1635 self.stress = stress
1636
1637 - def __str__(self):
1638 return(""" 1639 Sanity check failed on %s Love numbers: 1640 1641 %s 1642 1643 Real parts must have larger magnitudes than imaginary parts. Real parts must 1644 be positive. Imaginary parts must be negative. Satellite specification was 1645 read in from the following file: 1646 1647 %s 1648 1649 """ % (stress.__name__, stress.love, stress.satellite.sourcefile.name))
1650 -class LoveExcessiveDeltaError(InvalidSatelliteParamError):
1651 """Raised when S{Delta} > 10^9 for any of the ice layers, at which point 1652 the Love number code becomes unreliable."""
1653 - def __init__(self, stress, layer_n):
1654 self.stress = stress 1655 self.n = layer_n
1656
1657 - def __str__(self):
1658 return(""" 1659 The value of Delta (= lame_mu/(viscosity*omega), a measure of how viscously 1660 a layer responds to an applied forcing frequency) was found to be too large: 1661 1662 Delta = %g 1663 1664 in layer %d (LAYER_ID = %s), under the %s forcing. 1665 1666 specified in the input file: 1667 1668 %s 1669 1670 For the Love number code to produce reliable results, Delta should be less than 1671 about 10^9. Either reduce the forcing period or increase the viscosity. If 1672 you want an infinite forcing period (i.e. if you want all stresses due to the 1673 forcing to relax to zero) use: 1674 1675 FORCING_PERIOD = infinity 1676 """ % (self.stress.Delta(self.n),\ 1677 self.n,\ 1678 self.stress.satellite.layers[self.n].layer_id,\ 1679 self.stress.__name__,\ 1680 self.stress.satellite.sourcefile.name))
1681
1682 -class GravitationallyUnstableSatelliteError(InvalidSatelliteParamError):
1683 """ 1684 Raised if the density of layers is found not to decrease as you move toward the 1685 surface from the center of the satellite. 1686 """
1687 - def __init__(self, sat, layer_n):
1688 """Overrides the base InvalidSatelliteParamError.__init__() function. 1689 1690 We also need to know which layers have a gravitationally unstable arrangement. 1691 """ 1692 self.sat = sat 1693 self.n = layer_n
1694
1695 - def __str__(self):
1696 return(""" 1697 You have specified a satellite which is gravitationally unstable because 1698 one of the layers is on top of another having lower density: 1699 1700 Density of layer %d (%s) = %g [kg m^-3] 1701 Density of layer %d (%s) = %g [kg m^-3] 1702 1703 This will cause the Love number code to blow up. Maybe you have a typo 1704 in this input file: 1705 1706 %s 1707 1708 Or maybe you ordered the layers incorrectly. Layer 0 is the core, and the 1709 layer number increases with increasing distance from the core (it's like the 1710 radial dimension in spherical coordinates, not like "depth"). 1711 """ % (self.n, self.sat.layers[self.n].layer_id, self.sat.layers[self.n].density,\ 1712 self.n-1, self.sat.layers[self.n-1].layer_id, self.sat.layers[self.n-1].density,\ 1713 self.sat.sourcefile.name))
1714
1715 -class NonNumberSatelliteParamError(InvalidSatelliteParamError):
1716 """Indicates that a non-numeric value was found for a numerical parameter."""
1717 - def __init__(self, sat, badname):
1718 self.sat = sat 1719 self.badname = badname
1720 - def __str__(self):
1721 return(""" 1722 Found a non-numeric value for a numeric parameter: 1723 1724 %s = %s 1725 1726 in the input file: 1727 1728 %s 1729 """ % (self.badname, self.sat.satParams[self.badname], self.sat.sourcefile.name))
1730
1731 -class LowLayerDensityError(InvalidSatelliteParamError):
1732 """Indicates that a layer has been assigned an unrealistically low density."""
1733 - def __init__(self, sat, layer_n):
1734 self.sat = sat 1735 self.n = layer_n
1736
1737 - def __str__(self):
1738 return(""" 1739 Found an unrealistically low layer density: 1740 1741 DENSITY_%d = %g 1742 1743 specified in the input file: 1744 1745 %s 1746 1747 Recall that all input parameters are in SI (meters, kilograms, seconds) units, 1748 and so, for example, water has a density of 1000 [kg m^-3], not 1.0 [g cm^-3]. 1749 """ % (self.n, float(self.sat.satParams['DENSITY_'+str(self.n)]), self.sat.sourcefile.name))
1750
1751 -class LowLayerThicknessError(InvalidSatelliteParamError):
1752 """Indicates that a layer has been given too small of a thickness"""
1753 - def __init__(self, sat, layer_n):
1754 self.sat = sat 1755 self.n = layer_n
1756
1757 - def __str__(self):
1758 return(""" 1759 Found a layer thickness of less than 100 meters: 1760 1761 THICKNESS_%d = %g 1762 1763 specified in the input file: 1764 1765 %s 1766 1767 Recall that all input parameters are in SI (meters, kilograms, seconds) units. 1768 """ % (self.n, float(self.sat.satParams['THICKENSS_'+str(self.n)]), self.sat.sourcefile.name))
1769
1770 -class NegativeLayerParamError(InvalidSatelliteParamError):
1771 """Indicates a layer has been given an unphysical material property."""
1772 - def __init__(self, sat, badparam):
1773 self.sat = sat 1774 self.badparam = badparam
1775
1776 - def __str__(self):
1777 return(""" 1778 Found an unexpectedly negative value for %s: 1779 1780 %s = %g 1781 1782 specified in the input file: 1783 1784 %s 1785 """ % (self.badparam, self.badparam, float(self.sat.satParams[self.badparam]), self.sat.sourcefile.name))
1786