--BEGIN;

/*

 Outdated and incomplete!
 Go look at:
 	pgsphere http://gborg.postgresql.org/project/pgsphere/
	PostGIS http://postgis.refractions.net/


 Trigonometric functions for PostgreSQL
 by Hans Schou <chlor@schou.dk> 2002.
 License: GPL

 Calculation of nautic distance on the earth.

 Earth radius:
  6378km @ equator
  6357km @ poles
  mean  6367.5km

 Nautic mile: 1 minute
  1852m
 E.rad: 3438.17nm = 6367.5 / 1.852
*/

/*

Finding Distance Using the Earth's Grid System
http://mathforum.org/library/drmath/view/54945.html
Assuming the earth is a sphere and you want the great circle distance 
between two points:

Let (x,y,z) be a point on a sphere of radius a. The spherical 
coordinates of (x,y,z) are (a,phi,theta), where phi is like latitude, 
but it is measured from the positive z-axis. The angle phi varies 
between 0 and pi. The angle theta is like longitude, but is measured 
from the positive x-axis (towards the positive y-axis).  The angle 
theta varies between 0 and 2*pi.

If you are given the spherical coordinates (a,phi,theta) of a point, 
the (x,y,z) coordinates are given by

x = a*cos(theta)*sin(phi)
y = a*sin(theta)*sin(phi}
z = a*cos(phi).

So, if you have a point on the earth's surface, you know the number a, 
which is the radius of the earth, you know or can easily calculate phi 
from the latitude, and you know or can easily calculate theta from the 
longitude.

So, suppose you are given the points (a,phi_1,theta_1) and 
(a,phi_2,theta_2). The first thing to do is to calculate the 
rectangular coordinates r_1=(x_1,y_1,z_1) and r_2=(x_2,y_2,z_2) of 
these points, as above. After that, letting alpha be the angle between 
the lines joining the center (0,0,0) of the sphere to r_1 and r_2, use 
the dot product of the vectors r_1 and r_2 to calculate alpha.

You will find that 

alpha= arccos(r_1.r_2/a^2) = arccos(cos(phi_1)*cos(phi_2)
            +cos(theta_1-theta_2)*sin(phi_1)*sin(phi_2)).

The great circle distance between r_1 and r_2 is 

d = a*arccos(alpha).  

-Doctor Jerry,  The Math Forum
 Check out our web site!  http://mathforum.org/dr.math/
*/

drop sequence wrest_wID_seq;
drop table wrest;
create table wrest(
	wID		serial,
	name	text,		-- Name of wrest
	lat		float8,	-- Latitude in degrees decimal
	lon		float8,
	posX	float8,
	posY	float8,
	posZ	float8
);
--COMMENT ON TABLE wrest IS 'Diving Information database';

/*
 New York
  Latitude 40.75 north
  Longitude 74 west
  phi = (90 - 40.75)/180 * pi = 49.25/180 * pi = 0.8595
  the = (360 - 74)/180 * pi = 286/180 * pi = 4.991
 Copenhagen
  Latitude 56 north
  Longitude 13 east
  phi = (90 - 56)/180 * pi = 34/180 * pi = 0.5934119456781
  the = (360 - -13)/180 * pi = 373/180 * pi =  6.510078109939

  x = p * sin(phi) * cos(the)
  y = p * sin(phi) * sin(the)
  z = p * cos(phi)

  gam = ( x1*x2 + y1*y2 + z1*z2 ) / 1

  distance = gam * 3438.17mi
*/

DROP FUNCTION phi( float8 );
CREATE FUNCTION phi( float8 ) RETURNS FLOAT8 AS '
	DECLARE
		p_lat	ALIAS FOR $1;	-- Latitude
	BEGIN
		RETURN (90 - p_lat)/180 * 3.141567; --pi();
	END;
' language 'PLpgSQL';
 select phi( 56.0 );

DROP FUNCTION theta( float8 );
CREATE FUNCTION theta( float8 ) RETURNS FLOAT8 AS '
	DECLARE
		p_lon	ALIAS FOR $1;	-- Longitude
	BEGIN
		IF p_lon < 0 THEN
			RETURN (-1* (p_lon/180) * 3.141567); -- pi();
		ELSE
			RETURN (360 - p_lon)/180 * 3.141567; -- pi();
		END IF;
	END;
' language 'PLpgSQL';
 select theta( 13.0 );
 select theta( -13.0 );

DROP FUNCTION sphX( float8, float8 );
CREATE FUNCTION sphX( float8, float8 ) RETURNS FLOAT8 AS '
	DECLARE
		p_phi	ALIAS FOR $1;	-- phi - Latitude
		p_the	ALIAS FOR $2;	-- theta - Longitude
	BEGIN
		RETURN (sin(p_phi) * cos(p_the));
	END;
' language 'PLpgSQL';
select sphX( 0.0, 0.0 );
select sin(0.0);

DROP FUNCTION sphY( float8, float8 );
CREATE FUNCTION sphY( float8, float8 ) RETURNS FLOAT8 AS '
	DECLARE
		p_phi	ALIAS FOR $1;	-- phi - Latitude
		p_the	ALIAS FOR $2;	-- theta - Longitude
	BEGIN
		RETURN sin(p_phi) * sin(p_the);
	END;
' language 'PLpgSQL';

