расстояния

gistfile1.sql

                                                                                                                                                                   
      CREATE OR REPLACE FUNCTION distance(point1 point, point2 point) RETURNS float8 AS $$                                                                                                     
      DECLARE                                                                                                                                                                                  
        long1   float8;                                                                                                                                                                        
        long2   float8;                                                                                                                                                                        
        lat1    float8;                                                                                                                                                                        
        lat2    float8;                                                                                                                                                                        
        cl1     float8;                                                                                                                                                                        
        cl2     float8;                                                                                                                                                                        
        sl1     float8;                                                                                                                                                                        
        sl2     float8;                                                                                                                                                                        
        delta   float8;                                                                                                                                                                        
        cdelta  float8;                                                                                                                                                                        
        sdelta  float8;                                                                                                                                                                        
        y       float8;                                                                                                                                                                        
        x       float8;                                                                                                                                                                        
        ad      float8;                                                                                                                                                                        
        dist    float8;                                                                                                                                                                        
        rad     float8;                                                                                                                                                                        
      BEGIN                                                                                                                                                                                    
        rad := 6372795;                                                                                                                                                                        
        long1 := point1[0] * pi() / 180;                                                                                                                                                       
        long2 := point2[0] * pi() / 180;                                                                                                                                                       
                                                                                                                                                                                               
        lat1 := point1[1] * pi() / 180;                                                                                                                                                        
        lat2 := point2[1] * pi() / 180;                                                                                                                                                        
                                                                                                                                                                                               
        cl1 := cos(lat1);                                                                                                                                                                      
        cl2 := cos(lat2);                                                                                                                                                                      
        sl1 := sin(lat1);                                                                                                                                                                      
        sl2 := sin(lat2);                                                                                                                                                                      
        delta := long2 - long1;                                                                                                                                                                
        cdelta := cos(delta);                                                                                                                                                                  
        sdelta := sin(delta);                                                                                                                                                                  
                                                                                                                                                                                               
        y := ((cl2 * sdelta) ^ 2 + (cl1 * sl2 -sl1 * cl2 * cdelta) ^ 2) ^ 0.5;                                                                                                                 
        x := sl1 * sl2 + cl1 * cl2 * cdelta;                                                                                                                                                   
        ad := atan2(y,x);                                                                                                                                                                      
        dist := ad * rad;                                                                                                                                                                      
                                                                                                                                                                                               
        RETURN dist;                                                                                                                                                                           
      END;                                                                                                                                                                                     
      $$ LANGUAGE plpgsql;