-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathscatparams.hpp
254 lines (203 loc) · 8.77 KB
/
scatparams.hpp
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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
// scatparams.hpp
//
// This file develops a class that represents scattering parameters in
// the manner used by Sato and Fehler and provides computation methods
// to get quantities derived from those parameters (E.g., Sato's G
// functions, X functions, and P functions)
//
//
#ifndef SCATPARAMS_H_
#define SCATPARAMS_H_
//
#include "geom.hpp" /* For S2:: and Geometry::Pi */
#include "elastic.hpp" /* Elastic::HetSpec and
Elastic::Velocity */
//////
// CLASS: ::: ScatterParams :::
// ENCAPS: Scattering Parameters
//
// Encapsulates the set of parameters that specify a particular
// scattering regime. This includes such parameters as the scale
// length and scattering strength, etc.
//
// USAGE PREREQUISITES:
//
// Class member cm_omega needs to be set via SetFrequencyHertz()
// before a ScatterParams object can be constructed from a HetSpec
// and a Velocity alone. If the user attempts otherwise, an
// exception of type std::logic_error will be thrown.
//
class ScatterParams {
private:
// :::::::::::::::::::::::::::::::::::::::::::::::
// ::: Member Variables (ScatterParams Class) :::
// :::::::::::::::::::::::::::::::::::::::::::::::
//
// Scattering Parameters:
//
// (Form borrowed from PSPhonon)
//
// (In the future we may parameterize
// differently, but for now just do what
// Shearer and Earle did.)
//
// These values describe the heterogeneity spectrum, and are
// provided to the constructor via a HetSpec object:
Real nu; // density vs. velocity pert. scaling (see 4.48)
Real eps; // RMS velocity perturbation (fractional)
Real a; // correlation distance
Real kappa; // von Karman parameter for PSDF
// Values for the following are computed during construction from
// knowlege of the background velocities (provided as an argument to
// the constructor):
Real el; // S-wave wavenumber (=om/beta0)
Real gam0; // Pvel/Svel (often assumed = sqrt(3) )
private:
// ::::::::::::::::::::::::::::::::::::::::::::::::::::
// ::: Class-Private Members (ScatterParams Class) :::
// ::::::::::::::::::::::::::::::::::::::::::::::::::::
static Real cm_omega; // Class needs to be informed of phonon
static bool cm_omega_known; // angular frequency omega before el can
// be computed. Ideally, this should be
// set prior to constructing any
// ScatterParams object.
public:
// :::::::::::::::::::::::::::::::::::::::::::::
// ::: Set-Methods for Class-Private Members :::
// :::::::::::::::::::::::::::::::::::::::::::::
static void SetFrequencyHertz(Real freq) { // Use this to set the
cm_omega = 2.0 * freq * Geometry::Pi; // cm_omega member
cm_omega_known = true;
}
public:
// ::::::::::::::::::::::::::::::::::::::::::::::::
// ::: Custom Data Types (ScatterParams Class) :::
// ::::::::::::::::::::::::::::::::::::::::::::::::
// Define a few one-element enums to use as type-locked
// symbols. These willserve as initializer flags to trigger
// special-purpose constructors:
//
// For example, use as:
//
// ScatterParams SP(ScatterParams::SATO_TEST);
//
// to initialize a ScatterParams object with
// special-case constructor call.
//
enum SATO_TEST_e {SATO_TEST};
enum SHEARER_TEST_e {SHEARER_TEST};
enum SP_FORWARD_200_e {SP_FORWARD_200}; // Forward Scattering with
// MFP_S =~ 200,
// MFP_P =~ 600.
// :::::::::::::::::::::::::::::::::::::::::::
// ::: Constructors (ScatterParams Class) :::
// :::::::::::::::::::::::::::::::::::::::::::
public:
// Constructor takes background velocities and a heterogeneity
// spectrum, and computes the parameters needed for Sato and Fehler
// formulation:
ScatterParams(Elastic::Velocity vpvs,
Elastic::HetSpec hs) :
nu ( hs.nu() ),
eps ( hs.eps() ),
a ( hs.a() ),
kappa ( hs.kappa() ),
el ( cm_omega/vpvs.Vs() ), // Depends on frequency
gam0 ( vpvs.Vp()/vpvs.Vs()) {
if (!cm_omega_known)
{AckSpit();} // Barf if cm_omega not yet set
}
// This alternate constructor takes the heterogeneity spectrum and
// then directly takes the el and gam0 paramters. User won't
// typically use this constructor for general earth modelling, but
// could be used for more diagnostic purposes or for exploring the
// parameter space.
ScatterParams(Elastic::HetSpec hs,
Real el, Real gam0) :
nu ( hs.nu() ),
eps ( hs.eps() ),
a ( hs.a() ),
kappa ( hs.kappa() ),
el ( el ),
gam0 ( gam0 )
{}
// These ones represent special test cases:
ScatterParams(SATO_TEST_e dummy) : // Params intended to replicate
nu(0.8), eps(1.0), a(0.1), // figure 4.8 in Sato & Fehler.
kappa(0.5), el(1.0), gam0(1.7321)
{}
ScatterParams(SHEARER_TEST_e dummy) : // Params intended to be
nu(0.8), eps(1.0), a(1.0), // similar to those used by
kappa(0.5), el(1.0), // Shearer and Earle, 2004
gam0(1.7321) // (TODO: Find out what those
{} // parameters would be - these
// are just placeholder
// numbers.)
ScatterParams(SP_FORWARD_200_e dummy) : // Results in strong
nu(0.8), eps(0.01), a(4.0), // forward scattering and
kappa(0.8), el(2.17411), // MFP's (P,S) of ~600,
gam0(1.7301) // ~200.
{}
public:
// ::::::::::::::::::::::::::::::::::::::::::::
// ::: Member Access (ScatterParams Class) :::
// ::::::::::::::::::::::::::::::::::::::::::::
Real GetNu() const {return nu;}
Real GetEps() const {return eps;}
Real GetA() const {return a;}
Real GetKappa() const {return kappa;}
Real GetL() const {return el;}
Real GetGam0() const {return gam0;}
// :::::::::::::::::::::::::::::::::::::::::::::::::
// ::: Comparison Methods (ScatterParams Class) :::
// :::::::::::::::::::::::::::::::::::::::::::::::::
//
// Provides metrics for comparing one set
// of scatterparams to another, to
// facilitate parsimony in allocating the
// relatively memory intensive Scatterer
// objects at model build time.
//
Real CompareRoughly(const ScatterParams & other) const;
// Provides a VERY rough comparison between two ScatterParams
// which should be used only for the first design iteration of
// optimizing Scatterer allocation. The comparison is
// unsophisticated in that it only looks at the squared magnitude
// change in values, and sums them to an aggregate (possibly a
// weighted sum), but has no particular smarts about the affect
// each parameter change has on the actual resultant scattering
// behavior. This, it might in practice be over/under sensitive
// to particualar parameters, which may have performance or
// realism implications.
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ::: Computational Result Methods (ScatterParams Class) ::
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// These constant methods return
// functional mathematical results
// based on the parameters stored
// in this class.
//
// They do not change the internal
// state of the object.
//
void GSATO (S2::S2Point toa, // (A`la PSPhonon)
Real & gpp, Real & gps,
Real & gsp, Real & gss,
Real & spol) const;
void XSATO (S2::S2Point toa, // (A`la PSPhonon)
Real & xpp, Real & xps,
Real & xsp, Real & xss_psi,
Real & xss_zeta) const;
Real PSATO (Real m) const; // Sato P function, kinda like
// EXPSATO in PSPhonon, but we
// use von Karman autocor
private:
// :::::::::::::::::::::::::::::::::::::::::::::
// ::: Error Handling (ScatterParams Class) :::
// :::::::::::::::::::::::::::::::::::::::::::::
void AckSpit(); // Register a complaint if class used improperly
}; // class ScatterParams
////
///
#endif //#ifndef SCATPARAMS_H_
//