/******************************************************************************* ** PLANET: A program to determine the position of the Sun, Mercury, Venus, ** ** Mars, Jupiter, Saturn and Uranus. ** ** ** ** Reference Jean Meeus's book, ``Astronomical Formulae for Calculators'' ** ** ** ** Program by ** ** ** ** Fred T. Mendenhall ** ** ihnp4!inuxe!fred ** ** ** ** Copyright 1985 by Fred T. Mendenhall ** ** All rights reserved ** ** ** ** (This copyright prevents you from selling the program - ** ** the author grants anyone the right to copy and install the program on any ** ** machine it will run on) ** ******************************************************************************** ** Modified by Alan Paeth, awpaeth@watcgl, January 1987 for ** ** (1) non-interactive mode (uses current GMT) ** ** (2) output in simplified Yale format for use with star charter software ** ** (3) lines 365-380 rewritten as a simplier C expression (April '87) ** ******************************************************************************** ** Modified by Petri Launiainen, pl@intrin.FI, February 1987 for ** ** easier LOGFILE definition/modification in Makefile ** ******************************************************************************** ** Modified by Alan Paeth, September 1987, for command line flags. ** ** The program continues to provide reduced Yale output in the old style, ** ** but new versions of `starchart' accept either. The old format allows ** ** representation of magnitudes below -0.99 (useful for planets), but no ** ** spectral or stellar identification, which is not useful for planets. ** ******************************************************************************** ** Modified by Jim Cobb, March 1988, to compute moon positions as well. ** ** Also, to change slightly the computation of epli, the obliquity ** ** of the ecliptic. Formerly epli had a term ** ** ** ** +0.00256*cos(omeg*rad). ** ** ** ** This term is a compensation for computing the apparent position of the ** ** sun (Chapter 18 in Meeus). But for converting ecliptical coordinates ** ** to right ascension and declination we want the mean value ** ** (that is, we don't want this term). ** ******************************************************************************** ** Modified by Bob Leivian, Sept 1989, to compute local azimuth and altitude ** ** given a lat/lon, also fixed time so you can replace a single item ** ** without having to supply all (i.e., for today at 9:00 put -t 9) ** ** also ran the code through 'cb' to standardize the indention. ** ** I split the code in to more managible pieces. ** ** Also printed out time in local and UTC as well as JD. ** ** ** ** See "Time Notes" section below if you plan to rewrite the user interface. ** ******************************************************************************** ** Modified by Andreas Scherer, July 1992, obviously redoing a lot of Bob's ** ** work, in that the new AMIGA C-compilers are dealing time according to ** ** the ANSI standard (Jan 1st 1970!). At this place I included the ** ** JULDAY routine and some specialities of the SAS/C-compiler. ** ** Also I cleaned up the source code for better reading. Doesn't it look ** ** nicer that way? I didn't use a code beautifier! ** ** Most output texts now have their German counterparts. You have to ** ** #define EUTSCH, which gives the command line option `-DEUTSCH'. ** ** I also included some improvements, mostly concerning Cee, but some ** ** formulas needed a bit of speed up, like most of the polynomials, which ** ** are now evaluated by the Horner algorithm. ** *******************************************************************************/ #include #include #include /******************************************************************************* ** We need the Longitude, the Latitude, and the Time Zone. ** ** Default latitude longitude (of Tempe, AZ) -- (-111.94, 33.32, 7 MST). ** ** Replace with your favorite place. ** ** Now it's Lichtenfels, Bavaria. ** *******************************************************************************/ #define DEFLONG 11.1 #define DEFLAT 50.1 #define DEFZONE (-1.) /* CET */ #define RAD 0.01745329251994329651 /* PI/180 */ #define JD1900 2415020. /* Julian date for January 1st, 1900, 12h */ #define JD2000 2451545. /* Julian date for January 1st, 2000, 12h */ #define JJ100 36525. /* Days in a Julian century */ #ifndef PLANETFILE # define PLANETFILE "planet.star" #endif double longitude = DEFLONG; double latitude = DEFLAT; double GHA_Aries; int zone; FILE *logfile; char *progname; /******************************************************************************* ** Times Notes: ** ** ** ** This software would ideally use a generic UNIX system call to which would ** ** provide day and date info as integers, plus current time west of GMT. ** ** These values would become the default settings for each of the switches ** ** ** ** -m -d -y -t -z. ** ** ** ** Unfortunately, time on across many UNIX systems does not (to my knowledge) ** ** exist as a uniform subroutine call. (With ANSI they should. Try it!) ** ** ** ** Instead, we presume that the very low level "gettimeofday" (GMT in ** ** seconds, from 1 Jan 1970) exists on all UNIX installations, and use this ** ** subroutine to use the present time of day when no command line parameters ** ** are present, as defaults for the current time have not been filled in. ** ** ** ** When any flag-style command line parameters are present, hardcoded ** ** defaults are substitude for -m, -d, -y -t switches, which may then be ** ** overridden. The -z flag is filled in by a companion return value from ** ** "gettimeofday". This approach makes for poor (on account of hard coding) ** ** defaults. ** ** ** ** A worthwhile change would be to simplify the command line control ** ** structure and default mechanism by using a higher level system call, or ** ** (safer), write the routine to convert Unix GMT into day, month, year, so ** ** that it lives within this program. This would allow for more reasonable ** ** defaults, and still merely one UNIX system call to get GMT time. Perhaps ** ** sources from a UNIX system may be studied to do this correctly. ** ** ** ** Expertise in "time" on a specific UNIX system is *NOT* what is required ** ** -- what *IS* required is knowledge and access to enough independent UNIX ** ** systems to insure (to a high probability) that the code remains portable. ** ** In other words, only low-level system calls known to work across a large ** ** range of systems should be employed. ** *******************************************************************************/ /******************************************************************************* ** For getting the current GMT or UTC. ** *******************************************************************************/ #include /******************************************************************************* ** This string is especially for the SAS/C-compiler. The first three chars ** ** are the name of your local time zone, the next three represent the normal ** ** difference from UTC to your local time, and the last three characters ** ** determine, whether the so called `daylight saving time' is momentarily in ** ** effect, which means that one hour more is substracted from the UTC value. ** *******************************************************************************/ char *_TZ = "CET-02CDT"; double current_time(m,d,y,t,z) int *m,*d,*y,*z; double *t; { double julday(); int to_int(); #ifdef AMIGA /******************************************************************************* ** The old version didn't work with the new ANSI-C-compiler. time doesn't ** ** return the seconds from January 1st, 1978 (happy birthday, AMIGA) but ** ** the seconds from January 1st, 1970. We use it here solely for localtime. ** ** The julday-routine does all the computing. ** ** The first instruction sets up a handful of variables in `time.h', one of ** ** them is `timezone', the number of seconds, by which UTC is off. This is ** ** determined by the `_TZ'-string, initialized above. ** *******************************************************************************/ time_t amiga_now; struct tm *p; tzset(); amiga_now = time(&amiga_now); p = gmtime(&amiga_now); *y = p->tm_year + 1900; /* local year */ *m = p->tm_mon + 1; /* local month */ *d = p->tm_mday; /* local day of month */ *t = p->tm_hour + p->tm_min/60.; /* local time (0.-23.999) */ *z = to_int(timezone / 60.); /* zone difference in min */ return(julday(*d,*m,*y,*t)); #endif /* AMIGA */ #ifdef UNIX /******************************************************************************* ** I don't know whether the time-routine works in the ANSI-standard. ** ** It better should. Just in case it doesn't I hold the old version. ** *******************************************************************************/ struct timeval tv; struct timezone tz; struct tm *p; gettimeofday(&tv,&tz); p = localtime(&tv); *y = p->tm_year + 1900; /* local year */ *m = p->tm_mon + 1; /* local month */ *d = p->tm_mday; /* local day of month */ *t = p->tm_hour + p->tm_min/60.; /* local time (0.-23.999) */ *z = tz.tz_minuteswest; /* local time west of Greenwich */ return(julday(*d,*m,*y,*t + *z)); #endif /* UNIX */ #ifdef NO_TIME_AVAIL /******************************************************************************* ** No time is available, pick a default, i. e., Noon January 1st, 1989 ** ** Due to JULDAY you may use any other set of values, maybe your birthday? ** *******************************************************************************/ *y = 1989; /* local year */ *m = 1; /* local month */ *d = 1; /* local day of month */ *t = 12.; /* local time */ *z = to_int(DEFZONE * 60.); /* zone difference west of Greenwich */ return(julday(*d,*m,*y,*t + DEFZONE); #endif /* NO_TIME_AVAIL */ } /******************************************************************************* ** The user is allowed to change the default time settings from his system by ** ** some switches. If some values are changed, they are checked and then they ** ** may replace the original values. ** *******************************************************************************/ double commandtime(argc,argv) int argc; char **argv; { double time,timez,now; int mm,dd,yy,j; double htod(),current_time(),julday(); int to_int(); void die(); #ifdef EUTSCH static char *usage = "\nAufruf: planet {Zeit}" "\noder\tplanet Julianisches Datum" "\noder\tplanet [-t hh.mm -z Zone -y Jahr" " -m Monat -d Tag -lo Laenge -la Breite]"; #else static char *usage = "\nusage:planet {current time used}" "\nor\tplanet juliandate" "\nor\tplanet [-t hh.mm -z zone -y year" " -m mon -d day -lo longitude -la latitude]"; #endif /******************************************************************************* ** Get the current (system) time for defaults. ** ** You're getting the fraction of Julian days and the time zone in hours. ** *******************************************************************************/ now = current_time(&mm,&dd,&yy,&time,&zone); timez = zone / 60.; /******************************************************************************* ** No args - use current time. ** *******************************************************************************/ if(argc == 1) return(now); /******************************************************************************* ** One numeric arg - use it as the Julian date. ** *******************************************************************************/ if((argv[1][0] != '-')) { if(argc == 2) return(atof(argv[1])); else /******************************************************************************* ** Too many args without switches. ** *******************************************************************************/ #ifdef EUTSCH die("Keine Schalter - %s",usage); #else die("no switches - %s",usage); #endif } /******************************************************************************* ** Flags present, set up defaults, process any user overrides. ** *******************************************************************************/ else { for(j=1; j 12) #ifdef EUTSCH die("Zeitzone: GMT +- 12h (oestlich negativ)\n"); #else die("time zone: UTC +- 12h (east negativ)\n"); #endif zone = timez * 60.; break; /******************************************************************************* ** Change the current year. ** *******************************************************************************/ case 'y': yy = atoi(argv[++j]); break; /******************************************************************************* ** Change the current month of the year. ** *******************************************************************************/ case 'm': mm = atoi(argv[++j]); if(mm < 1 || mm > 12) #ifdef EUTSCH die("Monat zwischen 1 und 12\n"); #else die("month must be 1..12\n"); #endif break; /******************************************************************************* ** Change the current day. ** *******************************************************************************/ case 'd': dd = atoi(argv[++j]); if(dd < 1 || dd > 31) #ifdef EUTSCH die("Tag zwischen 1 und 31\n"); #else die("day must be 1..31\n"); #endif break; /******************************************************************************* ** Change the position of observation. ** *******************************************************************************/ case 'l': switch(argv[j][2]) { /******************************************************************************* ** Change the default longitude of the observer's position. ** *******************************************************************************/ case 'o': longitude = atof(argv[++j]); if(longitude < -180 || longitude > 180) #ifdef EUTSCH die("Laenge zwischen -180 und +180\n"); #else die("longitude must be +-180\n"); #endif break; /******************************************************************************* ** Change the default latitude of the observer's position. ** *******************************************************************************/ case 'a': latitude = atof(argv[++j]); if(latitude < -90 || latitude > 90) #ifdef EUTSCH die("Breite zwischen -90 und +90\n"); #else die("latitude must be +-90\n"); #endif break; /******************************************************************************* ** Tell me what you want. ** *******************************************************************************/ default: #ifdef EUTSCH die("Breite und Laenge angeben\n"); #else die("specify lat or lon\n"); #endif } break; /******************************************************************************* ** Something went wrong. ** *******************************************************************************/ default: #ifdef EUTSCH die("Unbekannter Schalter - %s",usage); #else die("unknown switch - %s",usage); #endif } /******************************************************************************* ** ??? ** *******************************************************************************/ if(j == argc) #ifdef EUTSCH die("Schalter zuviel - %s",argv[j-1]); #else die("trailing command line flag - %s",argv[j - 1]); #endif } } return(julday(dd,mm,yy,time + timez)); } /******************************************************************************* ** It's much easier to `clear the screen' than most people think. ** ** There might be some trouble on IBM's, but who cares? ** *******************************************************************************/ void cls() { putchar('\f'); } /******************************************************************************* ** Compute the value of a polynomial with degree 3 at time t. ** *******************************************************************************/ double poly(a0,a1,a2,a3,t) double a0,a1,a2,a3,t; { return(a0 + (a1 + (a2 + a3 * t) * t) * t); } /******************************************************************************* ** Return the integer part of a floating point number. ** *******************************************************************************/ int to_int(z) double z; { return((int)z); } int to_long(z) double z; { return((long)z); } double kepler(e,M) double e,M; { double corr=1.,E0=M,E1; while(corr > 0.000001) { corr = (M + e/RAD * sin(E0*RAD) - E0) / (1 - e * cos(E0*RAD)); E1 = E0 + corr; if(corr < 0) corr *= -1.; E0 = E1; } return(E1); } double truean(e,E) double e,E; { double nu; nu = 2. * atan(sqrt((1 + e) / (1 - e)) * tan(E * RAD / 2)) / RAD; if(nu < 0.) nu += 360.; if(fabs(nu - E) > 90.) if(nu > 180.) nu -= 180.; else nu += 180.; return(nu); } double longi(w2,i,u) double w2,i,u; { double x,y,l; y = cos(i*RAD) * sin(u*RAD); x = cos(u*RAD); l = atan2(y,x)/RAD; if(l<0.) l += 360.; return(l + w2); } double lati(u,i) double u,i; { return(asin(sin(u*RAD) * sin(i*RAD)) / RAD); } void speak(ra,dec,dis,mag,sym,str) double ra,dec,dis,mag; char *sym,*str; { double lha,sha,altitude,azimuth; int h,m,s,deg; double range(); int to_int(); if(ra < 0) ra += 360.; sha = 360. - ra; /******************************************************************************* ** Convert from degs to hours. ** *******************************************************************************/ ra /= 15.; h = to_int(ra); m = to_int( (ra - h) * 60 ); s = to_int(((ra - h) * 60 - m) * 60); printf(" %-12s%2dh %02dm %02ds ",str,h,m,s); if(logfile) fprintf(logfile,"%02d%02d%02d",h,m,s); deg = to_int(dec); m = to_int( (dec - deg) * 60 ); s = to_int(((dec - deg) * 60 - m) * 60); if(m < 0) m *= -1; if(s < 0) s *= -1; #ifdef AMIGA /******************************************************************************* ** Its printf doesn't know about signed. Well, I didn't try it yet. ** ***************************************************************************** */ if(deg > 0) printf(" +%2ddeg %02dm %02ds ",deg,m,s); else printf(" %3ddeg %02dm %02ds ",deg,m,s); #else printf(" %+3ddeg %02dm %02ds ",deg,m,s); #endif /* AMIGA */ if(dis > 0.0001) printf(" %6.3f",dis); else printf(" "); if(logfile) if(mag < 0.) fprintf(logfile,"%+03d%02d-%2d%s %s\n", deg,m,-to_int(10 * mag),sym,str); else fprintf(logfile,"%+03d%02d%3d%s %s\n", deg,m,to_int(100 * mag),sym,str); /******************************************************************************* ** If we have a valid angle. ** *******************************************************************************/ if(GHA_Aries > 0) { /******************************************************************************* ** Now tell where it is at a given place on earth. ** *******************************************************************************/ lha = range(GHA_Aries + sha + longitude); altitude = sin(latitude * RAD) * sin(dec * RAD) + cos(latitude * RAD) * cos(dec * RAD) * cos(lha * RAD); altitude = asin(altitude) / RAD; azimuth = sin(lha * RAD) / (cos(lha * RAD) * sin(latitude * RAD) - tan(dec * RAD) * cos(latitude * RAD)); /******************************************************************************* ** Correct for quadrant. ** *******************************************************************************/ if(lha <= 180.) { if(azimuth > 0) azimuth = 180. + atan(azimuth) / RAD; else azimuth = 360. + atan(azimuth) / RAD; } else { if(azimuth <= 0.) azimuth = 180. + atan(azimuth) / RAD; else azimuth = atan(azimuth) / RAD; } printf(" %5.1f %5.1f\n",altitude,azimuth); } else printf("\n"); } /******************************************************************************* ** Helio_trans converts heliocentric ecliptical coordinates to geocentric ** ** right ascension and declination. It also outputs the converted value. ** *******************************************************************************/ void helio_trans(r,b,ll,Stheta,Sr,epli,mag,sym,str) double r,b,ll,Stheta,Sr,epli,mag; char *sym,*str; { double N,D,lambda,delta,beta,rcosb=r*cos(b*RAD),rsinb=r*sin(b*RAD); double range(); void geo_trans(); N = rcosb * sin((ll - Stheta)*RAD); D = rcosb * cos((ll - Stheta)*RAD) + Sr; lambda = atan2(N,D) / RAD; if(lambda<0.) lambda += 360.; lambda = range(lambda + Stheta); delta = sqrt(N*N + D*D + rsinb * rsinb); beta = asin(rsinb / delta) / RAD; geo_trans(lambda,beta,epli,delta,mag,sym,str); } /******************************************************************************* ** geo_trans converts geocentric ecliptical coordinates to geocentric right ** ** ascension and declination. It also outputs the converted value. Note ** ** that the output coordinates are referred to the mean equinox of date. If ** ** these coordinates are to be used in conjunction with a star database, use ** ** the epoch program to convert the planet's coordinates to the epoch for the ** ** star database. ** *******************************************************************************/ void geo_trans(lambda,beta,epli,delta,mag,sym,str) double lambda,beta,epli,delta,mag; char *sym,*str; { double N,D,RA,DEC; double sinlam=sin(lambda*RAD),cosepl=cos(epli*RAD),sinepl=sin(epli*RAD); void speak(); N = sinlam * cosepl - tan(beta*RAD) * sinepl; D = cos(lambda*RAD); RA = atan2(N,D)/RAD; DEC = asin(sin(beta*RAD)*cosepl + cos(beta*RAD)*sinepl*sinlam)/RAD; speak(RA,DEC,delta,mag,sym,str); } /******************************************************************************* ** htod(str) reads floating point strings of the form {+-}hh.{mm} thus ** ** allowing for fractional values in base 60 (i.e., degrees/minutes). ** *******************************************************************************/ double htod(s) char *s; { double sign; int full,frac; char *t; void die(); t = s-1; while(*++t) { if(((*t<'0') || (*t>'9')) && (*t!='.') && (*t!='+') && (*t!='-')) #ifdef EUTSCH die("Falsches Zeichen im numerischen Argument tt.mm"); #else die("non-digit in dd.mm style numeric argument"); #endif } if(s == NULL) return(0.); full = frac = 0; sign = 1.; if(*s == '-') { sign = -1.; s++; } else if(*s == '+') s++; while(*s && *s != '.') full = 10*full + *s++ - '0'; if(*s++ == '.') { if(*s) frac = 10*(*s++ - '0'); if(*s) frac += *s++ - '0'; if(frac>59) frac = 59; } return(sign * ((double)full + ((double)frac) / 60.)); } /******************************************************************************* ** In case of emergency you need a safety exit. ** *******************************************************************************/ void die(a,b) char *a,*b; { fprintf(stderr,"%s: ",progname); fprintf(stderr,a,b); fprintf(stderr,"\n"); exit(1); } /******************************************************************************* ** The main routine. I placed it to the end. ** *******************************************************************************/ void main(argc,argv) int argc; char **argv; { double jd,T0,epli,temp; long int_jd; int y,m,d,h,min,lh,lmin; char *str; double planet_pos(),sidereal_time(),commandtime(); int to_int(); void cls(),speak(),die(),moon_pos(),caldat(),caltim(); progname = argv[0]; if((logfile = fopen(PLANETFILE,"w")) == 0) #ifdef EUTSCH die("Fehler beim Oeffnen der Logdatei %s\n",PLANETFILE); #else die("fail on opening logfile %s\n",PLANETFILE); #endif /******************************************************************************* ** Calculate Julian date in fractions of days. ** *******************************************************************************/ jd = commandtime(argc,argv); /******************************************************************************* ** Convert the Julian date to calendar day. ** *******************************************************************************/ int_jd = to_int(jd +.5001); caldat(int_jd,&m,&d,&y); /******************************************************************************* ** Julian day to UTC. ** *******************************************************************************/ temp = jd - .4999; caltim(temp - floor(temp),&h,&min); cls(); #ifdef EUTSCH printf("\n Am %d.%02d.%4d um %d Uhr %02d GMT (Julianisch: %.3f)\n", d,m,y,h,min,jd); #else printf("\n At %d/%d %d %d:%02d UTC (Julian: %.3f)\n", m,d,y,h,min,jd); #endif /******************************************************************************* ** UTC to Local. ** *******************************************************************************/ temp -= (zone / 1440.); caltim(temp - floor(temp),&lh,&lmin); int_jd = to_int(jd + .5001 - (zone / 1440.)); caldat(int_jd,&m,&d,&y); #ifdef EUTSCH str = "keine Zeitangabe"; #else if(lh > 12) { lh -= 12; str = "pm"; } else { str = "am"; if(lh == 0) lh = 12; } #endif #ifdef EUTSCH printf(" Lokal %d.%02d. %d Uhr %02d (", d,m,lh,lmin); #else printf(" local %d/%d %d:%02d%s (", m,d,lh,lmin,str); #endif /******************************************************************************* ** Compute exact local time for a given longitude. ** *******************************************************************************/ temp = jd - 0.5 + (longitude / 360.); caltim(temp - floor(temp),&lh,&lmin); GHA_Aries = 15. * sidereal_time(jd); #ifdef EUTSCH printf("%d:%02d wahr lokal bei Brt: %1.1f Lng: %1.1f, GHA Aries: %1.1f)\n", lh,lmin,latitude,longitude,GHA_Aries); #else printf("%d:%02d true local at lat:%1.1f lon:%1.1f, GHA Aries:%1.1f)\n", lh,lmin,latitude,longitude,GHA_Aries); #endif T0 = (jd - JD1900) / JJ100; #ifdef EUTSCH printf("\n Objekt Rektaszension Deklination Abstand Höhe Azimuth\n"); #else printf("\n OBJECT RA DEC DISTANCE altitude azimuth\n"); #endif epli = planet_pos(jd,T0); moon_pos(T0,epli); putchar('\n'); if(logfile) close(logfile); logfile = 0; /******************************************************************************* ** Show local position of some stars also. ** *******************************************************************************/ speak(101.1714,-16.7013,0.,-1.46," ","Sirius"); #ifdef EUTSCH speak(213.7956,19.2360,0.,-0.04," ","Arktur"); speak( 88.6508, 7.4057,0., 0.50," ","Beteigeuze"); #else speak(213.7956,19.2360,0.,-0.04," ","Arcturus"); speak( 88.6508, 7.4057,0., 0.50," ","Betelgeuse"); #endif putchar('\n'); } /******************************************************************************* ** End of the main program and end of this file. ** *******************************************************************************/