#ifdef ournix #include "ournix.h" #endif char sccsID[] = "@(#) satelite.c V1.3 Copyright Julian H. Stacey, Munich, 17th April 1989, 2015\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 #include #include 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 The old airport East of the city) 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 See also: http://www.berklix.com/~jhs/maps/#coords http://www.berklix.com/~jhs/src/bsd/fixes/FreeBSD/ports/gen/astro/xearth/files/ */ #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 530 m over NN * http://www.muenchen.de/volltextsuche.html?query=530%20m * 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 http://maps.google.com/?ie=UTF8&z=17&ll=48.129492,11.568764&spn=0.005757,0.008905&t=h&om=1 */ /* http://www.meteoblue.com/en/germany/weather-m%C3%BCnchen 48.14 N, 11.57 E, 519m asl */ #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,\nAzimuth (North=0, South East=+135): 169.81. satelite 19.2 Inclination: 34.24,\nAzimuth (North=0, South East=+135): 169.81. satelite 28.2 Inclination: 32.42,\nAzimuth (North=0, South East=+135): 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:\t\t\t\t%.2f %s horizon\n", 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 (North=0, South East=+135):\t%.2f\n", azimuth + 0.005 ) ; exit(0) ; }