
-- This script defines several stored functions to demonstrate how that can be done.

-- I want to create a view that, for each observation, lists the observed star and the altitude
-- the star was above the horizon at the time it was observed. Researchers may want to filter
-- out observations done at too low an altitude since such observations might be unreliable.

-- The following three functions compute the local sidereal time from the location of an
-- observer, the date of an observation, and the universal time when that observation is made.
-- Reference: "Astronomical Algorithms", second edition, by Jean Meeus, copyright 1998,
-- published by Willmann-Bell, ISBN=0-9943396-61-1.

-- Compute the Julian day number for 0h UTC on the Gregorian date given.
CREATE FUNCTION julian_day(
    incoming_year IN INTEGER, incoming_month IN INTEGER, incoming_day IN INTEGER) RETURNS FLOAT AS $$
DECLARE
  the_year  INTEGER := incoming_year;
  the_month INTEGER := incoming_month;
  the_day   INTEGER := incoming_day;
  a   INTEGER;
  b   INTEGER;
BEGIN
  IF the_month <= 2 THEN
    the_year  := the_year - 1;
    the_month := the_month + 12;
  END IF;
  a := the_year/100;
  b := 2 - a + (a/4);
  RETURN floor(365.25 * (the_year + 4716)) + floor(30.6001 * (the_month + 1)) + the_day + B - 1524.5;
END;
$$ LANGUAGE plpgsql;


-- Compute the sidereal time at the prime meridian for the date/time given. The return value is
-- given as an angle in degrees (24 hours = 360 degrees).
CREATE FUNCTION sidereal_time(datetime IN TIMESTAMP WITH TIME ZONE) RETURNS FLOAT AS $$
DECLARE
  t       FLOAT;
  theta_0 FLOAT;
  theta   FLOAT;
  the_year   INTEGER := EXTRACT(YEAR   FROM datetime);
  the_month  INTEGER := EXTRACT(MONTH  FROM datetime);
  the_day    INTEGER := EXTRACT(DAY    FROM datetime);
  the_hour   INTEGER := EXTRACT(HOUR   FROM datetime);
  the_minute INTEGER := EXTRACT(MINUTE FROM datetime);
  the_second INTEGER := EXTRACT(SECOND FROM datetime);
BEGIN
  t := (julian_day(the_year, the_month, the_day) - 2451545.0)/36525.0;
  theta_0 := 100.46061837 + 36000.770053608*t + 0.000387933*t*t + (t*t*t)/38710000.0;
  theta := 1.00273790935*360.0*((the_second + 60.0*(the_minute + 60.0*the_hour))/86400);
  theta := theta + theta_0;
  WHILE theta > 360.0 LOOP
    theta := theta - 360.0;
  END LOOP;
  RETURN theta;
END;
$$ LANGUAGE plpgsql;


-- Compute the sidereal time at a specific longitude on the earth. The return value is given in
-- degrees (24 hours = 360 degrees).
CREATE FUNCTION local_sidereal_time(
    datetime IN TIMESTAMP WITH TIME ZONE, longitude IN FLOAT) RETURNS FLOAT AS $$
DECLARE
  prime_time FLOAT := sidereal_time(datetime);
BEGIN
  RETURN prime_time + longitude;
END;
$$ LANGUAGE plpgsql;


-- Compute the altitude of a spot on the sky at a place and time. The return value is given in
-- degrees.
CREATE FUNCTION altitude(
        -- The time...
        datetime IN TIMESTAMP WITH TIME ZONE,

        -- The place...
        longitude IN FLOAT, latitude IN FLOAT,

        -- The position on the sky...
        ra IN FLOAT, declination IN FLOAT) RETURNS FLOAT AS $$
DECLARE
  LST float := local_sidereal_time(datetime, longitude);
  right_ascension FLOAT := 360.0 * (ra/24.0);

  pole_zenith FLOAT;  -- Distance from pole to zenith.
  pole_star   FLOAT;  -- Distance from pole to star.
  zenith_star FLOAT;  -- Distance from zenith to star.
  RA_diff     FLOAT;  -- Different between RA and local sidereal time.
BEGIN
  -- This computation involves spherical trigonometry (rule of cosines).
  pole_zenith := 90.0 - latitude;
  pole_star   := 90.0 - declination;
  RA_diff     := right_ascension - LST;

  -- Convert to radians.
  pole_zenith := (pole_zenith/360.0) * 2.0 * pi();
  pole_star   := (pole_star/360.0) * 2.0 * pi();
  RA_diff     := (RA_diff/360.0) * 2.0 * pi();

  -- Compute the interesting angle.
  zenith_star := acos(cos(pole_zenith) * cos(pole_star) +
                      sin(pole_zenith) * sin(pole_star) * cos(RA_diff));

  -- Convert back to degrees.
  zenith_star := 360.0 * zenith_star/(2.0 * pi());
  
  RETURN 90.0 - zenith_star;
END;
$$ LANGUAGE plpgsql;


-- For debugging purposes. Let us evalute the functions and see if the results look okay.
SELECT id, julian_day(CAST(EXTRACT(YEAR  FROM datetime) AS INTEGER),
                      CAST(EXTRACT(MONTH FROM datetime) AS INTEGER),
                      CAST(EXTRACT(DAY   FROM datetime) AS INTEGER)), sidereal_time(datetime)
FROM observation;


-- Create the desired view. Researchers can join this view against the observation table to
-- look up other details about an observation once it passes altitude criteria.
CREATE VIEW star_altitude AS
  SELECT ob.id as observation_id, ob.name, ob.const, altitude(
    ob.datetime, l.longitude, l.latitude, s.ra, s.dec) AS altitude
  FROM observation ob
    JOIN location l ON ob.location_id = l.id
    JOIN star s ON (ob.name = s.name AND ob.const = s.const);


-- Show the entire view just so we can see what it looks like.
SELECT * FROM star_altitude;
   
-- Show the entire view again, but pretty up the output a bit.
SELECT observation_id, name, const, round(CAST(altitude AS NUMERIC), 2) AS altitude FROM star_altitude;

-- This demonstrates joining the view to the observation table.
SELECT *
FROM star_altitude JOIN observation ON observation_id = observation.id;
