1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

import numpy as num 

 

from pyrocko.guts import Float 

import pyrocko.orthodrome as od 

 

from . import compound_engine as ce 

from .base import SandboxSource, SourceProcessor 

 

 

d2r = num.pi / 180. 

r2d = 180. / num.pi 

km = 1e3 

 

 

class EllipsoidSource(SandboxSource): 

 

__implements__ = 'CompoundModel' 

 

length_x = Float.T( 

help='Length of x semi-axis of ellisoid before rotation in [m]', 

default=1.) 

length_y = Float.T( 

help='Length of y semi-axis of ellisoid before rotation in [m]', 

default=1.) 

length_z = Float.T( 

help='Length of z semi-axis of ellisoid before rotation in [m]', 

default=1.) 

rotation_x = Float.T( 

help='Clockwise rotation of ellipsoid around x-axis in [deg]', 

default=0.) 

rotation_y = Float.T( 

help='Clockwise rotation of ellipsoid around y-axis in [deg]', 

default=0.) 

rotation_z = Float.T( 

help='Clockwise rotation of ellipsoid around z-axis in [deg]', 

default=0.) 

 

mu = Float.T( 

help='Shear modulus, 2. Lame constant in [GPa]', 

default=8.,) 

lamda = Float.T( 

help='Lame constant in [GPa]', 

default=8.) 

cavity_pressure = Float.T( 

help='Pressure on the cavity walls in [GPa]', 

default=.5) 

 

@property 

def volume(self): 

K = self.lamda + 2 * self.mu / 3 # Bulk Modulus 

V = 4./3 * num.pi * self.length_x * self.length_y * self.length_z 

V = (self.cavity_pressure * V) / K 

return V 

 

def ECMParameters(self): 

if self._sandbox: 

east_shift, north_shift = od.latlon_to_ne_numpy( 

self._sandbox.frame.llLat, self._sandbox.frame.llLon, 

self.lat, self.lon) 

else: 

east_shift = 0. 

north_shift = 0. 

 

params = { 

'x0': self.easting + east_shift, 

'y0': self.northing + north_shift, 

'z0': self.depth, 

'rotx': self.rotation_x, 

'roty': self.rotation_y, 

'rotz': self.rotation_z, 

'ax': self.length_x, 

'ay': self.length_y, 

'az': self.length_z, 

'P': self.cavity_pressure * 1e9, 

'mu': self.mu * 1e9, 

'lamda': self.lamda * 1e9 

} 

return params 

 

 

class PointCompoundSource(SandboxSource): 

 

__implements__ = 'CompoundModel' 

 

rotation_x = Float.T( 

help='Clockwise rotation of ellipsoid around x-axis in [deg]', 

default=0.) 

rotation_y = Float.T( 

help='Clockwise rotation of ellipsoid around y-axis in [deg]', 

default=0.) 

rotation_z = Float.T( 

help='Clockwise rotation of ellipsoid around z-axis in [deg]', 

default=0.) 

dVx = Float.T( 

help='Volume change in x-plane in [m3]', 

default=1.) 

dVy = Float.T( 

help='Volume change in y-plane in [m3]', 

default=1.) 

dVz = Float.T( 

help='Volume change in z-plane in [m3]', 

default=1.) 

nu = Float.T( 

help='Poisson\'s ratio, typically 0.25', 

default=0.25) 

 

@property 

def volume(self): 

# After Nikkhoo, M et al. 2017 

return (self.dVx + self.dVy + self.dVz) / 1.8 

 

def pointCDMParameters(self): 

if self._sandbox: 

east_shift, north_shift = od.latlon_to_ne_numpy( 

self._sandbox.frame.llLat, self._sandbox.frame.llLon, 

self.lat, self.lon) 

else: 

east_shift = 0. 

north_shift = 0. 

 

params = { 

'x0': self.easting + east_shift, 

'y0': self.northing + north_shift, 

'z0': self.depth, 

'rotx': self.rotation_x, 

'roty': self.rotation_y, 

'rotz': self.rotation_z, 

'dVx': self.dVx**3, 

'dVy': self.dVy**3, 

'dVz': self.dVz**3, 

'nu': self.nu 

} 

return params 

 

 

class CompoundModelProcessor(SourceProcessor): 

 

__implements__ = 'CompoundModel' 

 

def process(self, sources, sandbox, nthreads=0): 

result = { 

'processor_profile': dict(), 

'displacement.e': num.zeros(sandbox.frame.npixel), 

'displacement.n': num.zeros(sandbox.frame.npixel), 

'displacement.d': num.zeros(sandbox.frame.npixel), 

} 

 

coords = sandbox.frame.coordinatesMeter 

 

for src in sources: 

if isinstance(src, EllipsoidSource): 

res = ce.ECM(coords, **src.ECMParameters()) 

elif isinstance(src, PointCompoundSource): 

res = ce.pointCDM(coords, **src.pointCDMParameters()) 

else: 

raise AttributeError('Source of wrong type!') 

result['displacement.e'] += res[0] 

result['displacement.n'] += res[1] 

result['displacement.d'] += res[2] 

 

return result