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
159 import re
160
161
162 import sys
163
164
165 import scipy
166 import numpy
167
168
169 import physcon as pc
170
171
172
173
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
213 line = re.sub("\s*%s.*" % (comment), "", line)
214
215
216 if re.match("^\s*$", line):
217 continue
218
219
220 if not (re.match("[^\s]*.*=.*[^\s]", line)):
221 raise NameValueFileParseError(nvf, line)
222
223 try:
224 (name,value) = re.split("=",line,1)
225 except ValueError:
226 raise NameValueFileParseError(nvf, line)
227
228
229 name = name.strip()
230 value = value.strip()
231
232
233
234
235 if params.has_key(name):
236 raise NameValueFileDuplicateNameError(nvf, name)
237 else:
238 params[name] = value
239
240 return params
241
242
243
244
245
246
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
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
343 self.sourcefile = satFile
344
345
346 try:
347 self.satParams = nvf2dict(self.sourcefile, comment='#')
348 except NameValueFileError:
349 raise
350
351
352
353
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
389
390
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
399 self.layers = []
400 n = 0
401 while (n < self.num_layers):
402 self.layers.append(SatLayer(self,layer_n=n))
403
404
405
406
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
423
424
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
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
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
446 """Calculate the satellite's surface gravitational acceleration in [m s^-2]."""
447 return(pc.G*self.mass() / (self.radius()**2))
448
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
454 """Calculate the orbital mean motion of the satellite [rad s^-1]."""
455 return(2.0*scipy.pi / (self.orbit_period()))
456
457
458
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
519
520
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
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
616
617
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
666
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
687
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
704
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
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
761 """Calculate the bulk modulus (S{kappa}) of the layer [Pa]."""
762 return(self.lame_lambda + (2.0/3.0) * self.lame_mu)
763
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
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
773 """Calculate the velocity of a compression wave in the layer [m s^-1]"""
774 return(scipy.sqrt(self.bulk_modulus()/self.density))
775
776
777
778
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
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
801
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
840
841 omega = 0.0
842 satellite = None
843 love = LoveNum(0,0,0,0,0,0)
844
845
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
873 """
874 Calculate the forcing period based on the forcing frequency (omega).
875 """
876 return(2*scipy.pi/self.omega)
877
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
895
896
897
898
899
900
901
902
903
904
905
906 if self.omega == 0.0:
907 self.calcLoveInfinitePeriod()
908 else:
909 self.calcLoveWahr4LayerExternal()
910
911
912
913
914
915
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
925
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
937 self.love = LoveNum(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
938
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
964 import random
965
966 import os
967
968
969
970
971
972
973
974
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
982
983
984
985
986 try:
987 lovetmp = "lovetmp-%s" % hex(int(random.random()*2**32))[2:-1]
988 os.mkdir(lovetmp)
989
990
991
992
993
994
995
996
997
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
1046
1047 os.popen('calcLoveWahr4Layer')
1048
1049
1050
1051
1052 luv = open("out.love", 'r').readlines()[-1].strip().split()
1053
1054 (h2_real, h2_imag, k2_real, k2_imag, l2_real, l2_imag) = luv[3:5]+luv[6:8]+luv[9:11]
1055
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
1060
1061 raise
1062
1063 finally:
1064
1065
1066
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
1075
1076
1077
1078 pass
1079
1080
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
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
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
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
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
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
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
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
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
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
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
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
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
1304
1305 self.satellite.layers[0].lame_mu = self.satellite.layers[0].lame_mu/1000
1306
1307
1308
1309
1310
1311 self.omega = 4.0*scipy.pi/self.satellite.nsr_period
1312 self.calcLove()
1313
1314
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
1355
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
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
1416
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
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
1456 tau = numpy.array([ [0.0, 0.0],\
1457 [0.0, 0.0] ])
1458
1459
1460
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
1471
1472
1473
1474
1475
1476
1477
1478
1479
1481 """Base class for errors within the SatStress module."""
1482 pass
1483
1485 """Base class for errors related to NAME=VALUE style input files."""
1486
1488 """Indicates a poorly formatted NAME=VALUE files."""
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
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
1514
1516 """Indicates multiple copies of the same name in an input file."""
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
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
1535
1536
1537
1539 """Indicates a problem with the Satellite initialization."""
1540 pass
1541
1543 """Indicates a required parameter was not found in the input file."""
1545 self.sat = sat
1546 self.missingname = missingname
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
1555 """Raised when a required parameter is not found in the input file."""
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
1566 """Raised when satellite orbital eccentricity is > 0.25"""
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
1580 """Raised if the satellite's NSR period is less than zero"""
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
1595 """
1596 Raised if the satellite's parent planet is less than 10x as massive as the
1597 satellite."""
1598
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
1615 """
1616 Raised if the number of layers specified in the satellite definition file
1617 is incompatible with the Love number code.
1618 """
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
1633 """Raised if the Love numbers are found to be suspicious."""
1635 self.stress = stress
1636
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))
1651 """Raised when S{Delta} > 10^9 for any of the ice layers, at which point
1652 the Love number code becomes unreliable."""
1654 self.stress = stress
1655 self.n = layer_n
1656
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
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 """
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
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
1716 """Indicates that a non-numeric value was found for a numerical parameter."""
1718 self.sat = sat
1719 self.badname = badname
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
1732 """Indicates that a layer has been assigned an unrealistically low density."""
1734 self.sat = sat
1735 self.n = layer_n
1736
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
1752 """Indicates that a layer has been given too small of a thickness"""
1754 self.sat = sat
1755 self.n = layer_n
1756
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
1771 """Indicates a layer has been given an unphysical material property."""
1773 self.sat = sat
1774 self.badparam = badparam
1775
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