mapper.py.html
Source file:   mapper.py
Directory:   /var/www/html/clawpack/links/awr10/radial-ocean-island
Converted:   Thu Jul 29 2010 at 13:49:17   using clawcode2html
This documentation file will not reflect any later changes in the source file.

 

from pylab import *
from pyclaw.data import Data


def latlong(d,theta,phi,Rearth):
    """
    Take a point at distance d from the North pole with longitude theta
    and return longitude and latitude of point after rotating pole down 
    to latitude phi.
    Applying this function to a set of points on a circle of radius d 
    will result in points that are equi-distant (distance d on the earth) 
    from the point at latitude phi, longitude 0.  

    Used to construct a radially symmetric ocean on the earth and place 
    gauges, islands, etc.
    """

    # Convert phi to radians:
    theta = theta * pi/180.
    phi = phi * pi/180.

    alpha = pi/2. - d/Rearth   # latitude of original point
    # x,y,z coordinates of this point on the earth:
    x1 = cos(alpha)*cos(theta)
    y1 = cos(alpha)*sin(theta)
    z1 = sin(alpha)

    # rotate so centered at latitude phi:
    x2 = x1*sin(phi) + z1*cos(phi)
    y2 = y1
    z2 = -x1*cos(phi) + z1*sin(phi)

    # compute longitude xhat and latitude yhat:
    xhat = -arctan(y2/x2)
    yhat = arcsin(z2)

    # convert to degrees:
    xhat = xhat * 180./pi
    yhat = yhat * 180./pi

    return xhat,yhat

def plot_ocean_and_shelf(d1=1580e3, d2=1645e3, phi=40.):

    theta = linspace(0, 360., 200)
    d = d1*ones(theta.shape)
    xhat, yhat = latlong(d, theta, phi)

    clf()
    plot(xhat,yhat,'b')
    hold(True)

    d = d2*ones(theta.shape)
    xhat, yhat = latlong(d, theta, phi)

    plot(xhat,yhat,'r')
    legend(['Continental shelf', 'Shoreline'],loc='lower left')

    (xi1,yi1) = latlong(1600.e3,220.,40.)
    plot([xi1],[yi1],'ko')
    (xi2,yi2) = latlong(1600.e3,260.,40.)
    plot([xi2],[yi2],'ko')
    axis('scaled')
    xlim([-20, 20])
    ylim([15, 60])