Sorry for that long delay with this review. Had real difficult weeks of my life :(

Now it seems to be over, so getting back to normal.

Below is my opinions of the code i see in the branch

bb-10.1-anel-MDEV-13467-gis-feature-st_sphere-v2 

Hope it's the recent version.

+  Adapter for functions that takes two or three arguments.
+*/
+
+class Create_func_arg2_or_3 : public Create_func
+{

Do you think we need the separate Create_fund_arg2_or_3 class?
If you don't see any other class to inherit from it,
i'd recommend to leave the Create_func_distance_sphere alone.
No need to have two constructors in the Item_func_sphere_distance
too. Just use the list argument as it's done for the
Item_func_json_contains (in 10.2).


Anyway i belive we should get rid of duplicating code
about the 'param1' and 'param2'. They can be popped/checked
with no if (arg_count == 2) condition.

+  double distance= 0.0;
+  null_value= (args[0]->null_value || args[1]->null_value);
+  if (null_value)
+  {
+    // Returned item must be set to Null
+    DBUG_ASSERT(maybe_null);
   don't think it's needed.
+    return 0;
Should be return 0.0 i belive. Or goto handle_errors.
Thet goes to all but one 'return' statements in this function.

+  }

> ...
> if (!arg1 || !arg2)

This last condition should never happen since we have not-null_value.
So i'd remove that 'if' section completely.
And anyway it should do the 'return' in that branch.

+    null_value= args[2]->null_value;
+    if (null_value)
+    {
+      return 0;
+    }

seems nicer to me this way
    if (args[2]->null_value)
    {
      null_value= true;
      goto handle_errors;
    }

+  if (!(g1= Geometry::construct(&buffer1, arg1->ptr(), arg1->length())) ||
+      !(g2= Geometry::construct(&buffer2, arg2->ptr(), arg2->length())))
+  {
+     goto handle_errors;

we should launch the ER_GIS_INVALID_DATA here.

+    // Generate error message in case different geometry is used?
Yes, i think the explainig message here is good to have.

+double Item_func_sphere_distance::spherical_distance(Geometry *g1,
+                                                     Geometry *g2,
+                                                     const double sphere_radius)
+double Item_func_sphere_distance::spherical_distance_points(Geometry *g1,
+                                                            Geometry *g2,
+                                                            const double r)

Why these two are member static functions, not simply static?
If you plan any 'external' use for them, please add comments about
what exacly they do.

+double Item_func_sphere_distance::spherical_distance(Geometry *g1,
+                                                     Geometry *g2,
+                                                     const double sphere_radius)
+{
+  double res= 0.0;
+  switch (g1->get_class_info()->m_type_id)
+  {
+    case Geometry::wkb_point:
+      res= spherical_distance_points(g1, g2, sphere_radius);
+      break;
+
+    case Geometry::wkb_multipoint:
+      res= spherical_distance_points(g1, g2, sphere_radius);
+      break;

Why do we need this function?
We confirmed already that the type is either wkb_point or wkb_multipolygon,
so just can call spherical_distance_points. No?

>>>>>


+class Item_func_sphere_distance: public Item_real_func
+{
+  double sphere_radius= 6370986.0; // Default radius equals Earth radius

Some compilers don't like this syntax, and we didn't use it yet, so
i'd recomment to move this to the constructor or even better to the ::val_real().

+double Gis_point::calculate_haversine(const Geometry *g,
+                                            const double sphere_radius)
+{
+  DBUG_ASSERT(sphere_radius > 0);
+  double x1r= 0.0;
+  double x2r= 0.0;
+  double y1r= 0.0;
+  double y2r= 0.0;

That lot of unnecessary initializations seems excessive.

+  if (g->get_class_info()->m_type_id == Geometry::wkb_multipoint)
+  {
+    const char point_size= 4 + WKB_HEADER_SIZE + POINT_DATA_SIZE+1; //1 for the type
+    char point_temp[point_size];
+    memset(point_temp+4, Geometry::wkb_point, 1);
+    memcpy(point_temp+5, static_cast<const Gis_multi_point *>(g)->get_data_ptr()+5, 4);
+    memcpy(point_temp+4+WKB_HEADER_SIZE, g->get_data_ptr()+4+WKB_HEADER_SIZE,
+           POINT_DATA_SIZE);
+    point_temp[point_size-1]= '\0';
+    Geometry_buffer gbuff;
+    Geometry *gg= Geometry::construct(&gbuff, point_temp, point_size-1);

There is the Gis_multi_point::geometry_n() function, should be used here to get the point.

I don't like that lot of static_cast-s.
Isn't it nicer to have just the function of
double calculate_harvesine_points(Geometry *p1, Geometry *p2, double radius)
double calculate_harvesine_multipoint_point(Geometry *mp, Geometry *p, double radius)
double calculate_harvesine_multipoint_multipoint(Geometry *mp1, Geometry *mp2, double radius)
 which would call each other?

The ::get_xy_radian() method looks unneeded to me.
It's enough to have the to_radian(double d) function.

BTW it's possible to do the multipoint/multipoint calculations faster than n^2.
Not necessary to do it right now, just not to remember.

Best regards.

HF