distance between 2 points

[sparkz_alot]

This is a small part of a bigger project I'm working on to Pythonize (a real word I just made up) the formulas in Jean Meeus' book "Astronomical Algorithms" copywrite 1991, pages 80 and 81 and calculates the geodesic distance between two points on the earths surface. The two co-ordinates referenced at the end are, per Mr. Meeus for the Observatoire de Paris (France) and the US Naval Observatory (Washington DC) with an accuracy of 50 meters. Comments and improvements are welcome.

#! /usr/bin/python3
# -*- coding: utf-8 -*-

from math import atan, cos, pi, sin, sqrt


def calc_geo_dist(lon1, lat1, lon2, lat2):
    a = 6378.137    # in km, semi-major axis
    # b = 6356.752    # polar radius in km, rounded  approx 6356752.314140347 meters per GRS 80
    rad = pi/180
    flat = 1/298.257    # ok with nasa    actual 0.003352810681183637418
    e = ((lon1 - lon2) / 2) * rad
    f = ((lat1+lat2) / 2) * rad
    g = ((lat1-lat2) / 2) * rad
    s = sin(g) * sin(g) * cos(e) * cos(e) + cos(f) * cos(f) * sin(e) * sin(e)
    c = cos(g) * cos(g) * cos(e) * cos(e) + sin(f) * sin(f) * sin(e) * sin(e)
    w = atan(sqrt(s/c))
    r = sqrt(s * c) / w
    d = 2 * w * a
    h1 = (3 * r - 1) / (2 * c)
    h2 = (3 * r + 1) / (2 * s)
    s1 = d * (1 + flat * h1 * sin(f) * sin(f) * cos(g) * cos(g) - flat * h2 * cos(f) * cos(f) * sin(g) * sin(g))
    return s1

while True:
    try:
        my_long1 = float(input('\nEnter longitude in decimal degrees: '))
        if -180 <= my_long1 <= 180:
            break
        else:
            print('\n WARNING: Value must be between -180 and 180 ')
    except ValueError as err:
        print('\n ERROR: Enter value between -180 and 180 {}'.format(err))

while True:
    try:
        my_lat1 = float(input('Enter latitude in decimal degrees: '))
        if -90 <= my_lat1 <= 90:
            break
        else:
            print('\n WARNING: Value must be between -90 and 90 ')
    except ValueError:
        print('\n ERROR: Enter value between -90 and 90 ')

distance = input('\nDo you require the distance between two locations (y or n)? ')
if distance == 'y':
    my_long2 = float(input('For distance, enter second longitude in decimal degrees: '))
    my_lat2 = float(input('Enter second latitude in decimal degrees: '))

if distance == 'y':
    dbp = calc_geo_dist(my_long1, my_lat1, my_long2, my_lat2)
    print('The distance between the two points = {} km'.format(str(dbp)))
# lon1 = -2.3372222, lat1 = 48.8363889; lon2 = 77.0655556, lat 2 = 38.9213889. Should be 6181.63 km