	program TestEq2Hor
	implicit real*8 (a-h,o-z)

c Example of input data / Przykladowe dane
	fi = 0.7d0
	D  =-1.1d0
	t  = 5.23d0

c Direct transformation / Transformacja prosta
	call Eq2Hor(D,t,fi,Alt,Az)
	print'(1x,2f17.13,5x,a,f10.6)',D,t,'Alt ->',Alt

c Reverse transformation / Transformacja odwrotna
	call Hor2Eq(Alt,Az,fi,D,t)
	print'(1x,2f17.13,5x,a,f10.6)',D,t,' Az ->',Az

c Output of this program
c -1.1000000000000  5.2300000000000     Alt -> -0.414212
c -1.1000000000000  5.2300000000000      Az ->  5.838039
	end


	subroutine Eq2Hor(D,t,fi,Alt,Az)

c Converts equatorial (D,t - declination, hour angle) into
c horizontal (Alt,Az - altitude, azimuth) coordinates at
c geographical latitude fi [all angles expressed in radians]

c Przelicza wspolrzedne rownikowe (deklinacja D i kat godzinny t)
c na horyzontalne (wysokosc Alt i azymut Az) na szerokosci
c geograficznej fi [wszystkie katy sa wyrazone w radianach]

	implicit real*8 (a-h,o-z)
	data TwoPi/6.2831853071795864769d0/
	sinfi = dSIN(fi)
	cosfi = dSQRT(1 - sinfi*sinfi)
	Az = dATAN2(dSIN(t),sinfi*dCOS(t) - dTAN(D)*cosfi)
	Alt= dASIN(dSIN(D)*sinfi + dCOS(D)*cosfi*dCOS(t))
	Az = dMOD(Az + TwoPi,TwoPi)
	end


	subroutine Hor2Eq(Alt,Az,fi,D,t)

c Converts horizontal (Alt,Az - altitude, azimuth) into equatorial 
c (D,t - declination, hour angle) coordinates at geographical
c latitude fi [all angles expressed in radians]

c Przelicza wspolrzedne horyzontalne (wysokosc Alt i azymut Az)
c na rownikowe (deklinacja D i kat godzinny t) na szerokosci
c geograficznej fi [wszystkie katy sa wyrazone w radianach]

	implicit real*8 (a-h,o-z)
	call Eq2Hor(-Alt,Az,fi,DNeg,t)
	D = -DNeg
	end
