1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6from pyrocko.gui.snuffling import Snuffling, Param, Switch, Choice 

7from pyrocko.gui.marker import PhaseMarker 

8from pyrocko import gf 

9from pyrocko import cake 

10import numpy as num 

11 

12 

13class CakePhase(Snuffling): 

14 ''' 

15 <html> 

16 <head> 

17 <style type="text/css"> 

18 body { margin-left:10px }; 

19 </style> 

20 </head> 

21 <body> 

22 <h1 align="center">Theoretical Phase Arrivals</h1> 

23 <p> 

24 This snuffling uses pyrocko's 

25 <a href="http://emolch.github.io/pyrocko/v0.3/cake_doc.html">Cake</a> 

26 module to calculate seismic rays for layered earth models. </p> 

27 <p> 

28 <b>Parameters:</b><br /> 

29 <b>&middot; Global shift</b> - Add time onset to phases. <br /> 

30 <b>&middot; Add Model</b> - Add a model to drop down menu. <br /> 

31 <b>&middot; Add Phase</b> - Add a phase definition. 

32 (GUI reset required)<br /> 

33 </p> 

34 <p> 

35 Instructions and information on Cake's syntax of seismic rays can be 

36 found in the 

37 <a href="http://emolch.github.io/pyrocko/ 

38 v0.3/cake_doc.html#cmdoption-cake--phase">Cake documentation</a>. 

39 </p> 

40 </body> 

41 </html> 

42 ''' 

43 

44 def setup(self): 

45 self.set_name('Cake Phase') 

46 

47 # self._phase_names = ('PmP ~S ~P ~P(moho)s ~P(sill-top)s' 

48 # ' ~P(sill-bottom)s Pdiff').split() 

49 self._phase_names = ('~P Pg Sg pP p P Pdiff PKP PcP PcS PKIKP pPKIKP' 

50 ' SSP PPS SPP PSP SP PS ~PS ~SP Pn s S Sn PP PPP' 

51 ' ScS Sdiff SS SSS SKS SKIKS').split() 

52 

53 for iphase, name in enumerate(self._phase_names): 

54 self.add_parameter(Switch(name, 'wantphase_%i' % iphase, 

55 iphase == 0)) 

56 

57 model_names = cake.builtin_models() 

58 model_names = [ 

59 'Cake builtin: %s' % model_name for model_name in model_names] 

60 

61 self._engine = gf.LocalEngine(use_config=True) 

62 store_ids = self._engine.get_store_ids() 

63 

64 for store_id in store_ids: 

65 model_names.append('GF Store: %s' % store_id) 

66 

67 self._models = model_names 

68 

69 self.model_choice = Choice('Model', 'chosen_model', 

70 'ak135-f-continental.m', self._models) 

71 

72 self.add_parameter(self.model_choice) 

73 

74 self.add_parameter(Param('Global shift', 'tshift', 0., -20., 20.)) 

75 self.add_parameter(Switch('Use station depth', 

76 'use_station_depth', False)) 

77 self.add_trigger('Add Phase', self.add_phase_definition) 

78 self.add_trigger('Add Model', self.add_model_to_choice) 

79 self.add_trigger('Plot Model', self.plot_model) 

80 self.add_trigger('Plot Rays', self.plot_rays) 

81 

82 self._phases = {} 

83 self._model = None 

84 

85 def panel_visibility_changed(self, bool): 

86 pass 

87 

88 def wanted_phases(self): 

89 try: 

90 wanted = [] 

91 for iphase, name in enumerate(self._phase_names): 

92 if getattr(self, 'wantphase_%i' % iphase): 

93 if name in self._phases: 

94 phases = self._phases[name] 

95 else: 

96 if name.startswith('~'): 

97 phases = [cake.PhaseDef(name[1:])] 

98 else: 

99 phases = cake.PhaseDef.classic(name) 

100 

101 self._phases[name] = phases 

102 for pha in phases: 

103 pha.name = name 

104 

105 wanted.extend(phases) 

106 except (cake.UnknownClassicPhase, cake.PhaseDefParseError) as e: 

107 self.fail(str(e)) 

108 

109 return wanted 

110 

111 def call(self, plot_rays=False): 

112 

113 self.cleanup() 

114 wanted = self.wanted_phases() 

115 

