根据坐标计算两点距离的Oracle函数

CREATE OR REPLACE FUNCTION ITS_KK.GPSDistance(lon1 in FLOAT, lat1 in FLOAT, lon2 in FLOAT,
            lat2 in FLOAT) RETURN FLOAT
AS  v_addr FLOAT;
    a1     FLOAT;
    b      FLOAT;
    f      FLOAT;
    rlat1  FLOAT;
    rlat2  FLOAT;
    L      FLOAT;
    U1     FLOAT;
    U2     FLOAT;
    sinU1  FLOAT;
    cosU1  FLOAT;
    sinU2  FLOAT;
    cosU2  FLOAT;
    lambda FLOAT;
    lambdap FLOAT;
    iterLimit INT;
    sinLambda FLOAT;
    cosLambda FLOAT;
    sinSigma  FLOAT;
    cosSigma  FLOAT;
    Sigma     FLOAT;
    sinAlpha  FLOAT;
    cosSqAlpha FLOAT;
    cos2SigmaM FLOAT;
    c         FLOAT;
   
   
    uSq       FLOAT;
    aA        FLOAT;
    bB        FLOAT;
    deltaSigma FLOAT;
    s         FLOAT;
begin

   a1 := 6378137.0;
   b := 6356752.3142;
   f :=  1 / 298.257223563;
   rlat1 := Rad(lat1);
   rlat2 := Rad(lat2);
   L := Rad(lon2-lon1);
   U1 := atan((1-f)*tan(rlat1));
   U2 := atan((1-f)*tan(rlat2));
   sinU1 := sin(U1);
   cosU1 := cos(U1);
   sinU2 := sin(U2);
   cosU2 := cos(U2);
   lambda := L;
   lambdap := 0.0;
   iterLimit := 100;
  
   while(abs(lambda - lambdap)>0.000000000001 and iterLimit > 0) LOOP
       sinLambda := sin(lambda) ;
       cosLambda := cos(lambda) ;
       sinSigma  := sqrt(power((cosU2 * sinLambda), 2)+power((cosU1 * sinU2 - sinU1 * cosU2 * cosLambda), 2));
          if(sinSigma=0.0)
               then return 0.0;
          else
           cosSigma := sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
           Sigma := atan2(sinSigma, cosSigma);
           sinAlpha := cosU1 * cosU2 * sinLambda / sinSigma ;
           cosSqAlpha := 1.0 - power(sinAlpha, 2);
           cos2SigmaM := cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha  ;
           select nvl(cos2SigmaM,0) into cos2SigmaM from dual;
          
           C := f / 16.0 * cosSqAlpha * (4.0 + f * (4.0 - 3.0 * cosSqAlpha)) ;
           lambdap := lambda;
           lambda := L + (1 - C) * f * sinAlpha * (Sigma + C * sinSigma * (cos2SigmaM +
           C * cosSigma * (-1 + 2.0 * power(cos2SigmaM , 2)))) ;
           iterLimit := iterLimit - 1  ;
           
       end if;
   end LOOP;
  
    if(iterLimit = 0)
  then return 111;
     else
         uSq := cosSqAlpha*(power(a1,2)-power(b,2))/(power(b,2));
         aA  := 1 + uSq / 16384.0 * (4096.0 + uSq *
              (-768.0 + uSq * (320.0 - 175.0 * uSq)));
         bB  := uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq)))  ;
         deltaSigma := bB * sinSigma * (cos2SigmaM + bB / 4.0 * (cosSigma
          * (-1 + 2 * power(cos2SigmaM,2)) -bB / 6.0 * cos2SigmaM
          * (-3 + 4 * power(sinSigma , 2)) * (-3 + 4 * power(cos2SigmaM , 2))));
         
         s := b * aA * (Sigma - deltaSigma) ;
        
        
         return s;
  end if;
 
END GPSDistance;
/

 

 

 

 

 

CREATE OR REPLACE FUNCTION ITS_KK.RAD(lon in NUMBER) RETURN NUMBER
AS  v_addr NUMBER;
begin
  select acos(-1)*lon/180.0 into v_addr from dual ;
 
  return v_addr;
END RAD;
/