#ifdef	ournix
#include "ournix.h"
#endif
char sccsID[] = "@(#) satelite.c V1.1 Copyright Julian H. Stacey, Munich, 17th April 1989\n" ;

/* Calculates compass bearing and azimuth from a specific point on earth,
|  to point at a geo-stationary satellite in orbit above the equator.
| 
|  Written in C, runs on unix & msdos.
|
|  Equations, & orbit radius & Munich position from Christopher Sharp.
|  Spelling: this file is deliberately mis-spelt with just one 'l', to conform
|  to msdos file naming limitations.
|  
|  On unix compile as		cc satelite.c -lm
|  On msdos MSC V4 compile using cc.bat, that make no specific reference to
|	slibfp.obj as cl.exe includes it automatically.
|  
|  JJ Outstanding Work:
|	I have not checked whether trig functions are always passed parameters
|	within their operational range.
|	It currently works fine for Munich Germany,
|	But Australians or Americans , being in different quadrants,
|	may wish to check themselves.
|	Please email me to let me know if OK or not.
|
| Not needed here but useful info:
|	- earth axis is at 23.5 degree offset from 
|		right angle of line to sun.
|	- www.astrium-space.com running galileo say altitude 23616 Km
|	  from earth's surface
*/

#include	<stdio.h>
#include	<math.h>
extern	double atof() ;

/* Constants ---------------------------------- { */
#define	PYE		(4.0 * atan((double)1.0))	/* 3.14159 */
#define	DEGS_PER_RAD	(180.0 / PYE )
#define	RADS_PER_DEG	(PYE / 180.0 )

#define	SAT_RADIUS 42160.0 
				/* Geostationary orbit in km from centre of 
				earth (> 22000 miles above earth according to
				dixons leaflet */

#define	EARTH_RAD 6367.47 /* km earth radius (avge of polar & equatorial radii*/

/* I think a nautical mile is 1 minute of arc IE 
	2 * PYE * EARTH_RAD / ( 360 * 60 * 60 ) = 1.85 KM 
 	1 second of arc = 1/60th, ie ~30.9 metres.
 Munich Riem Airport Coordinates:
	Flughafenbezugspunkt Lage nach Koordinaten:
	48 degrees 07 minutes 54 seconds North,
	11 degrees 41 minutes 57 seconds East,
	Hoehe 528m ueber NN 		<-- not relevant but interesting
*/
#define	RIEM_LAT  (48.0 + ( 7.0/60.0) + (54.0/(60.0 * 60.0)))
		/* (54/60 + 7) / 60 + 48 = 48.131667 */
#define	RIEM_LONG (11.0 + (41.0/60.0) + (57.0/(60.0 * 60.0)))
		/* ((( 57 / 60 ) + 41 ) / 60 ) + 11 =  11.699167

/* City centre (north cathedral tower (Frauenturm) 48,8,23 North 11,34,28 East
 *	48 degrees, 8 minutes, 23 seconds North
 *        11 degrees, 34 minutes, 28 seconds East
 *        48.139722, 11.574444, height 530m over NN    
 * Sundials:
 *	12:00 Greenwich = 13:00 BST = 14:00 CEST.
 *	GMT zone zero extends from the Meridian to 15 degrees West.
 *	We would need to be a full 15 (=360/24) degrees East for a sundial
 *	to show 1 hour (1/24 day) ahead of Greenwich. 
 *	At 12:00 GMT when a German wrist watch in summer reports 14:00, 
 *	A sundial in Munich would show roughly 12:46    ( 60 x ( 11.5/15 ) )
 */
#define	CITY_LAT	48.139722
#define	CITY_LONG	11.574444
#define	HOME_LAT	CITY_LAT
#define	HOME_LONG	CITY_LONG

/* Marienplatz From Geoff GPS N 48,08,29.4	E 11,34,58.5
*/
#define MARIENPLATZ_LAT 	( 48.0 + 8/60 + 29.4/(60*60))
#define	MARIENPLATZ_LONG	( 11.0 + 34/6 + 58.5/(60*60))

