geoclaw.util.haversine behaves differently than 4.x routine gcdist

[rjleveque]

Both compute the great circle distance, but the old routine took x1,y1,x1,y2 and computed distance from (x1,y1) to (x2,y2) while the new one takes x,y and computes distance from (x[0],x[1]) to (y[0],y[1]) so y has a different meaning than before. Since (x,y) are how we denote the 2d coordinates in general I prefer the previous interpretation.

Also the new version requires x and y to be of length 2 and would have to be called repeatedly to fill in an array, as in the below example where it is desired to compute the distance of each point in arrays x,y to a scalar center (x0,y0).

(This is from the radial ocean example I am trying to convert to 5.x form from http://depts.washington.edu/clawpack/links/awr11/)

def qinit(x,y):
    """
    Gaussian hump:
    """
    from numpy import where
    x0 = 0.0
    y0 = 40.0
    d = gcdist(x0,y0,x,y,Rearth)
    ze = -d**2/2.e9
    z = where(ze>-100., 20.e0*exp(ze), 0.)
    return z

[mandli]

I would favor making it so this was vectorizable and consistently named (probably my fault on that). Besides that do the routines actually return different values though?

[rjleveque]

The formulas looked ok but I haven't checked carefully or tested the output.