~ubuntu-branches/ubuntu/wily/gsw/wily

« back to all changes in this revision

Viewing changes to gsw/gibbs/earth.py

  • Committer: Package Import Robot
  • Author(s): Alastair McKinstry
  • Date: 2012-12-17 19:12:21 UTC
  • Revision ID: package-import@ubuntu.com-20121217191221-k7d89tsyikdycbra
Tags: upstream-3.0.1
ImportĀ upstreamĀ versionĀ 3.0.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
# -*- coding: utf-8 -*-
 
2
 
 
3
from __future__ import division
 
4
 
 
5
import numpy as np
 
6
 
 
7
from constants import gamma, earth_radius, OMEGA
 
8
from gsw.utilities import match_args_return
 
9
from conversions import z_from_p
 
10
 
 
11
 
 
12
__all__ = ['f',
 
13
           'grav',
 
14
           'distance']
 
15
 
 
16
RAD = np.pi / 180.0
 
17
 
 
18
 
 
19
def f(lat):
 
20
    r"""Calculates the Coriolis parameter (f) defined by:
 
21
        f = 2*omega*sin(lat)
 
22
    where,
 
23
        omega = 7.292115e-5 (Groten, 2004) [radians s :sup:`-1`]
 
24
 
 
25
    Parameters
 
26
    ----------
 
27
    lat : array_like
 
28
          latitude [degrees north]
 
29
 
 
30
    Returns
 
31
    -------
 
32
    f : array_like
 
33
        Coriolis paramter  [s :sup:`-1`]
 
34
 
 
35
    References
 
36
    ----------
 
37
    .. [1] Groten, E., 2004: Fundamental Parameters and Current (2004) Best
 
38
    Estimates of the Parameters of Common Relevance to Astronomy, Geodesy, and
 
39
    Geodynamics. Journal of Geodesy, 77, pp. 724-797.
 
40
 
 
41
    .. [2] IOC, SCOR and IAPSO, 2010: The international thermodynamic equation
 
42
    of seawater -  2010: Calculation and use of thermodynamic properties.
 
43
    Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
 
44
    UNESCO (English), 196 pp.
 
45
 
 
46
    Modifications:
 
47
    1993-04-20. Phil Morgan
 
48
    2010-07-28. Paul Barker
 
49
    """
 
50
 
 
51
    lat = np.asanyarray(lat)
 
52
    return 2 * OMEGA * np.sin(lat * RAD)
 
53
 
 
54
 
 
55
@match_args_return
 
56
def grav(lat, p=0):
 
57
    r"""Calculates acceleration due to gravity as a function of latitude and as
 
58
    a function of pressure in the ocean.
 
59
 
 
60
    Parameters
 
61
    ----------
 
62
    lat : array_like
 
63
          latitude in decimal degrees north [-90...+90]
 
64
    p : number or array_like. Default p = 0
 
65
        pressure [dbar]
 
66
 
 
67
    Returns
 
68
    -------
 
69
    g : array_like
 
70
        gravity [m s :sup:`2`]
 
71
 
 
72
    See Also
 
73
    --------
 
74
    TODO
 
75
 
 
76
    Notes
 
77
    -----
 
78
    In the ocean z is negative.
 
79
 
 
80
    Examples
 
81
    --------
 
82
    >>> import gsw
 
83
    >>> lat = [-90, -60, -30, 0]
 
84
    >>> p = 0
 
85
    >>> gsw.grav(lat, p)
 
86
    array([ 9.83218621,  9.81917886,  9.79324926,  9.780327  ])
 
87
    >>> gsw.grav(45)
 
88
    9.8061998770458008
 
89
 
 
90
    References
 
91
    ----------
 
92
    .. [1] IOC, SCOR and IAPSO, 2010: The international thermodynamic equation
 
93
    of seawater -  2010: Calculation and use of thermodynamic properties.
 
94
    Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
 
95
    UNESCO (English), 196 pp.
 
96
 
 
97
    .. [2] Moritz (2000) Goedetic reference system 1980. J. Geodesy, 74,
 
98
    128-133.
 
99
 
 
100
    .. [3] Saunders, P.M., and N.P. Fofonoff (1976) Conversion of pressure to
 
101
    depth in the ocean. Deep-Sea Res.,pp. 109 - 111.
 
102
 
 
103
    Modifications:
 
104
    2011-03-29. Trevor McDougall & Paul Barker
 
105
    """
 
106
 
 
107
    X = np.sin(lat * RAD)
 
108
    sin2 = X ** 2
 
109
    gs = 9.780327 * (1.0 + (5.2792e-3 + (2.32e-5 * sin2)) * sin2)
 
110
    z = z_from_p(p, lat)
 
111
    # z is the height corresponding to p.
 