DROP FUNCTION sphZ( float8 );
CREATE FUNCTION sphZ( float8 ) RETURNS FLOAT8 AS '
	DECLARE
		p_phi	ALIAS FOR $1;	-- phi - Latitude
	BEGIN
		RETURN cos(p_phi);
	END;
' language 'PLpgSQL';

DROP FUNCTION GreatCircle( float8, float8, float8, float8, float8, float8 );
CREATE FUNCTION GreatCircle( float8, float8, float8, float8, float8, float8 ) RETURNS FLOAT8 AS '
	DECLARE
		p_x1	ALIAS FOR $1;
		p_y1	ALIAS FOR $2;
		p_z1	ALIAS FOR $3;
		p_x2	ALIAS FOR $4;
		p_y2	ALIAS FOR $5;
		p_z2	ALIAS FOR $6;
	BEGIN
		RETURN acos( p_x1*p_x2 + p_y1*p_y2 + p_z1*p_z2 );
	END;
' language 'PLpgSQL';

DROP FUNCTION GC( float8, float8, float8, float8 );
CREATE FUNCTION GC( float8, float8, float8, float8 ) RETURNS FLOAT8 AS '
	DECLARE
		p_lat1	ALIAS FOR $1;
		p_lon1	ALIAS FOR $2;
		p_lat2	ALIAS FOR $3;
		p_lon2	ALIAS FOR $4;
		phi1		float8;
		the1		float8;
		phi2		float8;
		the2		float8;
		sinphi1		float8;
		sinphi2		float8;
	BEGIN
		phi1 = phi(p_lat1);
		the1 = theta(p_lon1);
		phi2 = phi(p_lat2);
		the2 = theta(p_lon2);
		sinphi1 = sin(phi1);
		sinphi2 = sin(phi2);
		RETURN acos(
			 	sinphi1*cos(the1)*sinphi2*cos(the2)
			+ sinphi1*sin(the1)*sinphi2*sin(the2)
			+ cos(phi1)*cos(phi2)
			);
	END;
' language 'PLpgSQL';

-- select GC( 40.75, 74, 51.5, 0 ) * 3438.17 ;

DROP FUNCTION AddWrest( text, float8, float8 );
CREATE FUNCTION AddWrest( text, float8, float8 ) RETURNS BOOL AS '
	DECLARE
		p_name	ALIAS FOR $1;	-- Name of wrest
		p_lat		ALIAS FOR $2;
		p_lon		ALIAS FOR $3;
	BEGIN
		INSERT INTO wrest( name, lat, lon, posX, posY, posZ )
			VALUES( p_name, p_lat, p_lon,
				sphX( phi(p_lat), theta(p_lon) ),
				sphY( phi(p_lat), theta(p_lon) ),
				sphZ( phi(p_lat) )
			);
		RETURN ''T'';
	END;
' language 'PLpgSQL';

DROP FUNCTION DistID( int4, int4 );
CREATE FUNCTION DistID( int4, int4 ) RETURNS FLOAT8 AS '
	DECLARE
		p_id1	ALIAS FOR $1;
		p_id2	ALIAS FOR $2;
		p1 RECORD;
		p2 RECORD;
	BEGIN
		SELECT posX, posY, posZ INTO p1 FROM wrest WHERE wID = p_id1;
		SELECT posX, posY, posZ INTO p2 FROM wrest WHERE wID = p_id2;
		RETURN acos(
			p1.posX*p2.posX + p1.posY*p2.posY + p1.posZ*p2.posZ
		);
	END;
' language 'PLpgSQL';
 
SELECT AddWrest('Copenhagen', 56.00, -13.0 );
SELECT AddWrest('New York', 40.75, 74.0 );
SELECT AddWrest('London', 51.50, 0.0 );
SELECT AddWrest('0,0', 0.0, 0.0 );
SELECT AddWrest('Cape Town', -35.00, -19.0 );
SELECT AddWrest('Berlin', 53.00, -14.0 );
SELECT AddWrest('Stockholm', 58.00, -19.0 );
SELECT AddWrest('North pole', 90, 0.0 );
SELECT AddWrest('Tokyo', 35.00, -140.0 );
SELECT AddWrest('Wrangel island', 70.00, -180 );
SELECT AddWrest('Thule', 75.0, 75.0 );

--SELECT DistID( 1, 2 ) * 3438.17;
--SELECT DistID( 1, 3 ) * 3438.17;
--SELECT DistID( 3, 2 ) * 3438.17;

SELECT
  wID,
  name AS "From Copenhagen",
  trunc( DistID(1, wID)*3438.17 ) AS Distance
 FROM wrest
 WHERE DistID(1, wID) < 3000/3438.17
 ORDER BY DistID(1, wID)
 LIMIT 10
 OFFSET 1	-- Do not show own position
;

SELECT
  wID,
  name AS "From North pole",
  trunc( DistID(8, wID)*3438.17 ) AS Distance
 FROM wrest
 WHERE DistID(8, wID) < 4000/3438.17
 ORDER BY distance
 LIMIT 10
 OFFSET 1	-- Do not show own position
;

SELECT
  wID,
  name AS "From Stockholm",
  trunc( DistID(7, wID)*3438.17 ) AS "Distance nm"
 FROM wrest
 WHERE 3000/3438.17 > DistID(7,wID) -- "Distance nm"
 ORDER BY "Distance nm"
 LIMIT 10
 OFFSET 1	-- Do not show own position
;

COMMIT;