116 if not wanted: 

117 return 

118 

119 event, stations = self.get_active_event_and_stations() 

120 

121 if not stations: 

122 self.fail('No station information available.') 

123 

124 self.update_model() 

125 model = self._model[1] 

126 

127 depth = event.depth 

128 if depth is None: 

129 depth = 0.0 

130 

131 allrays = [] 

132 alldists = [] 

133 for station in stations: 

134 dist = event.distance_to(station) 

135 alldists.append(dist) 

136 

137 if self.use_station_depth: 

138 rdepth = station.depth 

139 else: 

140 rdepth = 0.0 

141 

142 multi_dists = [] 

143 nmax = 1 

144 for i in range(0, nmax): 

145 multi_dists.append(dist*cake.m2d + 360.*i) 

146 multi_dists.append((i+1)*360. - dist*cake.m2d) 

147 

148 rays = model.arrivals( 

149 phases=wanted, 

150 distances=multi_dists, 

151 zstart=depth, 

152 zstop=rdepth) 

153 

154 for ray in rays: 

155 time = ray.t 

156 name = ray.given_phase().name 

157 incidence_angle = ray.incidence_angle() 

158 takeoff_angle = ray.takeoff_angle() 

159 

160 time += event.time + self.tshift 

161 m = PhaseMarker([(station.network, station.station, '*', '*')], 

162 time, time, 2, 

163 phasename=name, 

164 event=event, 

165 incidence_angle=incidence_angle, 

166 takeoff_angle=takeoff_angle) 

167 self.add_marker(m) 

168 

169 allrays.extend(rays) 

170 

171 if plot_rays: 

172 fig = self.figure(name='Ray Paths') 

173 from pyrocko import cake_plot 

174 cake_plot.my_rays_plot(model, None, allrays, depth, 0.0, 

175 num.array(alldists)*cake.m2d, 

176 axes=fig.gca()) 

177 

178 fig.canvas.draw() 

179 

180 def update_model(self): 

181 if not self._model or self._model[0] != self.chosen_model: 

182 if self.chosen_model.startswith('Cake builtin: '): 

183 load_model = cake.load_model( 

184 self.chosen_model.split(': ', 1)[1]) 

185 

186 elif self.chosen_model.startswith('GF Store: '): 

187 store_id = self.chosen_model.split(': ', 1)[1] 

188 

189 load_model = self._engine.get_store(store_id)\ 

190 .config.earthmodel_1d 

191 else: 

192 load_model = cake.load_model(self.chosen_model) 

193 

194 self._model = (self.chosen_model, load_model) 

195 

196 def update_model_choices(self): 

197 self.set_parameter_choices('chosen_model', self._models) 

198 

199 def add_model_to_choice(self): 

200 ''' 

201 Called from trigger 'Add Model'. 

202 

203 Adds another choice to the drop down 'Model' menu. 

204 ''' 

205 

206 in_model = self.input_filename('Load Model') 

207 if in_model not in self._models: 

208 self._models.append(in_model) 

209 self.update_model_choices() 

210 

211 self.set_parameter('chosen_model', in_model) 

212 self.call() 

213 

214 def add_phase_definition(self): 

215 ''' 

216 Called from trigger 'Add Phase Definition'. 

217 

218 Adds another phase option. 

219 Requires a reset of the GUI. 

220 ''' 

221 phase_def = str(self.input_dialog('Add New Phase', 

222 'Enter Phase Definition')) 

223 self._phase_names.append(phase_def) 

224 

225 self.add_parameter( 

226 Switch(phase_def, 

227 'wantphase_%s' % str(len(self._phase_names)-1), True)) 

228 self.reset_gui(reloaded=True) 

229 self.call() 

230 

231 def plot_model(self): 

232 self.update_model() 

233 

234 from pyrocko import cake_plot 

235 

236 fig = self.figure(name='Model: %s' % self._model[0]) 

237 

238 cake_plot.my_model_plot(self._model[1], axes=fig.gca()) 

239 

240 fig.canvas.draw() 

241 

242 def plot_rays(self): 

243 self.call(plot_rays=True) 

244 

245 

246def __snufflings__(): 

247 ''' 

248 Returns a list of snufflings to be exported by this module. 

249 ''' 

250 

251 return [CakePhase()]