This is dms.c in view mode; [Download] [Up]
/********************************************************/
/* dms.c */
/* */
/* Utility programs for unit conversions and calendar */
/********************************************************/
/***** description
*
* $Id: dms.c,v 1.7 1993/04/30 18:19:54 craig Exp $
*
*/
/***** modification history
*
* $Log: dms.c,v $
* Revision 1.7 1993/04/30 18:19:54 craig
* Changed the output calls from sending the output to stdout to
* sending the output to a named output file (which can be stdout).
*
* Revision 1.6 1993/04/27 18:18:42 craig
* corrected for an annoying rounding error in hms that would
* have 60.000 seconds displayed.
*
* Revision 1.5 1993/04/27 16:53:21 craig
* Changed showrd to not print if prtflg is not set.
*
* Changed showcor so that the position corrections are made even
* of the prtflg is not set.
*
* Revision 1.4 1993/04/22 19:42:55 craig
* minor changes.
*
* Revision 1.3 1993/04/21 21:12:56 craig
* Changed the path for the satellite.h include.
* Changed ecnsts to pcnsts.
*
* Revision 1.2 1993/04/21 15:13:09 craig
* First working version. Ran through indent and converted to ansi.
* Added hooks to work with the satellite programs.
*
*
*/
/***** include files *****/
#include <string.h>
#include <math.h>
#include "aaproto.h"
#include "satellite.h"
/***** global variables *****/
/* from cnstinit.c */
extern int prtflg;
extern FILE *outfile;
extern double TDT;
extern double UT;
extern struct MCONSTANTS mcnsts;
extern struct PCONSTANTS pcnsts;
#define RTOH (12.0/mcnsts.pi)
/***** local global variables *****/
static char *months[12] = {
"January",
"February",
"March",
"April",
"May",
"June",
"July",
"August",
"September",
"October",
"November",
"December"
};
static char *days[7] = {
"Sunday",
"Monday",
"Tuesday",
"Wednesday",
"Thursday",
"Friday",
"Saturday"
};
static int month = 1;
static long cyear = 1986;
static short yerend;
static double day = 1.0;
static int hours = 0;
static int minutes = 0;
static double seconds = 0.0;
/********************************************************/
/* showrd */
/* */
/* Display Right Ascension and Declination */
/* from input equatorial rectangular unit vector. */
/* Output vector pol[] contains R.A., Dec., and radius. */
/********************************************************/
int showrd (char *msg, double *p, double *pol)
{
int i;
double x, y, r;
r = 0.0;
for (i = 0; i < 3; i++)
{
x = p[i];
r += x * x;
}
r = sqrt (r);
x = matan2 (p[1], p[0]);
pol[0] = x;
y = asin (p[2] / r);
pol[1] = y;
pol[2] = r;
if (prtflg)
{
fprintf (outfile, "%s R.A. ", msg);
hms (outfile, x);
fprintf (outfile, "Dec. ");
dms (outfile, y);
fprintf (outfile, "\n");
}
return (0);
}
/*********************************************************/
/* showcor */
/* */
/* Display magnitude of correction vector in arc seconds */
/*********************************************************/
int showcor (char *strng, double *p, double *dp)
{
int i;
double p1[3], dr, dd;
for (i = 0; i < 3; i++)
{
p1[i] = p[i] + dp[i];
}
deltap (p, p1, &dr, &dd);
if (prtflg)
{
fprintf (outfile, "%s dRA %.3fs dDec %.2f\"\n",
strng, mcnsts.ra2sec * dr / 15.0, mcnsts.ra2sec * dd);
}
return (0);
}
/****************************************/
/* dms */
/* */
/* Radians to degrees, minutes, seconds */
/****************************************/
int dms (FILE *out, double x)
{
double s;
int d, m;
s = x * mcnsts.ra2de;
if (s < 0.0)
{
fprintf (out, " -");
s = -s;
}
else
{
fprintf (out, " ");
}
d = s;
s -= d;
s *= 60;
m = s;
s -= m;
s *= 60;
fprintf (out, "%3dd %02d\' %05.2f\" ", d, m, s);
return (0);
}
/****************************************/
/* hms */
/* */
/* Radians to hours, minutes, seconds */
/****************************************/
int hms (FILE *out, double x)
{
int h, m;
double s;
s = x * RTOH;
if (s < 0.0)
{
s += 24.0;
}
h = (int) s;
s -= (double) h;
s *= 60.;
m = (int) s;
s -= (double) m;
s *= 60.;
/* correct for rounding */
if (60. - s <= .0005)
{
m++;
s = 0.;
if (m == 60)
{
h++;
m = 0;
}
}
fprintf (out, "%3dh %02dm %06.3fs ", h, m, s);
return (0);
}
/* routines from julian.c
*
* julian.c calculates Julian day number from calendar
* date, and the date and day of the week from Julian day.
* The Julian date is double precision floating point
* with the origin used by astronomers. The calendar output
* converts fractions of a day into hours, minutes, and seconds
* There is no year 0. Enter B.C. years as negative; i.e.,
* 2 B.C. = -2.
*
* The approximate range of dates handled is 4713 B.C. to
* 54,078 A.D. This should be adequate for most applications.
*
* B.C. dates are calculated by extending the Gregorian sequence
* of leap years and century years into the past. This seems
* the only sensible definition, but I don't know if it is
* the official one.
*
* Note that the astronomical Julian day starts at noon on
* the previous calendar day. Thus at midnight in the morning
* of the present calendar day the Julian date ends in .5;
* it rolls over to tomorrow at noon today.
*
* The month finding algorithm is attributed to Meeus.
*
* - Steve Moshier
*/
/***********/
/* getdate */
/***********/
double getdate (void)
{
double J;
/* Get operator to type in a date. */
getlong ("Calendar date: Year", &cyear);
if ((cyear > 53994) || (cyear < -4713))
{
printf ("Year out of range.\n");
J = 0.0;
goto pdate;
}
if (cyear == 0)
{
printf ("There is no year 0.\n");
J = 0.0;
goto pdate;
}
getint ("Month (1-12)", &month);
getdouble ("Day.fraction", &day);
/* Find the Julian day. */
J = caltoj (cyear, month, day);
/* printf( "Julian day %.1f\n", J ); */
pdate:
/* Convert back to calendar date. */
/* jtocal( J ); */
return (J);
}
/*****************************************************/
/* caltoj */
/* */
/* Calculate Julian day from Gregorian calendar date */
/*****************************************************/
double caltoj (long year, int month, double day)
{
int BC;
long y, a, b, c, e, m;
double J;
BC = 0;
/* The origin should be chosen to be a century year that is also a
leap year. We pick 4801 B.C. */
y = year + 4800;
if (year < 0)
{
BC = 1;
y += 1;
}
/* The following magic arithmetic calculates a sequence whose
successive terms differ by the correct number of days per calendar
month. It starts at 122 = March; January and February come after
December. */
m = month;
if (m <= 2)
{
m += 12;
y -= 1;
}
e = (306 * (m + 1)) / 10;
a = y / 100; /* number of centuries */
if (year <= 1582L)
{
if (year == 1582L)
{
if (month < 10)
{
goto julius;
}
if (month > 10)
{
goto gregor;
}
if (day >= 15)
{
goto gregor;
}
}
julius:
fprintf (outfile, " Julian Calendar assumed.\n");
b = -38;
}
else
{ /* -number of century years that
are not leap years */
gregor:
b = (a / 4) - a;
}
c = ((long) pcnsts.dapcen * y) / 100; /* Julian calendar years and leap
years */
/* Add up these terms, plus offset from J 0 to 1 Jan 4801 B.C. Also
fudge for the 122 days from the month algorithm. */
J = b + c + e + day - 32167.5;
return (J);
}
/***************************************************/
/* jtocal */
/* */
/* Calculate month, day, and year from Julian date */
/***************************************************/
int jtocal (FILE *out, double J)
{
int month, day;
int BC;
long year, a, c, d, x, y, jd;
double dd;
if (J < 1721425.5) /* January 1.0, 1 A.D. */
{
BC = 1;
}
else
{
BC = 0;
}
jd = J + 0.5; /* round Julian date up to integer */
/* Find the number of Gregorian centuries since March 1, 4801 B.C. */
a = (100 * jd + 3204500L) / 3652425L;
/* Transform to Julian calendar by adding in Gregorian century years
that are not leap years. Subtract 97 days to shift origin of JD to
March 1. Add 122 days for magic arithmetic algorithm. Add four
years to ensure the first leap year is detected. */
c = jd + 1486;
if (jd >= 2299160.5)
{
c += a - a / 4;
}
else
{
c += 38;
}
/* Offset 122 days, which is where the magic arithmetic month formula
sequence starts (March 1 = 4 * 30.6 = 122.4). */
d = (100 * c - 12210L) / (long) pcnsts.dapcen;
/* Days in that many whole Julian years */
x = ((long) pcnsts.dapcen * d) / 100L;
/* Find month and day. */
y = ((c - x) * 100L) / 3061L;
day = c - x - ((306L * y) / 10L);
month = y - 1;
if (y > 13)
{
month -= 12;
}
/* Get the year right. */
year = d - 4715;
if (month > 2)
{
year -= 1;
}
/* Day of the week. */
a = (jd + 1) % 7;
/* Fractional part of day. */
dd = day + J - jd + 0.5;
/* post the year. */
cyear = year;
if (BC)
{
year = -year + 1;
cyear = -year;
fprintf (out, "%ld B.C. ", year);
}
else
{
fprintf (out, "%ld ", year);
}
day = dd;
fprintf (out, "%s %d %s", months[month - 1], day, days[a]);
/* Flag last or first day of year */
if (((month == 1) && (day == 1))
|| ((month == 12) && (day == 31)))
{
yerend = 1;
}
else
{
yerend = 0;
}
/* Display fraction of calendar day as clock time. */
a = dd;
dd = dd - a;
hms (out, 2.0 * mcnsts.pi * dd);
if (J == TDT)
{
fprintf (out, "TDT\n"); /* Indicate Terrestrial Dynamical Time */
}
else if (J == UT)
{
fprintf (out, "UT\n"); /* Universal Time */
}
else
{
fprintf (out, "\n");
}
return (0);
}
/********************************/
/* mod360 */
/* */
/* Reduce x modulo 360 degrees */
/********************************/
double mod360 (double x)
{
long k;
double y;
k = x / 360.0;
y = x - k * 360.0;
while (y < 0.0)
{
y += 360.0;
}
while (y > 360.0)
{
y -= 360.0;
}
return (y);
}
/************************/
/* modtp */
/* */
/* Reduce x modulo 2 pi */
/************************/
double modtp (double x)
{
double y;
y = floor (x / mcnsts.twopi);
y = x - y * mcnsts.twopi;
while (y < 0.0)
{
y += mcnsts.twopi;
}
while (y >= mcnsts.twopi)
{
y -= mcnsts.twopi;
}
return (y);
}
/********************************************************/
/* gethms */
/* */
/* Get operator to type in hours, minutes, and seconds */
/********************************************************/
double gethms (void)
{
double t;
getint ("Time: Hours", &hours);
getint ("Minutes", &minutes);
getdouble ("Seconds", &seconds);
t = (3600.0 * hours + 60.0 * minutes + seconds) / pcnsts.secpda;
return (t);
}
/***********/
/* getlong */
/***********/
int getlong (char *msg, long *num)
{
char s[40];
printf ("%s ( %ld ) ? ", msg, *num);
if (strlen (gets (s)))
{
sscanf (s, "%ld", num);
}
return (0);
}
/**********/
/* getint */
/**********/
int getint (char *msg, int *num)
{
char s[40];
printf ("%s ( %d ) ? ", msg, *num);
if (strlen (gets (s)))
{
sscanf (s, "%d", num);
}
return (0);
}
/*************/
/* getdouble */
/*************/
int getdouble (char *msg, double *num)
{
char s[40];
printf ("%s ( %f ) ? ", msg, *num);
if (strlen (gets (s)))
{
sscanf (s, "%lf", num);
}
return (0);
}
/*************/
/* getstring */
/*************/
int getstring (char *msg, char *num, int maxsize)
{
char s[40];
printf ("%s ( %s ) ? ", msg, num);
if (strlen (gets (s)))
{
strncpy (num, s, maxsize);
}
return (0);
}
/****************************************************************/
/* deltap */
/* */
/* Convert change in rectangular coordinatates to change */
/* in right ascension and declination. */
/* For changes greater than about 0.1 degree, the */
/* coordinates are converted directly to R.A. and Dec. */
/* and the results subtracted. For small changes, */
/* the change is calculated to first order by differentiating */
/* tan(R.A.) = y/x */
/* to obtain */
/* dR.A./cos**2(R.A.) = dy/x - y dx/x**2 */
/* where */
/* cos**2(R.A.) = 1/(1 + (y/x)**2). */
/* */
/* The change in declination arcsin(z/R) is */
/* d asin(u) = du/sqrt(1-u**2) */
/* where u = z/R. */
/* */
/* p0 is the initial object - earth vector and */
/* p1 is the vector after motion or aberration. */
/****************************************************************/
int deltap (double *p0, double *p1, double *dr, double *dd)
{
int i;
double dp[3], A, B, P, Q, x, y, z;
P = 0.0;
Q = 0.0;
z = 0.0;
for (i = 0; i < 3; i++)
{
x = p0[i];
y = p1[i];
P += x * x;
Q += y * y;
y = y - x;
dp[i] = y;
z += y * y;
}
A = sqrt (P);
B = sqrt (Q);
if ((A < 1.e-7) || (B < 1.e-7) || (z / (P + Q)) > 5.e-7)
{
P = matan2 (p0[1], p0[0]);
Q = matan2 (p1[1], p1[0]);
Q = Q - P;
while (Q < -mcnsts.pi)
{
Q += 2.0 * mcnsts.pi;
}
while (Q > mcnsts.pi)
{
Q -= 2.0 * mcnsts.pi;
}
*dr = Q;
P = asin (p0[2] / A);
Q = asin (p1[2] / B);
*dd = Q - P;
return (0);
}
x = p0[0];
y = p0[1];
if (x == 0.0)
{
*dr = 1.0e38;
}
else
{
Q = y / x;
Q = (dp[1] - dp[0] * y / x) / (x * (1.0 + Q * Q));
*dr = Q;
}
x = p0[2] / A;
P = sqrt (1.0 - x * x);
*dd = (p1[2] / B - x) / P;
return (0);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.