#!/usr/bin/python -tt
#=======================================================================
# General Documentation
"""Single-function module.
See function docstring for description.
"""
#-----------------------------------------------------------------------
# Additional Documentation
#
# RCS Revision Code:
# $Id: curl_2d.py,v 1.11 2004/05/06 23:33:37 jlin Exp $
#
# Modification History:
# - 01 Mar 2004: Original by Johnny Lin, Computation Institute,
# University of Chicago. Passed passably reasonable tests.
# - 04 May 2004: Add non-scalar R_sphere capability. Passed
# passably reasonable tests.
#
# Notes:
# - Written for Python 2.2.
# - Module docstrings can be tested using the doctest module. To
# test, execute "python curl_2d.py". Additional tests can be
# found in the test sub-directory of the distribution (see the
# package home page below for information).
# - Module has a single overall function, within which there are
# several private level-1 nested functions that define different
# methods of computing curl.
# - See import statements throughout module for packages and modules
# required.
#
# Copyright (c) 2004 by Johnny Lin. For licensing, distribution
# conditions, contact information, and additional documentation see
# the URL http://www.johnny-lin.com/py_pkgs/gemath/doc/.
#=======================================================================
#----------------------- Overall Module Imports ------------------------
#- The "from __future__ import" line as the first line in the module
# enables nested-scope access for Python versions prior to 2.2.
# From 2.2 and on nested-scope access is automatic:
from __future__ import nested_scopes
#- Set module version to package version:
import gemath_version
__version__ = gemath_version.version
__author__ = gemath_version.author
__date__ = gemath_version.date
__credits__ = gemath_version.credits
del gemath_version
#------------- Overall Function: Documentation and Import -------------
def curl_2d( x, y, Fx, Fy, missing=1e+20 \
, algorithm='default', R_sphere=6.37122e+6):
"""Curl of a vector F on a 2-D "rectangular" grid.
The 2-D grid F is defined on is rectangular, meaning that while
the grid spacings do not have to be even, each grid box is
rectangular and grid boundaries are rectilinear and parallel to
each other. If the 2-D grid is on a sphere, we assume that lines
of constant longitude and latitude form grids that are rectangular
(though in reality this is not strictly true).
Since F does not have a z-component, the curl of F is positive
pointing as a normal out of the plane defined by the basis for
x and y.
Positional Input Arguments:
* x: x-coordinate of each F. 1-D Numeric array with number of
elements equal to the number of columns in array F. Typically
x is in km, or on a sphere the longitude in deg. Floating or
integer type.
* y: y-coordinate of each F. 1-D Numeric array with number of
elements equal to the number of rows in array F. Typically
y is in km, or on a sphere the latitude in deg. Floating or
integer type.
* Fx: x-component of vector quantity F. 2-D Numeric array of
shape (len(y), len(x)). Floating or integer type.
* Fy: y-component of vector quantity F. 2-D Numeric array of
same shape and size as Fx. Floating or integer type.
Keyword Input Arguments:
* missing: If Fx and/or Fy has missing values, this is the
missing value value. Scalar. Default is 1e+20.
* algorithm: Name of the algorithm to use. String scalar.
Default is 'default'. Possible values include:
+ 'default': Default method. Algorithm 'order1_cartesian'
is used.
+ 'default_spherical': If 'spherepack' can be used, chooses
that method. Otherwise 'order1_spherical' is used. Either
way this option assumes the domain is on a sphere of radius
equal to keyword R_sphere.
+ 'order1_cartesian': First-order finite-differencing (back-
ward and forward differencing used at the endpoints, and
centered differencing used everywhere else). If abscissa
intervals are irregular, differencing will be correspon-
dingly asymmetric. Grid is assumed Cartesian. The curl
returned is in units of F divided by the units of x and y
(x and y must be in the same units).
+ 'order1_spherical': First-order finite-differencing (back-
ward and forward differencing used at the endpoints, and
centered differencing used everywhere else). If abscissa
intervals are irregular, differencing will be correspon-
dingly asymmetric. Grid is assumed spherical and x and y
are longitude and latitude (in deg). The curl returned is
in units of F divided by units of R_sphere. Note because
of singularities in the algorithm near the poles, the
curl for all points north of 88N or south of 88S latitude
are set to missing.
+ 'spherepack': The NCAR package SPHEREPACK 3.0 is used for
calculating the curl. The spatial domain is the global
sphere; x and y are longitude and latitude (in deg),
respectively; x must be evenly spaced and y must be gauss-
ian or evenly spaced. The domain must cover the entire
sphere and x and y must be longitude and latitude, respec-
tively, in deg. There can be no missing values. The
algorithm also assumes that x and y are monotonically in-
creasing from index 0. The curl returned is in units of F
divided by m. Thus, if F is velocity in m/s, the curl of
F is in units 1/s (and is the vorticity). Note that in
this function's implementation, the final curl output using
algorithm 'spherepack' is adjusted for R_sphere not equal
to the default mean earth radius, if appropriate, so this
algorithm can be used on other planets.
Note that SPHEREPACK operates on single precision floating
point (Float32) values. This function silently converts
input to Float32 for the computations, as needed (the input
parameters are not altered, however).
* R_sphere: If the grid is the surface of a sphere, this is
the radius of the sphere. This keyword is only used when
the algorithm is 'order1_spherical' or 'spherepack' (either
chosen explicitly by the algorithm keyword or automatically
when 'default_spherical' is used). Default value is the
"mean" radius of the Earth, in meters. "Mean" is in quotes
because this value changes depending on what you consider the
mean; the default in this function is the same Earth radius
used in module sphere. Note that if the level is a distance
z above the "mean" radius of the earth, R_sphere should be
set to the "mean" radius plus z.
Value can be a scalar or a Numeric or MA array of the same
size and shape as Fx. If keyword is an array, the value of
R_sphere for each element corresponds to the same location
as in the Fx and Fy arrays; there should also be no missing
values in R_sphere (function doesn't check for this, however).
Output:
* Curl of F. Units depend on algorithm used to calculate curl
(see above discussion of keyword algorithm). Numeric array of
same shape as Fx and Fy inputs.
If there are missing values, those elements in the output that
used missing values in the calculation of the curl are set to
the value in keyword missing. If there are missing values in
the output due to math errors and keyword missing is set to
None, output will fill those missing values with the MA default
value of 1e+20.
References:
* Glickman, T. S. (Ed.) (2000): "Curl," Glossary of Meteorology,
Boston, MA: American Meteorological Society, ISBN 1-878220-
34-9.
* Haltiner, G. J, and R. T. Williams (1980): Numerical Prediction
and Dynamic Meteorology, New York: John Wiley & Sons, pp. 8-9.
* Holton, J. R. (1992): An Introduction to Dynamic Meteorology,
San Diego, CA: Academic Press, p. 482.
* Kauffman, B. G., and W. G. Large (2002): The CCSM Coupler,
Version 5.0, User's Guide, Source Code Reference, and Scientific
Description. Boulder, CO: NCAR. Value for earth radius is
given at URL:
http://www.ccsm.ucar.edu/models/ccsm2.0/cpl5/users_guide/node10.html
* Overview of the CDAT Interface to the NCAR SPHEREPACK 3.0.
SPHEREPACK is by J. C. Adams and P. N. Swarztrauber. More
information on SPHEREPACK can be found at:
http://www.scd.ucar.edu/css/software/spherepack/
while information on the CDAT implementation used here is at:
http://esg.llnl.gov/cdat/getting_started/cdat.htm
Example of a synthetic dataset of winds. The curl is the vorti-
city. Images of the wind field and the vorticity are online at
http://www.johnny-lin.com/py_pkgs/gemath/doc/test_curl_2d_add.html.
Velocity at y greater than 60 and less than -60 are set to 0, to
prevent spectral algorithms from giving errors at the poles. Tiny
values less than 1e-10 are also set to 0:
(a) Import statements and create data:
>>> import Numeric as N
>>> from curl_2d import curl_2d
>>> nx = 181
>>> ny = 91
>>> x = N.arange(nx) * 2.0 - 180.0
>>> y = N.arange(ny) * 2.0 - 90.0
>>> y_2d = N.reshape( N.repeat(y, nx), (ny, nx) )
>>> x_2d = N.reshape( N.repeat(N.reshape(x,(1,nx)), ny), (ny, nx) )
>>> u = N.sin(x_2d*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> v = N.sin((x_2d-30.)*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> N.putmask( u, N.where(N.absolute(u) < 1e-10, 1, 0), 0. )
>>> N.putmask( v, N.where(N.absolute(v) < 1e-10, 1, 0), 0. )
>>> N.putmask( u, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )
>>> N.putmask( v, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )
(b) Use order 1 Cartesian algorithm: If x and y are in meters,
and u and v are in m/s, the curl is in 1/s:
>>> curl = curl_2d(x, y, u, v, algorithm='order1_cartesian')
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-0.007251593', '-0.003628007', '0', '0.003628007']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['0', '0', '0', '0']
(c) Use spectral method from SPHEREPACK (need to first remove the
repeating dateline point from the domain). Here we take x and
y in deg, while u and v are in m/s. If so the curl is in 1/s.
These results as much less than for example (b) above since
the domain is much larger:
>>> x = x[0:-1]
>>> u = u[:,0:-1]
>>> v = v[:,0:-1]
>>> curl = curl_2d(x, y, u, v, algorithm='spherepack')
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-6.436402e-08', '-3.220152e-08', '1.33852e-13', '3.220165e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['-2.299736e-20', '-4.806512e-21', '-9.283145e-22', '-1.368445e-20']
(d) Use order 1 spherical algorithm: x and y are in deg and u and
v are in m/s. The curl is in 1/s. Note that these results are
nearly identical to those from SPHEREPACK, except near the poles
the values are set to missing:
>>> curl = curl_2d(x, y, u, v, algorithm='order1_spherical')
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-6.517317e-08', '-3.260645e-08', '0', '3.260645e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['1e+20', '1e+20', '1e+20', '1e+20']
(e) Use "default" algorithm, which for this domain chooses the
order 1 cartesian algorithm:
>>> curl = curl_2d(x, y, u, v, algorithm='default')
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-0.007251593', '-0.003628007', '0', '0.003628007']
(f) Use "default_spherical" algorithm with missing data (equal to
1e+20), which because of the missing data will choose the order
1 spherical algorithm. We also do the calculation on Mars, and
illustrate passing in the radius as a constant as well as an
array of the same size and shape as u and v:
>>> u[30:46,90:114] = 1e+20
>>> curl = curl_2d( x, y, u, v, algorithm='default_spherical' \
, R_sphere=3390e+3)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-1.224875e-07', '-6.128107e-08', '1e+20', '1e+20']
>>> R_mars = N.zeros(u.shape, typecode=N.Float) + 3390e+3
>>> curl = curl_2d( x, y, u, v, algorithm='default_spherical' \
, R_sphere=R_mars )
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-1.224875e-07', '-6.128107e-08', '1e+20', '1e+20']
"""
import MA
import Numeric as N
from has_close import has_close
from can_use_sphere import can_use_sphere
#------- Nested Function: Curl By Order 1 Cartesian Finite Diff. ------
def _order1_cartesian_curl():
"""
Calculate curl using Cartesian first-order differencing.
An MA array is returned.
Algorithm: First-order differencing (interior points use
centered differencing, and end points use forward or backward
differencing, as applicable) in Cartesian coordinates (see
Glickman [2000], p. 194).
"""
from deriv import deriv
dFy_dx_N = N.zeros( (len(y), len(x)), typecode=N.Float )
dFx_dy_N = N.zeros( (len(y), len(x)), typecode=N.Float )
for iy in range(len(y)):
dFy_dx_N[iy,:] = deriv( x, N.ravel(Fy[iy,:]) \
, missing=missing, algorithm='order1')
for ix in range(len(x)):
dFx_dy_N[:,ix] = deriv( y, N.ravel(Fx[:,ix]) \
, missing=missing, algorithm='order1')
dFy_dx = MA.masked_values(dFy_dx_N, missing, copy=0)
dFx_dy = MA.masked_values(dFx_dy_N, missing, copy=0)
return dFy_dx - dFx_dy
#------- Nested Function: Curl By Order 1 Spherical Finite Diff. ------
def _order1_spherical_curl():
"""
Calculate curl using spherical first-order differencing.
An MA array is returned.
Algorithm: First-order differencing (interior points use
centered differencing, and end points use forward or backward
differencing, as applicable) in spherical coordinates (see
Holton 1992).
Note because of singularities in the algorithm near the poles,
the return value for all points north of 88N or south of 88S
latitude are set to mask true.
Key to some variables:
* xrad is x in radians (x is assumed in deg)
* yrad is y in radians (y is assumed in deg)
* ddFy is d(Fy) / d(xrad) as a masked array. The version with
a "_N" suffix means the Numeric version.
* ddFx is d(Fx * cos(yrad)) / d(yrad) as a masked array. The
version with a "_N" suffix means the Numeric version.
"""
from deriv import deriv
xrad = x * N.pi/180.0
yrad = y * N.pi/180.0
#- Derivative preliminaries:
ddFy_N = N.zeros( (len(yrad), len(xrad)), typecode=N.Float )
ddFx_N = N.zeros( (len(yrad), len(xrad)), typecode=N.Float )
for iy in range(len(yrad)):
ddFy_N[iy,:] = deriv( xrad, N.ravel(Fy[iy,:]) \
, missing=missing, algorithm='order1')
for ix in range(len(xrad)):
tmp_MA = MA.masked_values(N.ravel(Fx[:,ix]), missing, copy=0)
tmp_N = MA.filled(tmp_MA * N.cos(yrad), missing)
ddFx_N[:,ix] = deriv( yrad, tmp_N \
, missing=missing, algorithm='order1')
ddFy = MA.masked_values(ddFy_N, missing, copy=0)
ddFx = MA.masked_values(ddFx_N, missing, copy=0)
#- Calculate the curl:
tmpr = ( ddFy - ddFx ) / \
( R_sphere * N.reshape( N.repeat(N.cos(yrad),len(xrad)) \
, (len(yrad),len(xrad)) ) )
#- Make points "near" poles be missing and return:
np_mask_lat = 88.0
sp_mask_lat = -88.0
yarr = N.reshape( N.repeat(y,len(x)), (len(y),len(x)) )
np_mask = MA.make_mask( N.where(yarr > np_mask_lat, 1, 0) )
sp_mask = MA.make_mask( N.where(yarr < sp_mask_lat, 1, 0) )
tmpr_mask = tmpr.mask()
tmpr_mask = MA.mask_or(tmpr_mask, np_mask)
tmpr_mask = MA.mask_or(tmpr_mask, sp_mask)
return MA.masked_array(tmpr, mask=tmpr_mask)
#--------- Nested Function: Curl By SPHEREPACK Spectral Method --------
def _spherepack_curl():
"""
Calculate curl using SPHEREPACK spectral methods.
If the value of R_sphere is not the same as the earth radius
used in package sphere, adjust curl to reflect R_sphere. A
Numeric array is returned.
SPHEREPACK outputs a number of diagnostic messages to
stdout (and perhaps to stderr). The one that's irritating
is the one saying data will be converted to Float32. Thus,
we silently ensure that calculations are used only on Float32
data, to surpress the warnings. I did this because redirect-
ing stdout/stderr turned out to be too difficult; just
altering the sys attributes didn't work.
"""
import sphere, MA
sph_obj = sphere.Sphere( x.astype(MA.Float32) \
, y.astype(MA.Float32) )
if N.allclose(N.array(R_sphere), sphere.radius):
curl = sph_obj.vrt( Fx.astype(MA.Float32) \
, Fy.astype(MA.Float32) )
else:
curl = sph_obj.vrt( Fx.astype(MA.Float32) \
, Fy.astype(MA.Float32) ) \
* sphere.radius / R_sphere
return curl
#-------------- Overall Function: Settings, Body, Return --------------
#- Choose algorithm to compute curl:
if algorithm == 'default':
_calculate_curl = _order1_cartesian_curl
elif algorithm == 'default_spherical':
if (can_use_sphere(x,y)[0] == 1) and \
(not has_close(Fx, missing)) and \
(not has_close(Fy, missing)):
_calculate_curl = _spherepack_curl
else:
_calculate_curl = _order1_spherical_curl
elif algorithm == 'order1_cartesian':
_calculate_curl = _order1_cartesian_curl
elif algorithm == 'order1_spherical':
_calculate_curl = _order1_spherical_curl
elif algorithm == 'spherepack':
if has_close(Fx, missing) or has_close(Fy, missing):
raise ValueError, "curl_2d: has missing values"
else:
_calculate_curl = _spherepack_curl
else:
raise ValueError, "curl_2d: bad algorithm"
#- Calculate curl and return from function:
return MA.filled( _calculate_curl(), missing )
#-------------------------- Main: Test Module -------------------------
#- Define additional examples for doctest to use:
__test__ = { 'Additional Examples':
"""
>>> import Numeric as N
>>> from curl_2d import curl_2d
>>> nx = 181
>>> ny = 91
>>> x = N.arange(nx) * 2.0 - 180.0
>>> y = N.arange(ny) * 2.0 - 90.0
>>> y_2d = N.reshape( N.repeat(y, nx), (ny, nx) )
>>> x_2d = N.reshape( N.repeat(N.reshape(x,(1,nx)), ny), (ny, nx) )
>>> u = N.sin(x_2d*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> v = N.sin((x_2d-30.)*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> x = x[0:-1]
>>> u = u[:,0:-1]
>>> v = v[:,0:-1]
(i) Use SPHEREPACK on Mars:
>>> curl = curl_2d(x, y, u, v, algorithm='spherepack', R_sphere=3390e+3)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-1.213008e-07', '-6.068711e-08', '2.165924e-13', '6.068734e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['-1.150181e-20', '1.700353e-20', '1.227383e-20', '-2.283483e-20']
>>> ['%.7g' % curl[-1,i] for i in range(10,14)]
['1.596643e-20', '5.265175e-21', '4.54957e-21', '-1.17195e-20']
>>> ['%.7g' % curl[70,i] for i in range(130,134)]
['1.474248e-06', '1.470224e-06', '1.464409e-06', '1.456811e-06']
(ii) Use order 1 spherical on Mars:
>>> curl = curl_2d(x,y,u,v,algorithm='order1_spherical',R_sphere=3390e+3)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-1.224875e-07', '-6.128107e-08', '-9.383406e-23', '6.128107e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['1e+20', '1e+20', '1e+20', '1e+20']
>>> ['%.7g' % curl[-1,i] for i in range(10,14)]
['1e+20', '1e+20', '1e+20', '1e+20']
>>> ['%.7g' % curl[70,i] for i in range(130,134)]
['1.413254e-06', '1.408895e-06', '1.402819e-06', '1.395035e-06']
(iii) Use SPHEREPACK with a bad spatial domain (no poles):
>>> y2 = y[1:-1]
>>> curl = curl_2d(x, y2, u, v, algorithm='spherepack')
Traceback (most recent call last):
...
ValueError
(iv) Use SPHEREPACK with a bad spatial domain (repeating longitude):
>>> import MA
>>> x2 = MA.concatenate([x, [x[0]+360.]])
>>> curl = curl_2d(x, y2, u, v, algorithm='spherepack')
Traceback (most recent call last):
...
ValueError
(v) Use SPHEREPACK with missing values:
>>> u2 = u
>>> u2[3:5,5:8] = 1e+20
>>> curl = curl_2d(x, y2, u, v, algorithm='spherepack')
Traceback (most recent call last):
...
ValueError: curl_2d: has missing values
""",
'Additional Examples for Non-Scalar R_sphere: Data Set 1':
"""
>>> import Numeric as N
>>> from curl_2d import curl_2d
>>> nx = 181
>>> ny = 91
>>> x = N.arange(nx) * 2.0 - 180.0
>>> y = N.arange(ny) * 2.0 - 90.0
>>> y_2d = N.reshape( N.repeat(y, nx), (ny, nx) )
>>> x_2d = N.reshape( N.repeat(N.reshape(x,(1,nx)), ny), (ny, nx) )
>>> u = N.sin(x_2d*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> v = N.sin((x_2d-30.)*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> x = x[0:-1]
>>> u = u[:,0:-1]
>>> v = v[:,0:-1]
(i) Use SPHEREPACK on Mars:
>>> R_mars = N.zeros(u.shape, typecode=N.Float) + 3390e+3
>>> R_mars[45,88] = 3390e+3 + 1000.0
>>> curl = curl_2d(x, y, u, v, algorithm='spherepack', R_sphere=R_mars)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-1.21265e-07', '-6.068711e-08', '2.165924e-13', '6.068734e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['-1.150181e-20', '1.700353e-20', '1.227383e-20', '-2.283483e-20']
>>> ['%.7g' % curl[-1,i] for i in range(10,14)]
['1.596643e-20', '5.265175e-21', '4.54957e-21', '-1.17195e-20']
>>> ['%.7g' % curl[70,i] for i in range(130,134)]
['1.474248e-06', '1.470224e-06', '1.464409e-06', '1.456811e-06']
(ii) Use order 1 spherical on Mars:
>>> R_mars = N.zeros(u.shape, typecode=N.Float) + 3390e+3
>>> R_mars[0,8] = 3390e+3 + 2000.0
>>> curl = curl_2d(x,y,u,v,algorithm='order1_spherical',R_sphere=R_mars)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-1.224875e-07', '-6.128107e-08', '-9.383406e-23', '6.128107e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['1e+20', '1e+20', '1e+20', '1e+20']
>>> ['%.7g' % curl[-1,i] for i in range(10,14)]
['1e+20', '1e+20', '1e+20', '1e+20']
>>> ['%.7g' % curl[70,i] for i in range(130,134)]
['1.413254e-06', '1.408895e-06', '1.402819e-06', '1.395035e-06']
""",
'Additional Examples for Non-Scalar R_sphere: Data Set 2':
"""
(a) Import statements and create data:
>>> import MA
>>> import Numeric as N
>>> from curl_2d import curl_2d
>>> nx = 181
>>> ny = 91
>>> x = N.arange(nx) * 2.0 - 180.0
>>> y = N.arange(ny) * 2.0 - 90.0
>>> y_2d = N.reshape( N.repeat(y, nx), (ny, nx) )
>>> x_2d = N.reshape( N.repeat(N.reshape(x,(1,nx)), ny), (ny, nx) )
>>> u = N.sin(x_2d*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> v = N.sin((x_2d-30.)*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> N.putmask( u, N.where(N.absolute(u) < 1e-10, 1, 0), 0. )
>>> N.putmask( v, N.where(N.absolute(v) < 1e-10, 1, 0), 0. )
>>> N.putmask( u, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )
>>> N.putmask( v, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )
(b) Use order 1 Cartesian algorithm: If x and y are in meters,
and u and v are in m/s, the curl is in 1/s. Note that the
R_sphere keyword value has no effect with this algorithm:
>>> R_earth = MA.zeros(u.shape, typecode=N.Float) + 6371220.0
>>> R_earth[45,89] = 6371220.0 + 2000.0
>>> curl = curl_2d( x, y, u, v, algorithm='order1_cartesian' \
, R_sphere=R_earth )
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-0.007251593', '-0.003628007', '0', '0.003628007']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['0', '0', '0', '0']
(c) Use spectral method from SPHEREPACK (need to first remove the
repeating dateline point from the domain). Here we take x and
y in deg, while u and v are in m/s. If so the curl is in 1/s.
These results as much less than for example (b) above since
the domain is much larger. We do two cases, one where the array
is not all the same as the default Earth radius, and one where
it is all the same as the default Earth radius:
>>> x = x[0:-1]
>>> u = u[:,0:-1]
>>> v = v[:,0:-1]
>>> R_earth = MA.zeros(u.shape, typecode=N.Float) + 6371220.0
>>> R_earth[45,89] = 6371220.0 + 2000.0
>>> curl = curl_2d(x, y, u, v, algorithm='spherepack', R_sphere=R_earth)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-6.436402e-08', '-3.219142e-08', '1.33852e-13', '3.220165e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['-2.299736e-20', '-4.806512e-21', '-9.283145e-22', '-1.368445e-20']
>>> R_earth = MA.zeros(u.shape, typecode=N.Float) + 6371220.0
>>> curl = curl_2d(x, y, u, v, algorithm='spherepack', R_sphere=R_earth)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-6.436402e-08', '-3.220152e-08', '1.33852e-13', '3.220165e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['-2.299736e-20', '-4.806512e-21', '-9.283145e-22', '-1.368445e-20']
(d) Use order 1 spherical algorithm: x and y are in deg and u and
v are in m/s. The curl is in 1/s. Note that these results are
nearly identical to those from SPHEREPACK:
>>> R_earth = MA.zeros(u.shape, typecode=N.Float) + 6371220.0
>>> R_earth[45,89] = 6371220.0 + 2000.0
>>> curl = curl_2d( x, y, u, v, algorithm='order1_spherical' \
, R_sphere=R_earth )
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-6.517317e-08', '-3.259621e-08', '0', '3.260645e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['1e+20', '1e+20', '1e+20', '1e+20']
(e) Use "default" algorithm, which for this domain chooses the
order 1 cartesian algorithm. Thus, the R_sphere keyword
setting should have no effect:
>>> R_earth = MA.zeros(u.shape, typecode=N.Float) + 6371220.0
>>> R_earth[45,89] = 6371220.0 + 2000.0
>>> curl = curl_2d(x, y, u, v, algorithm='default', R_sphere=R_earth)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-0.007251593', '-0.003628007', '0', '0.003628007']
(f) Use "default_spherical" algorithm with missing data (equal to
1e+20), which because of the missing data will choose the order
1 spherical algorithm. We also do the calculation on Mars:
>>> R_mars = N.zeros(u.shape, typecode=N.Float) + 3390e+3
>>> R_mars[45,90] = 3390e+3 + 1500.0
>>> R_mars[45,88] = 3390e+3 + 5000.0
>>> u[30:46,90:114] = 1e+20
>>> curl = curl_2d( x, y, u, v, algorithm='default_spherical' \
, R_sphere=R_mars)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-1.223071e-07', '-6.128107e-08', '1e+20', '1e+20']
"""
}
#- Execute doctest if module is run from command line:
if __name__ == "__main__":
"""Test the module.
Tests the examples in all the module documentation strings, plus
__test__.
Note: To help ensure that module testing of this file works, the
parent directory to the current directory is added to sys.path.
"""
import doctest, sys, os
sys.path.append(os.pardir)
doctest.testmod(sys.modules[__name__])
# ===== end file =====