Zips Near Zips in pl/pgsql (postgres)

zips_near_zips.sql

/*
Requirements: 
1. A table of zip codes that contain latitude and longitude columns.

   Mine is called zip_codes and has the following relative columns:
         city, state, zip, latitude, longitude

   These datasets are available for sale online.

This data seems to be free from http://geocoder.ibegin.com/downloads.php

2. plpgsql language installed. See http://www.postgresql.org/docs/8.4/static/xplang-install.html for details.
3. PostgreSQL 8.4

*/
CREATE OR REPLACE FUNCTION zips_near_zip (zip_code varchar(5), max_miles int) RETURNS TABLE(zip_code varchar(5), miles_away float) as $$
	DECLARE 
	  earth_radius_in_miles FLOAT := 3963.1676;
	  latitude1 float;
	  longitude1 float;
	BEGIN

	select into latitude1, longitude1
		    latitude,  longitude from zip_codes
	where zip_codes.zip = zip_code;
	
	--RAISE NOTICE 'latitude1: %', latitude1;
	--RAISE NOTICE 'longitude1: %', longitude1;
	
	return query 
	SELECT zip_codes.zip, 
		Case 
			When zip_code=zip_codes.zip then 0
			else --distance between two points using latitude and longitude as coordinates on a sphere
				earth_radius_in_miles * acos(sin(radians(latitude1))*sin(radians(latitude)) 
				+ cos(radians(latitude1))*cos(radians(latitude))*cos(radians(longitude1-longitude)))
		End distance_apart
	FROM zip_codes
	Where max_miles >
		(
		Case 
			When zip_code=zip_codes.zip then 0
			else --distance between two points using latitude and longitude as coordinates on a sphere
				earth_radius_in_miles * acos(sin(radians(latitude1))*sin(radians(latitude)) 
				+ cos(radians(latitude1))*cos(radians(latitude))*cos(radians(longitude1-longitude)))
		End
		)
	order by distance_apart;
	
/* 
 *  Original Javascript Implementation found online at:
 *  http://www.movable-type.co.uk/scripts/latlong.html
  
	lat1 = lat1.toRad(); 
	lat2 = lat2.toRad();
	var dLon = (lon2-lon1).toRad();

	var y = Math.sin(dLon) * Math.cos(lat2);
	var x = Math.cos(lat1)*Math.sin(lat2) -
          Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
	return Math.atan2(y, x).toBrng();
*/

	END;
$$ LANGUAGE plpgsql;

-- you can retrieve nearby zip codes like this:
-- select * from zips_near_zip('93065',10) order by zip_code;

-- because it returns a table you can join to your tables:
-- select * from zips_near_zip('85297',100) join addresses on zip_code=zip;