/*
 *  photom: a package to perform photometric calibration
 *  Copyright (C) 2001  Michael William Richmond
 *
 *  Contact: Michael William Richmond
 *           Physics Department
 *           Rochester Institute of Technology
 *           85 Lomb Memorial Drive
 *           Rochester, NY  14623-5603
 *           E-mail: mwrsps@rit.edu
 *
 *  
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *  
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *  
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */


#include <stdio.h>
#include <math.h>
#include "photom_time.h"

	/* place the LST for the given Julian Date into the passed variable. 
	   uses the procedure 	
	   given by Duffett-Smith's "Practical Astronomy with your Pocket
	   Calculator", p. 17 and 20. The observatory longitude,
	   supplied in degrees west of Greenwich, 
	   must be converted to the form of "hours west of Greenwich".

	   the LST is set in the range (0.0-23.99). if all went well,
	   0 is returned; a return value of -1 indicates some error. */

int
get_lst(jd, longitude, lst)
double jd, longitude, *lst;
{
	double j, s, t, t0, u;

	/* 
	 * longitude must be converted from decimal degrees to 
	 *     hours West of Greenwich 
	 */
	longitude /= 15.0;

	j = floor(jd + 0.5) - 0.5;
	s = j - 2451545.0;
	t = s/36525.0;
	t0 = 6.697374558 + 2400.051336*t + 0.000025862*t*t;
	while (t0 > 24.0)
		t0 -= 24.0;
	while (t0 < 0.0)
		t0 += 24.0;
	u = (jd - j)*24.0;
	u *= 1.002737909;
	t0 += u;
	while (t0 > 24.0)
		t0 -= 24.0;
	while (t0 < 0.0)
		t0 += 24.0;


	t0 -= longitude;
	while (t0 > 24.0)
		t0 -= 24.0;
	while (t0 < 0.0)
		t0 += 24.0;

	*lst = t0;
	return(0);
}
	
	/* convert the given Local Sidereal Time on the given Julian Date to 
	   Universal Time. This routine always picks the smaller of the two 
	   possible UT values for LST values between (0h, 0m) and (0h, 4m). 

	   
	   If an error occurs, the routine returns -1.  Otherwise, the UT
	   in the range (0.0 - 23.9333) is placed in the passed parameter, 
	   and 0 is returned. */

int
lst2ut(jd, lst, longitude, ut)
double jd, lst, longitude, *ut;
{
	double j, s, t, t0, gst;
	
	/* longitude must be converted from degrees to hours West of Greenwich */
	longitude /= 15.0;

	/* first, convert LST to Greenwich mean ST */
	gst = lst + longitude;
	while (gst > 24.0)	
		gst -= 24.0;
	while (gst < 0.0)
		gst += 24.0;

	/* now convert from GMST to UT */
	j = floor(jd + 0.5) - 0.5;
	s = j - 2451545.0;
	t = s/36525.0;
	t0 = 6.697374558 + 2400.051336*t + 0.000025862*t*t;
	while (t0 > 24.0)
		t0 -= 24.0;
	while (t0 < 0.0)
		t0 += 24.0;
	
	gst -= t0;
	while (gst > 24.0)	
		gst -= 24.0;
	while (gst < 0.0)
		gst += 24.0;
	*ut = gst*0.9972695663;

	return(0);
}

	/* convert from UT (fractional hours) on the given Julian Date to GST 
	   (again, fractional hours). */

void
ut2gst(jd, ut, gst)
double jd, ut, *gst;
{
	double j, s, t, t0;

	j = floor(jd + 0.5) - 0.5;
	s = j - 2451545.0;
	t = s/36525.0;
	t0 = 6.697374558 + 2400.051336*t + 0.000025862*t*t;
	while (t0 > 24.0)
		t0 -= 24.0;
	while (t0 < 0.0)
		t0 += 24.0;
	t0 += ut*1.002737909;
	while (t0 > 24.0)
		t0 -= 24.0;
	while (t0 < 0.0)
		t0 += 24.0;
	*gst = t0;
}

	/* convert from UT (fractional hours) to Local Sidereal Time. return 0
	   if all went well, or -1 if there was a problem reading the longitude
	   from the CONFIG file. */

int
ut2lst(jd, ut, longitude, lst)
double jd, ut, longitude, *lst;
{
	double gst;

	/* longitude must be converted from degrees to hours West of Greenwich */
	longitude /= 15.0;

	ut2gst(jd, ut, &gst);
	*lst = gst - longitude;		/* longitude is in fractional hours already */
	if (*lst < 0.0)
		*lst += 24.0;

	return(0);
}



	/* convert from GST (fractional hours) on the given JD to UT
	   (also fractional hours). the output UT is in the range 
	   (0.0 - 23.9333); in case of ambiguity, the smaller UT value is
	   always chosen. */

void
gst2ut(jd, gst, ut)
double jd, gst, *ut;
{
	double j, s, t, t0;
	
	/* now convert from GST to UT */
	j = floor(jd + 0.5) - 0.5;
	s = j - 2451545.0;
	t = s/36525.0;
	t0 = 6.697374558 + 2400.051336*t + 0.000025862*t*t;
	while (t0 > 24.0)
		t0 -= 24.0;
	while (t0 < 0.0)
		t0 += 24.0;
	
	gst -= t0;
	while (gst > 24.0)	
		gst -= 24.0;
	while (gst < 0.0)
		gst += 24.0;
	*ut = gst*0.9972695663;
}