-
Notifications
You must be signed in to change notification settings - Fork 103
/
batohilm.html
115 lines (104 loc) · 10.1 KB
/
batohilm.html
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
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN"
"http://www.w3.org/TR/html4/strict.dtd">
<html>
<head>
<title>SHTOOLS - Gravity and magnetics routines</title>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<link rel="stylesheet" type="text/css" href="../../css/sh.css">
<link rel="icon" type="image/vnd.microsoft.icon" href="../../images/favicon.ico">
</head>
<body>
<div class="main">
<p class="centeredimage"><img src="../../images/logo.jpg" width=894 height=135 alt="SHTOOLS --- Tools for working with spherical harmonics"></p>
<table class="menu">
<tbody>
<tr>
<td><a href="../../../index.html">HOME</a></td>
<td><a href="https://github.com/SHTOOLS/SHTOOLS/releases">DOWNLOAD</a></td>
<td class="selected"><a href="../../documentation.html">DOCUMENTATION</a></td>
<td><a href="../../faq.html">FAQ</a> </td>
</tr>
</tbody>
</table>
<p class="dir">
> <a href="../../../index.html" class="dir">Home</a> > <a href="../../documentation.html" class="dir">Documentation</a> > <a href="../../f95-routines.html" class="dir">Fortran 95</a> > <a href="../../gravmag.html" class="dir">Gravity and Magnetics</a></p>
<h1 id="batohilm">BAtoHilm</h1>
<p>Calculate iteratively the relief along an interface of constant density contrast that corresponds to a given Bouguer anomaly using the algorithm of Wieczorek and Phillips (1998).</p>
<h1 id="usage">Usage</h1>
<p>call BAtoHilm (<code>cilm</code>, <code>ba</code>, <code>grid</code>, <code>lmax</code>, <code>nmax</code>, <code>mass</code>, <code>r0</code>, <code>rho</code>, <code>gridtype</code>, <code>w</code>, <code>plx</code>, <code>zero</code>, <code>filtertype</code>, <code>filterdeg</code>, <code>lmaxcalc</code>, <code>exitstatus</code>)</p>
<h1 id="parameters">Parameters</h1>
<dl>
<dt><code>cilm</code> : output, real*8, dimension (2, <code>lmaxcalc</code>+1, <code>lmaxcalc</code>+1)</dt>
<dd>An estimate of the real spherical harmonic coefficients (geodesy normalized) of relief along an interface with density contrast <code>rho</code> that satisfies the Bouguer anomaly <code>ba</code>. The degree zero term corresponds to the mean radius of the relief.
</dd>
<dt><code>ba</code> : input, real*8, dimension (2, <code>lmaxcalc</code>+1, <code>lmaxcalc</code>+1)</dt>
<dd>The real spherical harmonic coefficients of the Bouguer anomaly referenced to a spherical interface <code>r0</code>.
</dd>
<dt><code>grid</code> : input, real*8, dimension (<code>lmax</code>+1, 2*<code>lmax</code>+1) for <code>gridtype</code> 1, (2*<code>lmax</code>+2, 2*<code>lmax</code>+2) for <code>gridtype</code> 2, (2*<code>lmax</code>+2, 4*<code>lmax</code>+4) for <code>gridtype</code> 3</dt>
<dd>The initial estimate for the radii of the interface evaluated on a grid corresponding to a function of maximum spherical harmonic degree <code>lmax</code>. This is calculated by a call to either <code>MakeGridGLQ</code> or <code>MakeGridDH</code>. This grid must contain the degree-0 average radius of the interface.
</dd>
<dt><code>lmax</code> : input, integer</dt>
<dd>The spherical harmonic bandwidth of the input relief <code>grid</code>, which determines the dimensions of <code>grid</code>. If <code>lmaxcalc</code> is not set, this determines also the maximum spherical harmonic degree of the output spherical harmonic coefficients of the relief and the input spherical harmonics of the Bouguer anomaly.
</dd>
<dt><code>nmax</code> : input, integer</dt>
<dd>The maximum order used in the Taylor-series expansion used in calculating the potential coefficients.
</dd>
<dt><code>mass</code> : input, real*8</dt>
<dd>The mass of the planet in kg.
</dd>
<dt><code>r0</code> : input, real*8</dt>
<dd>The reference radius of the Bouguer anomaly <code>ba</code>.
</dd>
<dt><code>rho</code> : input, real*8</dt>
<dd>The density contrast of the relief in kg/m^3.
</dd>
<dt><code>gridtype</code> : input, integer</dt>
<dd>1 = Gauss-Legendre grids, calculated using <code>SHGLQ</code> and <code>MakeGridGLQ</code>. 2 = Equally sampled Driscoll-Healy grids, <code>n</code> by <code>n</code>, calculated using <code>MakeGridDH</code>. 3 = Equally spaced Driscoll-Healy grids, <code>n</code> by 2<code>n</code>, calculated using <code>MakeGridDH</code>.
</dd>
<dt><code>w</code> : optional, input, real*8, dimension (<code>lmax</code>+1)</dt>
<dd>The weights used in the Gauss-Legendre quadrature. These are calculated from a call to <code>SHGLQ</code>. If present, one of <code>plx</code> or <code>zero</code> must also be present.
</dd>
<dt><code>plx</code> : optional, input, real*8, dimension (<code>lmax</code>+1, (<code>lmax</code>+1)*(<code>lmax</code>+2)/2)</dt>
<dd>An array of the associated Legendre functions calculated at the nodes used in the Gauss-Legendre quadrature. These are determined from a call to <code>SHGLQ</code>.
</dd>
<dt><code>zero</code> : optional, input, real*8, dimension (<code>lmax</code>+1)</dt>
<dd>The nodes used in the Gauss-Legendre quadrature over latitude, calculated by a call to <code>SHGLQ</code>.
</dd>
<dt><code>filtertype</code> : optional, input, integer, default = 0</dt>
<dd>Apply a filter when calculating the relief in order to minimize the destabilizing effects of downward continuation which amplify uncertainties in the Bouguer anomaly. If 0, no filtering is applied. If 1, use the minimum amplitude filter <code>DownContFilterMA</code>. If 2, use the minimum curvature filter <code>DownContFilterMC</code>.
</dd>
<dt><code>filterdeg</code> : optional, input, integer</dt>
<dd>The spherical harmonic degree for which the filter is 0.5.
</dd>
<dt><code>lmaxcalc</code> : optional, input, integer, default = <code>lmax</code></dt>
<dd>The maximum degree that will be calculated in the spherical harmonic expansions.
</dd>
<dt><code>exitstatus</code> : output, optional, integer</dt>
<dd>If present, instead of executing a STOP when an error is encountered, the variable exitstatus will be returned describing the error. 0 = No errors; 1 = Improper dimensions of input array; 2 = Improper bounds for input variable; 3 = Error allocating memory; 4 = File IO error.
</dd>
</dl>
<h1 id="description">Description</h1>
<p>I<code>BAtoHilm</code> is used to solve iteratively for the relief along an interface that corresponds to a given Bouguer anomaly. This is equation 18 of Wieczorek and Phillips (1998) which implicitly takes into consideration the finite-amplitude correction. Each iteration takes as input a guess for the relief (specified by <code>grid</code>) and outputs the iteratively improved spherical harmonic coefficients of this relief. These coefficients can then be re-expanded and re-input into this routine as the next guess. For the initial guess, it is often sufficient to use the relief predicted using the first-order “mass sheet” approximation. The input relief <code>grid</code> can be of one of three type specified by <code>gridtype</code>: 1 for Gauss-Legendre grids, 2 for equally sampled Driscoll-Healy grids (<code>n</code> by <code>n</code>), and 3 for equally spaced Driscoll-Healy grid (<code>n</code> by 2<code>n</code>).</p>
<p>If the algorithm does not converge, one might want to try damping the initial estimate. Alternatively, iterations of the following form have proven successfulin in damping oscilations between successive iterations:</p>
<p><code>h3 = (h2+h1)/2</code> <code>h4 = f(h3)</code></p>
<p>It is important to understand that as an intermediate step, this routine calculates the spherical harmonic coefficients of the relief raised to the nth power, i.e., grid**n. As such, if the input function is bandlimited to degree <code>L</code>, the resulting function will thus be bandlimited to degree <code>L*nmax</code>. This subroutine implicitly assumes that <code>lmax</code> is greater than or equal to <code>L*nmax</code>. If this is not the case, aliasing will occur. In practice, for accurate results, it is found that <code>lmax</code> needs only to be about twice the size of <code>L</code>, though this should be verified for each application. Thus, if the input function is considered to be bandlimited to degree <code>L</code>, the function should be evaluated on a grid corresponding to a maximum degree of about <code>2L</code>.</p>
<p>If the input grid is evaluated on the Gauss-Legendre points, it is necessary to specify the optional parameters <code>w</code> and<code>zero</code>, or <code>w</code> and <code>plx</code>, which are calculated by a call to <code>SHGLQ</code>. If memory is not an issue, the algorithm can be speeded up by inputing the optional array <code>plx</code> of precomputed associated Legendre functions on the Gauss-Legendre nodes. If <code>plx</code> is not specified, then it is necessary to input the optional array <code>zero</code> that contains the latitudinal Gauss-Legendre quadrature nodes.</p>
<p>This routine uses geodesy 4-pi normalized spherical harmonics that exclude the Condon-Shortley phase.</p>
<h1 id="references">References</h1>
<p>Wieczorek, M. A. and R. J. Phillips, Potential anomalies on a sphere: applications to the thickness of the lunar crust, J. Geophys. Res., 103, 1715-1724, 1998.</p>
<h1 id="see-also">See also</h1>
<p><a href="batohilmrhoh.html">batohilmrhoh</a>, <a href="cilmplus.html">cilmplus</a>, <a href="cilmplusrhoh.html">cilmplusrhoh</a>, <a href="shexpandglq.html">shexpandglq</a>, <a href="makegridglq.html">makegridglq</a>, <a href="shglq.html">shglq</a>, <a href="shexpanddh.html">shexpanddh</a>, <a href="makegriddh.html">makegriddh</a>, <a href="glqgridcoord.html">glqgridcoord</a>, <a href="downcontfilterma.html">downcontfilterma</a>, <a href="downcontfiltermc.html">downcontfiltermc</a></p>
<p class="dir">
> <a href="../../../index.html" class="dir">Home</a> > <a href="../../documentation.html" class="dir">Documentation</a> > <a href="../../f95-routines.html" class="dir">Fortran 95</a> > <a href="../../gravmag.html" class="dir">Gravity and Magnetics</a></p>
<table class="footer2">
<tbody>
<tr>
<td class="c1"><a href="https://lagrange.oca.eu/?lang=en">Laboratoire Lagrange</a></td>
<td class="c2"><a href="https://www.oca.eu/?lang=en">Observatoire de la Côte d'Azur</a></td>
<td class="c3">© 2017 <a href="https://github.com/SHTOOLS">SHTOOLS</a></td>
</tr>
</tbody>
</table>
</div>
</body>
</html>