Great Circles on Mercator's Chart

Manipulate[
Graphics[{map[[1]],Dynamic[phil[p1,p2]]},PlotLabel->Dynamic[Style["distance = "<>ToString[distance[p1,p2]]<>" statute miles","Label"]],ImageSize->{500,350}],
{{p1,{6949,-2008.52}},{-10780,-7000},{10780,7000},Locator},{{p2,{-2112,-343.5885}},{-10780,-7000},{10780,7000},Locator} ,
SynchronousInitialization->False,
Initialization:>{Quiet@Get["WorldPlot`"];
map=Quiet@WorldPlot`WorldPlot[{WorldPlot`World,RGBColor[.6,.5,.4]&},WorldPlot`WorldBackground->RGBColor[.7,.7,1],WorldPlot`WorldToGraphics->True,WorldPlot`WorldRange->{{-75,75},{-180,180}},WorldPlot`WorldProjection-> WorldPlot`Mercator];distance[p1_,p2_]:=Module[{lon1,lat1,lon2,lat2,places,L0,L1,L2,ϕ1,ϕ2,tani,sf,ϕL,a=3963.1906,points=20},
lon1=p1[[1]]180./10800.;
lon2=p2[[1]]180./10800.;
lat1=(-(π/2)+2 ArcTan[(E^(p1[[2]]0.00029342723004694836` ))])180./π;
lat2=(-(π/2)+2 ArcTan[(E^(p2[[2]]0.00029342723004694836` ))])180./π;
If[lat1==lat2,lat2=lat2+0.00029;];
places={{lon1,lat1},{lon2,lat2}};
places=Sort[places]2π/360.;
If[ϕ1==ϕ2,ϕ2=ϕ2+0.00029;];
L1=places[[2,1]];  ϕ1=places[[2,2]];
L2=places[[1,1]]; ϕ2=places[[1,2]];
If[L1==L2,L2=L2+0.00029;];
IntegerPart[a ArcCos[Cos[L1-L2] Cos[ϕ1] Cos[ϕ2]+Sin[ϕ1] Sin[ϕ2]]]];phil[p1_,p2_]:=Module[{lon1,lat1,lon2,lat2,places,L0,L1,L2,ϕ1,ϕ2,tani,sf,ϕL,a=3963.1906,points=20,distance},
lon1=p1[[1]]180./10800.;
lon2=p2[[1]]180./10800.;
lat1=(-(π/2)+2 ArcTan[(E^(p1[[2]]0.00029342723004694836` ))])180./π;
lat2=(-(π/2)+2 ArcTan[(E^(p2[[2]]0.00029342723004694836` ))])180./π;
If[lat1==lat2,lat2=lat2+0.00029;];
places={{lon1,lat1},{lon2,lat2}};
places=Sort[places]2π/360.;
L1=places[[2,1]];  ϕ1=places[[2,2]];
L2=places[[1,1]]; ϕ2=places[[1,2]];
If[ϕ1==ϕ2,ϕ2=ϕ2+0.00029;];
If[L1==L2,L2=L2+0.00029;];
L0=(1/2) (L1+L2-2 ArcCot[Cot[((L1-L2)/2)] Csc[ϕ1+ϕ2] Sin[ϕ1-ϕ2]]);
If[L0==L1,L1=L1+0.00029;];
tani=-Csc[L0-L1] Tan[ϕ1];
If[Abs[L1-L2]≤π,
ϕL=Table[{ArcTan[Sin[t-L0]tani],t},{t,Min[L1,L2],Max[L1,L2],Abs[Max[L1,L2]-Min[L1,L2]]/(points-1)}];,
ϕL=Table[{ArcTan[Sin[t-L0]tani],t},
{t,Min[L1,Sign[L1]π],Max[L1,Sign[L1]π],Abs[L1-Sign[L1]π]/(points/2-1)}];
ϕL=Join[ϕL,Table[{ArcTan[Sin[t-L0]tani],t},
{t,Min[-Sign[L1]π,L2],Max[-Sign[L1]π,L2],Abs[L2+Sign[L1]π]/(points/2-1)}]];
];
Graphics[WorldPlot`WorldGraphics[{Thickness[0.005],Red,Line[(21600./(2π))ϕL]},WorldPlot`WorldRange->{{-75,75},{-180,180}},WorldPlot`WorldProjection-> WorldPlot`Mercator,WorldPlot`WorldGridStyle->GrayLevel[.25]]][[1]]
];}]

Originally designed for use in nautical navigation, Mercator's projection is a map of the Earth's surface on which straight lines, called rhumb lines, represent a course of constant compass heading. On the surface of a sphere, however, such rhumb lines are not the shortest distance between two points. Rather, the shortest distance between two points is a great circle, whose shape is shown as the red line on this map. Move the locators to different points and the map will update to show the path of the great circle and the distance between the points in statute miles.

Contributed by: John Snyder

Additional contributions by: Jeff Bryant, Pratik Desai, and Carl Woll