How to calculate smallcircle on sphere

Question

I use this python code to calculate a greatcircle: https://code.google.com/p/pyroms/source/browse/trunk/pyroms/greatcircle.py?r=39 After initializing the greatcircle I can calculate points on the line.

I would like to do the same with an arbitrary smallcircle. I would like to define the smallcircle using its center location (in latlon) and one location on the line (also latlon), so not using a diameter. Say as an input I would use a certain point spacing in meters, how would I then proceed to calculate the other locations on the smallcircle line?

Thanks!

Answer 1

Python:

# From: http://www.movable-type.co.uk/scripts/latlong.html

def distance(lon1, lat1, lon2, lat2):
    # http://code.activestate.com/recipes/576779-calculating-distance-between-two-geographic-points/
    # http://en.wikipedia.org/wiki/Haversine_formula
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    q = sin(dlat/2)**2 + (cos(lat1) * cos(lat2) * (sin(dlon/2)**2))
    return 2 * atan2(sqrt(q), sqrt(1-q))

from math import *

def SmallCircle(radius,clon,clat,llon,llat,meters):
    # clon = circle central longitude (deg)
    # clat = circle central latitude (deg)
    # llon = longitude of point on circle line (deg)
    # llat = latitude of point on circle line (deg)
    # meters = point spacing (meters). This is not the exact spacing, it will be used to calculate nrpoints.
    # radius = radius of sphere

    dlat = radians(llat-clat)
    dlon = radians(llon-clon)
    clat = radians(clat)
    clon = radians(clon)
    llat = radians(llat)
    llon = radians(llon)

    a = sin(dlat/2) * sin(dlat/2) + sin(dlon/2) * sin(dlon/2) * cos(clat) * cos(llat)
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    d = radius * c

    # bearing
    y = sin(dlon) * cos(llat)
    x = cos(clat)*sin(llat) - sin(clat)*cos(llat)*cos(dlon)
    brng = atan2(y, x)

    # circle
    angle = distance(clon,clat,llon,llat)
    circleradius = sin(angle) * radius
    circlelength = 2*pi*circleradius

    # loop circle
    list = []
    list.append([degrees(llon),degrees(llat)])
    nrpoints = int(round(circlelength / meters))
    for i in range(1, nrpoints):
        brng = brng + (float(i) / float(nrpoints))*2*pi
        lat = asin( sin(clat)*cos(d/radius) + cos(clat)*sin(d/radius)*cos(brng) )
        lon = clon + atan2(sin(brng)*sin(d/radius)*cos(clat), cos(d/radius)-sin(clat)*sin(lat))
        list.append([degrees(lon),degrees(lat)])

    return list