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

static float r2d = 180 / M_PI;
static float azimuth(float lat, float dec, float lha);
static float altitude(float lat, float dec, float lha);
static float dd2dm(float angle);

int main(int argc, char *argv[])
{

	float alt, azi, lat, start_lat, stop_lat, dec, start_dec, stop_dec, lha;
	int i, int_alt, frac_alt;

	/* start_lat, stop_lat, start_dec, stop_dec, */

	sscanf(argv[1], "%f", &start_lat);
	sscanf(argv[2], "%f", &stop_lat);
	sscanf(argv[3], "%f", &start_dec);
	sscanf(argv[4], "%f", &stop_dec);


	/* print file header information to make this a Rich Text Format (.rtf) document */
	printf("{\\rtf1\\ansi {\\fonttbl1 {\\f0 Courier}}\n");
	/* set paper size to A4 and set the print margins */
	printf("\\paperw11909 \\paperh16834 \\margl800 \\margr200 \\margt250 \\margb200 \n");


	for (lat = start_lat; lat <= stop_lat; lat = lat + 1.0)
	{

		/* print page header in 10 point bold */
		printf("{\\f0\\b\\fs20\\sl0 LATITUDE = %3.2f\\b0\\par\n", dd2dm(lat));
		/* rest in 6 point with reduced line spacing so that it all fits on 1 page of A4 */
		printf("\\fs12\\sl-120 DEC = ");

		/* print column headers */
		for (dec = start_dec; dec <= stop_dec; dec = dec + 0.5)
		{ printf("%5.2f      ", dd2dm(dec)); }

		printf("\\par\n");
		printf("LHA\\par\n");


		for (lha = 0; lha <= 109; lha = lha + 1)
		{

			printf("%3.0f", lha);

			for (dec = start_dec; dec <= stop_dec; dec = dec + 0.5)
			{

				alt = r2d * altitude(lat / r2d, dec / r2d, lha / r2d);
				azi = r2d * azimuth(lat / r2d, dec / r2d, lha / r2d);

				/* convert alt to degrees and minutes */
				alt = dd2dm(alt);

				/* round azi to nearest degree */
				azi = floor(azi + 0.49);

				if (alt <= 0)
					printf("           ");
				else
					printf("%7.2f%4.0f", alt, azi);
			}

			printf("\\par\n");
			if (remainder(lha + 1, 5) == 0) {printf("\\par\n");}

		}
		printf("\\page}\n");
	}
    printf("}\n");
    return 0;
}

float azimuth(float lat, float dec, float lha)
{
	float x, y, azi;
	x = sin(lha);
	y = tan(dec) * cos(lat) - sin(lat) * cos(lha);
	azi = atan2(-x,y);
	if (azi < 0) {azi = azi + 2 * M_PI;}
	return azi;
}

float altitude(float lat, float dec, float lha)
{
	float x, alt;
	x = sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(lha);
	alt = asin(x);
	return alt;
}

static float dd2dm(float angle)
{
	int int_angle, frac_angle, sgn;

	sgn = 1.0;
	if (angle < 0) {angle = -angle; sgn = -1;}
	int_angle = (int)floor(angle);
	frac_angle = (int)floor( (angle - int_angle) * 60 + 0.49);
	if (frac_angle == 60) {frac_angle = 0; int_angle++;}
	return ( sgn * ((float)int_angle + (float)frac_angle / 100 ));
}