/* Map of Holz Str 27d
	http://maps.google.com/?ie=UTF8&z=17&ll=48.129492,11.568764&spn=0.005757,0.008905&t=h&om=1 
*/

#define	ASTRA_LONG 19.2		/* East. Astra 1A satellite longitude */
/* ---------------------------------- } */

char **ARGV ;

wrong()
	{
	fprintf(stderr,
		"Syntax error, correct invocation is:\n  %s %s\n  (%s)\n",
		*ARGV,
		"[Earth_Latitude Earth_Longitude] [Satellite_Longitude] ",
     "+ve Latitude ==> N of equator, +ve longitude ==> E of Greenwich Meridian."
		) ;
		/* can either omit both earth co-ordinates,
		   or satellite longitude, or both sets of values,
		   in either of which cases, defaults are taken. */
	exit(1) ;
	}
			
main(argc,argv)
	int argc ; char **argv ;
	{
	double	earth_latitude ;	/* earth latitude +ve = N of equator */
	double	earth_longitude ;	/* earth longditude
					+ve = East of Greenwich meridian */
	double	sat_longitude ;		/* in degrees, 
					+ve = East of Greenwich meridian */
	double	dif_longitude ;
	double	azimuth,inclination ;
	double	x ;

	ARGV = argv ;
#ifdef	VSL	/* { */
#include	"../../include/vsl.h"
#endif		/* } */
	argc-- ;
	if ((argc == 2) || (argc > 3) ) wrong() ;

	/* Set default values */
	earth_latitude = HOME_LAT ;
	earth_longitude = HOME_LONG ;
	sat_longitude = ASTRA_LONG ;

	if ((argc == 3 ) || (argc==1))
		{
		if (argc == 3 )
			{
			earth_latitude = atof(*++argv) ;
			earth_longitude = atof(*++argv) ;
			}
		sat_longitude = atof(*++argv) ; 
			/* examples
			 "19.2" for 19.2 East 7 Astra inc. analogue
			 "28.2" for 28.2 East 3 Astra sats inc. digital
			  Note some of 19.2 now also transmits digital.
			So maybe its not necessary for me to put a new
			dish other end of flat ? 
		Run examples:
			satelite
			Inclination: 34.24 above horizon, Azimuth: 169.81.
			satelite 19.2
			Inclination: 34.24 above horizon, Azimuth: 169.81.
			satelite 28.2
			Inclination: 32.42 above horizon, Azimuth: 158.16.
			*/
		}

	if ( ( (earth_latitude > 90.0) || (earth_latitude < -90.0) ) ||
		( (earth_longitude > 180.0 ) || (earth_longitude < -180.0 ) ) ||
		( (sat_longitude > 180.0 ) || (sat_longitude < -180.0 ) ) )
		wrong() ;

	/* convert to radians */
	earth_latitude *= RADS_PER_DEG ;
	earth_longitude *= RADS_PER_DEG ;
	sat_longitude *= RADS_PER_DEG ;

	dif_longitude = sat_longitude - earth_longitude ;

	/* trig: cos(dif_longitude)==cos(dif_longitude) * cos(earth_latitude) */
	x = acos( cos(dif_longitude) * cos(earth_latitude) ) ;

	inclination = atan( ((SAT_RADIUS - EARTH_RAD)/(SAT_RADIUS + EARTH_RAD))
				/ tan(x/2.0) ) - (x/2.0) ;
	/* JJ lint reports for line above : division by 0. */
	inclination *= DEGS_PER_RAD ;
	printf("Inclination: %.2f %s horizon, ",
		inclination, (inclination > 0.0) ? "above" : "below");

	/* trig: sin(azimuth) == sin(dif_longitude) / sin(x) */
	azimuth = asin( (sin(dif_longitude)) / (sin(x)) ) ;
	azimuth *= DEGS_PER_RAD ;
	azimuth = 180 - azimuth ;
	printf("Azimuth: %.2f.\n", azimuth + 0.005 ) ;	
	exit(0) ;
	}