/* @(#)coordinate.c	17.1.1.1 (ES0-DMD) 01/25/02 17:50:05 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  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., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

#include <math.h>
double de2ra = 0.017453293;
double twopi = 6.28318531;
double dphi = -29.25694444;
double fmod();
void h2a(ra,dec,lst,daz,dal)
double ra,lst,dec,*daz,*dal;
{
double hha,ha,de,phi,sinphi,cosphi,tanphi,sinde,cosde;
double sinha,cosha,sinal,al,cosal;
double tazx,tazy,az,sinaz,cosaz;
hha = lst/3600. -ra/15.;
ha = hha*15.*de2ra;
de = dec*de2ra;

phi = dphi * de2ra;
sinphi = sin(phi);
cosphi = cos(phi);
tanphi = sinphi/cosphi;
sinde = sin(de);
cosde = cos(de);

sinha = sin(ha);
cosha = cos(ha);

sinal = sinphi * sinde + cosde * cosphi * cosha;
al = asin(sinal);
al = fmod(al,twopi);
if (al<= 0.) al += twopi;
cosal = cos(al);
*dal = al/de2ra;


tazx = cosde * sinha;
tazy = (sinphi*cosde*cosha-cosphi*sinde);
az = atan2(tazx,tazy);
az = fmod(az,twopi);
if (az <= 0.) az +=twopi;
cosaz = cos(az);
sinaz = sin(az);
*daz = az/de2ra;
}