112
    grav = gs * (1 - gamma * z)
 
113
 
 
114
    return grav
 
115
 
 
116
 
 
117
@match_args_return
 
118
def distance(lon, lat, p=0):
 
119
    r"""Calculates the distance in met res between successive points in the
 
120
    vectors lon and lat, computed using the Haversine formula on a spherical
 
121
    earth of radius 6,371 km, being the radius of a sphere having the same
 
122
    volume as Earth.  For a spherical Earth of radius 6,371,000 m, one nautical
 
123
    mile is 1,853.2488 m, thus one degree of latitude is 111,194.93 m.
 
124
 
 
125
    Haversine formula:
 
126
        R = earth's radius (mean radius = 6,371 km)
 
127
 
 
128
    .. math::
 
129
        a = \sin^2(\delta \text{lat}/2) +
 
130
            \cos(\text{lat}_1) \cos(\text{lat}_2) \sin^2(\delta \text{lon}/2)
 
131
 
 
132
        c = 2 \times \text{atan2}(\sqrt{a}, \sqrt{(1-a)})
 
133
 
 
134
        d = R \times c
 
135
 
 
136
    Parameters
 
137
    ----------
 
138
    lon : array_like
 
139
          decimal degrees east [0..+360] or [-180 ... +180]
 
140
    lat : array_like
 
141
          latitude in decimal degrees north [-90..+90]
 
142
    p : number or array_like. Default p = 0
 
143
        pressure [dbar]
 
144
 
 
145
    Returns
 
146
    -------
 
147
    dist: array_like
 
148
          distance between points on a spherical Earth at pressure (p) [m]
 
149
 
 
150
    See Also
 
151
    --------
 
152
    TODO
 
153
 
 
154
    Notes
 
155
    -----
 
156
    z is height and is negative in the oceanographic.
 
157
 
 
158
    Distances are probably good to better than 1\% of the "true" distance on
 
159
    the ellipsoidal earth.
 
160
 
 
161
    Examples
 
162
    --------
 
163
    >>> import gsw
 
164
    >>> lon = [159, 220]
 
165
    >>> lat = [-35, 35]
 
166
    >>> gsw.distance(lon, lat)
 
167
    array([[ 10030974.652916]])
 
168
    >>> p = [200, 1000]
 
169
    >>> gsw.distance(lon, lat, p)
 
170
    array([[ 10030661.63878009]])
 
171
    >>> p = [[200], [1000]]
 
172
    >>> gsw.distance(lon, lat, p)
 
173
    array([[ 10030661.63878009],
 
174
           [ 10029412.58776001]])
 
175
 
 
176
    References
 
177
    ----------
 
178
    .. [1] http://www.eos.ubc.ca/~rich/map.html
 
179
 
 
180
    Modifications:
 
181
    2000-11-06. Rich Pawlowicz
 
182
    2011-04-04. Paul Barker and Trevor McDougall
 
183
    """
 
184
    #FIXME? The argument handling seems much too complicated.
 
185
    # Maybe we can come up with some simple specifications of
 
186
    # what argument combinations are permitted, and handle everything
 
187
    # with broadcasting. - EF
 
188
 
 
189
    #FIXME: Eric what do you think? This assume p(stations, depth)
 
190
    lon, lat, = np.atleast_2d(lon), np.atleast_2d(lat)
 
191
 
 
192
    if (lon.size == 1) & (lat.size == 1):
 
193
        raise ValueError('more than one point is needed to compute distance')
 
194
    elif lon.ndim != lat.ndim:
 
195
        raise ValueError('lon, lat must have the same dimension')
 
196
 
 
197
    lon, lat, p = np.broadcast_arrays(lon, lat, p)
 
198
 
 
199
    dlon = np.diff(lon * RAD)
 
200
    dlat = np.diff(lat * RAD)
 
201
 
 
202
    a = ((np.sin(dlat / 2.)) ** 2 + np.cos(lat[:, :-1] * RAD) *
 
203
           np.cos(lat[:, 1:] * RAD) * (np.sin(dlon / 2.)) ** 2)
 
204
 
 
205
    # FIXME: Do we need np.ma here?
 
206
    angles = 2. * np.arctan2(np.ma.sqrt(a), np.ma.sqrt(1 - a))
 
207
 
 
208
    p_mid = 0.5 * (p[:, 0:-1] + p[:, 0:-1])
 
209
    lat_mid = 0.5 * (lat[:, :-1] + lat[:, 1:])
 
210
 
 
211
    z = z_from_p(p_mid, lat_mid)
 
212
 
 
213
    distance = (earth_radius + z) * angles
 
214
 
 
215
    return distance
 
216
 
 
217
 
 
218
if __name__ == '__main__':
 
219
    import doctest
 
220
    doctest.testmod()