calculate two point distance on earth

create_distance_sphere.sql

-- calculate two point distance on earth
-- ref:
-- [1] https://www.scribd.com/presentation/2569355/Geo-Distance-Search-with-MySQL
-- [2] http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates
-- [3] https://en.wikipedia.org/wiki/Circle_of_a_sphere
-- [4] https://en.wikipedia.org/wiki/Great-circle_distance
-- [5] https://www.bram.us/2018/03/01/mysql-st_distance_sphere-polyfill/
-- [6] https://stackoverflow.com/questions/35939853/st-distance-sphere-in-mysql-not-giving-accurate-distance-between-two-locations
-- [7] https://www.movable-type.co.uk/scripts/latlong.html
-- [8] https://dev.mysql.com/doc/refman/5.7/en/spatial-convenience-functions.html#function_st-distance-sphere
-- [9] https://www.oc.nps.edu/oc2902w/general/mapmodel.pdf

delimiter $$
drop function if exists ST_Distance_Sphere_R$$
-- Point('longitude', 'latitude')
create function ST_Distance_Sphere_R (pt1 point, pt2 point, R double)
       returns double
       comment '(pt1, pt2, R); pt(longitude, latitude)'
       deterministic
       begin
        declare phi1 double;
    	declare phi2 double;
	declare dphi double;
	declare dlamda double;
	declare a double;
	declare c double;
	
	set phi1 = radians(y(pt1));
	set phi2 = radians(y(pt2));
	set dphi = radians(y(pt1) - abs(y(pt2)));
	set dlamda = radians(x(pt1) - x(pt2));
	set a = power(sin(dphi / 2), 2) + cos(phi1)*cos(phi2)*power(sin(dlamda / 2), 2);
	set c = 2 * asin(sqrt(a));
	/* set c = 2 * atan2(sqrt(a), sqrt(a)); */
	return R * c;
      end$$

drop function if exists ST_Distance_Sphere$$
create function ST_Distance_Sphere (pt1 point, pt2 point)
       returns double
       comment '(pt1, pt2); pt(longitude, latitude); r=6370986'
       deterministic
       begin
	return ST_Distance_Sphere_R(pt1, pt2, 6370986);
       end$$
       
drop function if exists ST_Distance_Sphere_Vincenty_R$$
create function ST_Distance_Sphere_Vincenty_R (pt1 point, pt2 point, R double)
       returns double
       comment '(pt1, pt2, R); pt(longitude, latitude)'
       deterministic
       begin
        declare phi1 double;
    	declare phi2 double;
	declare dphi double;
	declare dlamda double;
	declare a double;
	declare b double;
	declare c double;

	declare sin_phi1 double;
	declare cos_phi1 double;
	declare sin_phi2 double;
	declare cos_phi2 double;
	declare cos_dlamda double;
	declare sin_dlamda double;

	set phi1 = radians(y(pt1));
	set phi2 = radians(y(pt2));
	set dphi = radians(y(pt1) - y(pt2));
	set dlamda = radians(x(pt1) - x(pt2));

	set sin_phi1 = sin(phi1);
	set cos_phi1 = cos(phi1);
	set sin_phi2 = sin(phi2);
	set cos_phi2 = cos(phi2);
	set cos_dlamda = cos(dlamda);
	set sin_dlamda = sin(dlamda);
	
	set a = sqrt(power(cos_phi2 * sin_dlamda, 2) + power(cos_phi1*sin_phi2 - sin_phi1*cos_phi2*cos_dlamda, 2));
	set b = sin_phi1 * sin_phi2 + cos_phi1 * cos_phi2 * cos_dlamda;
	set c = atan2(a, b);
	/* set c = 2 * atan2(sqrt(a), sqrt(a)); */
	return R * c;
      end$$
      
drop function if exists ST_Distance_Sphere_Vincenty$$
create function ST_Distance_Sphere_Vincenty (pt1 point, pt2 point)
       returns double
       comment '(pt1, pt2); pt(longitude, latitude); r=6370986'
       deterministic
       begin
	return ST_Distance_Sphere_Vincenty_R(pt1, pt2, 6370986);
       end$$

delimiter ;