Skip to content

Instantly share code, notes, and snippets.

@ajdawson
Created January 21, 2015 14:04
Show Gist options
  • Save ajdawson/b64d24dfac618b91974f to your computer and use it in GitHub Desktop.
Save ajdawson/b64d24dfac618b91974f to your computer and use it in GitHub Desktop.
Gaussian grid generation.
"""Tools for working with Gaussian grids."""
from __future__ import (absolute_import, division, print_function)
import functools
import numpy as np
import numpy.linalg as la
from numpy.polynomial.legendre import legcompanion, legder, legval
def __single_arg_fast_cache(func):
"""Caching decorator for functions of one argument."""
class CachingDict(dict):
def __missing__(self, key):
result = self[key] = func(key)
return result
@functools.wraps(func)
def __getitem__(self, *args, **kwargs):
return super(CachingDict, self).__getitem__(*args, **kwargs)
return CachingDict().__getitem__
@__single_arg_fast_cache
def gaussian_latitudes(n):
"""Construct latitudes and latitude bounds for a Gaussian grid.
Args:
* n:
The Gaussian grid number (half the number of latitudes in the
grid.
Returns:
A 2-tuple where the first element is a length `n` array of
latitudes (in degrees) and the second element is an `(n, 2)`
array of bounds.
"""
if abs(int(n)) != n:
raise ValueError('n must be a non-negative integer')
nlat = 2 * n
# Create the coefficients of the Legendre polynomial and construct the
# companion matrix:
cs = np.array([0] * nlat + [1], dtype=np.int)
cm = legcompanion(cs)
# Compute the eigenvalues of the companion matrix (the roots of the
# Legendre polynomial) taking advantage of the fact that the matrix is
# symmetric:
roots = la.eigvalsh(cm)
roots.sort()
# Improve the roots by one application of Newton's method, using the
# solved root as the initial guess:
fx = legval(roots, cs)
fpx = legval(roots, legder(cs))
roots -= fx / fpx
# The roots should exhibit symmetry, but with a sign change, so make sure
# this is the case:
roots = (roots - roots[::-1]) / 2.
# Compute the Gaussian weights for each interval:
fm = legval(roots, cs[1:])
fm /= np.abs(fm).max()
fpx /= np.abs(fpx).max()
weights = 1. / (fm * fpx)
# Weights should be symmetric and sum to two (unit weighting over the
# interval [-1, 1]):
weights = (weights + weights[::-1]) / 2.
weights *= 2. / weights.sum()
# Calculate the bounds from the weights, still on the interval [-1, 1]:
bounds1d = np.empty([nlat + 1])
bounds1d[0] = -1
bounds1d[1:-1] = -1 + weights[:-1].cumsum()
bounds1d[-1] = 1
# Convert the bounds to degrees of latitude on [-90, 90]:
bounds1d = np.rad2deg(np.arcsin(bounds1d))
bounds2d = np.empty([nlat, 2])
bounds2d[:, 0] = bounds1d[:-1]
bounds2d[:, 1] = bounds1d[1:]
# Convert the roots from the interval [-1, 1] to latitude values on the
# interval [-90, 90] degrees:
latitudes = np.rad2deg(np.arcsin(roots))
return latitudes, bounds2d
@guziy
Copy link

guziy commented Jun 8, 2017

Is it OK to put this file in my library?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment