/* * hjd - programme to calculate the JD, MJD, HJD dates * * (c) 2002 Richard Ogley. * * Version 0.1 06/02/2002 RNO. * * Does not rely on the Starlink heliocentric and barycentric position * and velocity routines. This is system independent. */ #include #include #include #include static double juldate(struct tm *greg_time); static double modjuldate(double jd); static double gethjd(double jd); static void heliocentric_ra_dec(double mjd, double *ra, double *dec); int main(void) { double jd; /* Julian date*/ double mjd; /* Modified Julian date */ double hjd; /* Heliocentric Julian date */ double hmjd; /* Heliocentric modified Julian date */ struct tm *ltime; /* Time and timezone structure */ time_t system_time; /* The system time in seconds from 1970 */ int year; /* The inputted year */ int month; /* The inputted month */ int day; /* The inputted day of the month */ int hour; /* The inputted hour */ int min; /* The inputted minutes */ int secs; /* The inputted seconds */ int early_date; /* A flag whether the date is before 1582 or not */ /* Defaults */ system_time = time(NULL); ltime = localtime(&system_time); if (ltime == NULL) { exit(EXIT_FAILURE); } /* Enter the time and date, or use the system time */ printf("\nEnter the time in the format: dd/mm/yyyy hh:mm:ss\n"); printf("Or enter a time of 0 to use the system time.\n"); (void) scanf("%d/%d/%d %d:%d:%d", &day, &month, &year, &hour, &min, &secs); /* Checks whether the user wants to input the system time */ if (day != 0) { ltime->tm_year = year - 1900; ltime->tm_mon = month - 1; ltime->tm_mday = day; ltime->tm_hour = hour; ltime->tm_min = min; ltime->tm_sec = secs; } /* Issue a warning if the date is before Oct 15 1582 */ early_date = 0; if (ltime->tm_year < -318) early_date = 1; else { if (ltime->tm_year == -318) { if (ltime->tm_mon < 9) early_date = 1; if (ltime->tm_mon == 9) { if (ltime->tm_mday < 16) early_date = 1; } } } if (early_date == 1) { printf("The Gregorian calender is inacurate before 15 Oct 1582.\n"); exit(EXIT_FAILURE); } /* Output the time and date specified, to check it is what the user wanted. */ printf("Time: %02i:%02i:%02i %i/%i/%i\n", ltime->tm_hour, ltime->tm_min, ltime->tm_sec, ltime->tm_mday, ltime->tm_mon+1, ltime->tm_year+1900); /* Calculate the Julian date */ jd = juldate(ltime); /* Calculate the modified Julian date */ mjd = modjuldate(jd); /* Calculate the modified heliocentric Julian date */ hmjd = gethjd(mjd); /* Calculate the heliocentric Julian date */ hjd = hmjd + 2400000.5; /* Print the results */ printf(" JD = %lf\n", jd); printf(" MJD = %f\n", mjd); printf("HMJD = %f\n", hmjd); printf(" HJD = %lf\n", hjd); exit(EXIT_SUCCESS); } /* * juldate * - procedure to calculate the Julian date * * Inputs * greg_time - a time structure of the form tm in * * Returns * Julian date * */ double juldate(struct tm *greg_time) { double date; /* A running variable calculating the Julian date */ int month; /* The month in the year */ int day; /* The day in the month */ int year; /* The years from 1900 */ /* Get the day, month and year from the greg_time structure and form them into 1/1/2000 rather than 1/0/100 from greg_time. */ day = greg_time->tm_mday; month = greg_time->tm_mon + 1; year = greg_time->tm_year + 1900; /* Calculate the Julian days */ date = day - 32076 + 1461*(year + 4800 + (month - 14)/12)/4 + 367*(month - 2 - (month - 14)/12*12)/12 - 3*((year + 4900 + (month - 14)/12)/100)/4; /* Add the fractional hours, mins and seconds */ date += (greg_time->tm_hour + 12.0)/24.0; date += (greg_time->tm_min)/1440.0; date += (greg_time->tm_sec)/86400.0; return date; } /* * modjuldate * - procedure to calculate the modified Julian date * * Inputs * jd - the Julian date * * Returns * the modified Julian date * */ double modjuldate(double jd) { return jd - 2400000.5; } /* * gethjd * - procedure to calculate the heliocentric modified Julian * date from the MJD value. Requests the co-ordinates of * the source in the procedure. * * Inputs * mjd - the modified Julian date * * Returns * the modified heliocentric Julian date * */ double gethjd(double mjd) { const double ausec = 499.01265; /* Time it takes light to travel 1 AU */ const double degtorad = 0.01745329251; /* Conversion factor from degrees to radians */ struct coordinates { double ra; /* The right ascension in degrees */ double dec; /* The declination in degrees */ double x; /* The X position in AU */ double y; /* The Y position in AU */ double z; /* The Z position in AU */ }; struct coordinates earth; /* Earth co-ordinates */ struct coordinates source; /* Source co-ordinates */ int ra_h; /* Inputted RA hours */ int ra_m; /* Inputted RA minutes */ float ra_s; /* Inputted RA seconds */ int dec_d; /* Inputted Dec degrees */ int dec_m; /* Inputted Dec minutes */ float dec_s; /* Inputted Dec seconds */ double correction_secs; /* Correction factor in seconds from MJD to HMJD */ double cel; /* intermediate calculation from spherical to cartesian co-ordinates */ double hmjd; /* The modified heliocentric Julian date */ /* Defaults */ ra_h = ra_m = 0; ra_s = dec_s = 0.0; dec_d = dec_m = 0; earth.ra = earth.dec = 0.0; /* Enter the co-ordinates of the source */ printf("Enter the co-ordinates of your source in J2000"); printf(" (hh mm ss.s dd mm ss.s)\n"); (void) scanf("%d %d %f %d %d %f", &ra_h, &ra_m, &ra_s, &dec_d, &dec_m, &dec_s); /* Calculate the RA and declination in decimal degree notation */ source.ra = (double)ra_h; source.ra += (ra_m / 60.0); source.ra += (ra_s / 3600.0); source.ra *= 15.0; source.dec = (double)dec_d; source.dec += (dec_m / 60.0); source.dec += (dec_s / 3600.0); /* Attempt to find the RA and Dec of the earth using the astronomical calculator book. */ heliocentric_ra_dec(mjd, &earth.ra, &earth.dec); /* Calculate the heliocentric co-ordinates as X, Y and Z terms */ cel = cos(earth.dec * degtorad); earth.x = cos(earth.ra * degtorad) * cel; earth.y = sin(earth.ra * degtorad) * cel; earth.z = sin(earth.dec * degtorad); /* Calculate the X,Y,Z co-ordinates of your source */ cel = cos(source.dec * degtorad); source.x = cos(source.ra * degtorad) * cel; source.y = sin(source.ra * degtorad) * cel; source.z = sin(source.dec * degtorad); /* Calculate the correction in seconds for the light travel time between the source and the earth vectors in a heliocentric reference frame. */ correction_secs = ausec * (earth.x * source.x + earth.y * source.y + earth.z * source.z); /* Modify the MJD in a heliocentric reference frame */ hmjd = mjd + (correction_secs / (24.0 * 3600)); return hmjd; } /* * heliocentric_ra_dec * - procedure to calculate heliocentric right ascension and * declination of the earth at a given date. * * Inputs * mjd - the modified Julian date * * Returns * *ra - the right ascension of the earth * *dec - the declination of the earth * */ void heliocentric_ra_dec(double mjd, double *ra, double *dec) { const double eccentricity = 0.016718; /* Eccentricity of the Earth's orbit*/ const double ecliptic_long = 278.833540; /* The longitude of the ecliptic at 1 Jan 1980 0:00 UT */ const double perigee_long = 282.596403; /* The longitude of perigee at 1 Jan 1980 00:00 UT */ const double deg_to_rad = M_PI / 180.0; /* A degrees to radians converion */ const double tropical_year = 365.24219572; /* The length of the tropical year in days */ const double obliquity = 23.441884; /* The obliquity of the orbit */ const double mjd_1980 = 44238.0; /* The MJD on 1 Jan 1980 00:00 UT */ double mean_anomoly; /* The mean anomoly of the sun */ double days_from_1980; /* The number of days from 1 Jan 1980 */ double solar_longitude; /* The longitude of the sun */ double number_of_deg; /* The number of degrees in longitude the sun has travelled */ double equation_of_centres; /* The value for the equation of centres */ double x; /* An X position */ double y; /* A Y position */ double beta; /* The ecliptic longitude */ double number_of_rotations; /* An integer number of solar orbits */ /* Calculate the number of days from 1 Jan 1980 */ days_from_1980 = (mjd - mjd_1980); /* Calculate the number of degrees around in the orbit travelled in this time */ number_of_deg = (360.0 / tropical_year) * days_from_1980; /* Adjust so the number of degrees is between 0 and 360 */ if ( (number_of_deg < 0.0) || (number_of_deg > 360.0)) { number_of_rotations = number_of_deg / 360.0; number_of_rotations = floor(number_of_rotations); number_of_deg -= number_of_rotations * 360.0; } /* Calculate the mean anomoly */ mean_anomoly = number_of_deg - perigee_long + ecliptic_long; /* Since the orbit is elliptical and not circular, calculate the equation of centres */ equation_of_centres = (360.0 / M_PI) * eccentricity * sin(mean_anomoly * deg_to_rad); /* Calculate the solar longitude */ solar_longitude = number_of_deg + equation_of_centres + ecliptic_long; if (solar_longitude > 360.0) solar_longitude -= 360.0; /* The ecliptic latitude is zero for the Sun. */ beta = 0.0; /* Calculate the RA and Dec of the sun */ *dec = asin( (sin(beta * deg_to_rad) * cos(obliquity * deg_to_rad)) + (cos(beta * deg_to_rad) * sin(obliquity * deg_to_rad) * sin(solar_longitude * deg_to_rad)) ); *dec /= deg_to_rad; x = cos(solar_longitude * deg_to_rad); y = (sin(solar_longitude * deg_to_rad) * cos(obliquity * deg_to_rad)) - (tan(beta * deg_to_rad) * sin(obliquity * deg_to_rad)); *ra = atan(y/x); *ra /= deg_to_rad; if (*ra < 0.0) *ra += 360.0; *ra /= 15.0; /* Convert from geocentric to heliocentric co-ordinates for the Earth*/ *dec *= -1.0; *ra -= 12.0; if (*ra < 0.0) { *ra += 24.0; } }