Calcule a distância entre dois pontos de longitude de latitude SQL

CREATE OR REPLACE FUNCTION calculate_distance(lat1 float, lon1 float, lat2 float, lon2 float, units varchar)
            RETURNS float AS $dist$
                DECLARE
                    dist float = 0;
                    radlat1 float;
                    radlat2 float;
                    theta float;
                    radtheta float;
                BEGIN
                    IF lat1 = lat2 OR lon1 = lon2
                        THEN RETURN dist;
                    ELSE
                        radlat1 = pi() * lat1 / 180;
                        radlat2 = pi() * lat2 / 180;
                        theta = lon1 - lon2;
                        radtheta = pi() * theta / 180;
                        dist = sin(radlat1) * sin(radlat2) + cos(radlat1) * cos(radlat2) * cos(radtheta);
                        IF dist > 1 THEN dist = 1; END IF;
                        dist = acos(dist);
                        dist = dist * 180 / pi();
                        dist = dist * 60 * 1.1515;
                        IF units = 'K' THEN dist = dist * 1.609344; END IF;
                        IF units = 'N' THEN dist = dist * 0.8684; END IF;
                        RETURN dist;
                    END IF;
                END;
            $dist$ LANGUAGE plpgsql;


//Usage:
SELECT calculate_distance(41.311081, 69.240562, 39.7681, 64.4556, 'M');
SELECT calculate_distance(41.311081, 69.240562, 39.7681, 64.4556, 'K');
SELECT calculate_distance(41.311081, 69.240562, 39.7681, 64.4556, 'N');
Mirik