来自"NOCOW"

解析几何，以地球中心为原点建立三维坐标系，求出两点的坐标，再求出距离，然后求出夹角，再乘R。

注意：C++可能要使用long double

---

华丽的分割线---------------------------------------------------------------------------------------------------------------------

来个具体版的： 求球面上两点的距离。

  建立三维坐标，设o(p,q),即经度为p，纬度为q，则
    x=r*cos(q)*sin(p);
     y=r*sin(q);（需根据南北赋为正负）
     z=r*cos(q)*cos(p);（需根据东西赋为正负）
  然后求出两点的直线距离，即弦长，然后利用arcsin，求解即可，注意精度问题，0.005，当然貌似π的取值如果用常数定义的话，只要位数多就没问题。

挂下丑陋的但过去的代码

uses math;
const r=3437.5;
var k1,k2,k3,k4,q,code,k,i,j,m,n:longint;
    s,st:string;
    p,ox,oy,oz,ox2,oy2,oz2,ans,x1,x2,y1,y2,dis:extended;
  function change(dob:extended):extended;
  begin
    exit(dob*pi/180);
  end;
begin
  for i:=1 to 3 do
  readln;
  //y1
  readln(s);
  st:=copy(s,1,pos('^',s)-1); val(st,p,code);  y1:=y1+p;  delete(s,1,pos('^',s));
  st:=copy(s,1,pos('''',s)-1);  val(st,p,code);  p:=p/60;  y1:=y1+p;  delete(s,1,pos('''',s));
  st:=copy(s,1,pos('"',s)-1);  val(st,p,code);  p:=p/3600;  y1:=y1+p;  delete(s,1,pos('"',s));
  if s[2]='S' then k2:=-1;
  //x1
  readln(s);
  q:=pos(' ',s);
  st:=copy(s,q+1,pos('^',s)-q-1);  val(st,p,code);  x1:=x1+p; delete(s,1,pos('^',s));
  st:=copy(s,1,pos('''',s)-1);  val(st,p,code);  p:=p/60;x1:=x1+p;  delete(s,1,pos('''',s));
  st:=copy(s,1,pos('"',s)-1);  val(st,p,code); p:=p/3600;x1:=x1+p;  delete(s,1,pos('"',s));
  if s[2]='W' then k1:=-1;
  //
  readln;
  //y2
  readln(s);
  st:=copy(s,1,pos('^',s)-1);  val(st,p,code);  y2:=y2+p;  delete(s,1,pos('^',s));
  st:=copy(s,1,pos('''',s)-1);  val(st,p,code);  p:=p/60;  y2:=y2+p;  delete(s,1,pos('''',s));
  st:=copy(s,1,pos('"',s)-1);  val(st,p,code); p:=p/3600;  y2:=y2+p;  delete(s,1,pos('"',s));
  if s[2]='S' then k4:=-1;
  //x2
  readln(s);
  q:=pos(' ',s);
  st:=copy(s,q+1,pos('^',s)-q-1);  val(st,p,code);  x2:=x2+p;  delete(s,1,pos('^',s));
  st:=copy(s,1,pos('''',s)-1);  val(st,p,code);  p:=p/60;  x2:=x2+p;  delete(s,1,pos('''',s));
  st:=copy(s,1,pos('"',s)-1);  val(st,p,code);  p:=p/3600;  x2:=x2+p;  delete(s,1,pos('"',s));
  if s[2]='W' then k3:=-1;
  //
  readln;
  //
  x1:=change(x1); y1:=change(y1); x2:=change(x2); y2:=change(y2);
  ox:=r*cos(y1)*sin(x1); if k1=-1 then ox:=-ox;
  ox2:=r*cos(y2)*sin(x2); if k3=-1 then ox2:=-ox2;
  oy:=r*sin(y1); if k2=-1 then oy:=-oy;
  oy2:=r*sin(y2); if k4=-1 then oy2:=-oy2;
  oz:=r*cos(y1)*cos(x1); oz2:=r*cos(y2)*cos(x2);
  dis:=sqrt(sqr(ox-ox2)+sqr(oy-oy2)+sqr(oz-oz2));
  ans:=arcsin((dis/2)/r)*2*r;
  writeln('The distance to the iceberg: ',ans:0:2,' miles. ');
  if 100.00-ans>0.005 then writeln('DANGER!')
end.

by Brad-vegetable