Contents:
Planetary mapping routines.
phi = 0 faces toward the observer phi = pi thus faces away from the observer theta=pi/2 is the zaxis or ‘north pole’ theta=pi/2 is the ‘south pole’
Give the chisquared to an input timeseries to get the wedge coefficients and inclination angle.
param – [inclination, phi0, wedge coefficients]
phi,theta – from MAKEGRID
rangle – rotation angles at which flux was measured
meas – measured flux
err – error on measurement (onesigma)
retmodel – whether to return the model, instead of the chisquared.
Make grids of phi and theta values with the specified number of points in each direction. Phi ranges from 0 to 2pi, and theta ranges from pi/2 to pi/2.
Returns (phi, theta)
INPUTS: 


EXAMPLE:  import maps
nlat, nlon = 60, 30
phi, theta = maps.makegrid(nlat, nlon)
# Make a small spot centered near, but not at, the equator:
equator_spot = maps.makespot(0, 0, 0.4, phi, theta)
# Make a larger spot centered near, but not at, the pole:
pole_spot = maps.makespot(1.2, 0, 0.7, phi, theta)
import maps
nlat, nlon = 60, 30
map = maps.map(nlat, nlon, i=0., deltaphi=0.)
phi = map.corners_latlon.mean(2)[:,1].reshape(nlon, nlat)
theta = map.corners_latlon.mean(2)[:,0].reshape(nlon, nlat)  np.pi/2.
# Make a small spot centered near, but not at, the equator:
equator_spot = maps.makespot(0, 0, 0.4, phi, theta)
# Make a larger spot centered near, but not at, the pole:
pole_spot = maps.makespot(1.2, 0, 0.7, phi, theta)

Make a spotmap.
EXAMPLE: import maps inc, rot = pi/2., 0. phi, theta = maps.makegrid(120,60) phi2, theta2 = maps.rotcoord(phi,theta,inc,rot) vmap = maps.visiblemap(phi,theta,inc,rot)
maps.makespot(phi,theta,inc=pi/2,rot=0.,lat=pi/4,long=1.*pi) maps.makespot(phi2,theta2,vmap=vmap,lat=pi/4,long=1.*pi)
Very handy spherical mapping object. :INPUTS:
 nlon, nlat : scalars
 If mod==’latlon’, these inputs specify the number of grid cells across map, in latitude and longitude.
 i : scalar
 the inclination, is in units of radians. Zero means we see the object equatoron; pi/2 means we see it poleon.
 type : str
 For now, only valid entry is ‘latlon’. Eventually, this could be expanded to allow other types of projection grids.
 deltaphi : scalar
 Rotation of map, specified in radians.
OUTPUT:  A mapclass object with various useful fields. Most of these fields refer to the coordinates (either Cartesian or spherical polar) or the projected radial velocities at the corners of specified grid cells, or the approximate projected areas of these grid cells. 

NOTES:  I have not been as careful as I should be in this code – my original goal was speed rather than exactitude. This means that some values are returned as ‘nan’, and the projected areas are only roughly correct. There’s plenty of room for improvement! 
Compute the area of a polygon whose vertices are at the points (x,y).
INPUTS: 


REFERENCE: 
Return a 2D map of the projected area. Note that you need to determine for yourself which parts of the areal map are actually visible to the observer.
Assumes phi and theta are grids straight out of meshgrid.
rotate coordinate system from local (planetary) coordinates into observeroriented coordinages. rangle and iangle are the rotation and inclination angles of the planet in radians.
returns (phi2,theta2)
phi2 will be in the range (0,2pi) and theta2 will be in the range (pi/2,pi/2)
Return a 2D boolean map that’s True for the planetary latitude, longitude values visible from an observer. For rangle (rotation angle) zero and iangle (inclination angle) equal to pi/2, phi=pi is toward the observer.
Return a set of basis functions for the flux from each wedge.
phi,theta are the observercentered grids from MAPS.MAKEGRID iangle is the inclination of the system (0 = poleon) rangle (seq.) is a sequence of (0,2pi) rotation values at which
the wedgebased flux is evaluated.
nwedge is the number of wedges phi0 defines the center of wedge zero.
fwedge – an optional sequence to set individual flux values for each wedge.