// SofaJplDemo.cs
// Demonstration program for SofaJpl DLL v. 2.1.2.
// Paul S. Hirose, 2018 August 26
// DOCUMENTATION
// The Visual Studio Object Browser displays the main SofaJpl documentation. After adding a
// reference to SofaJpl to your project, open Object Browser and expand the HirosePS. SofaJpl
// namespace.
// revision history
#if false
2018-09-07 Added a compatibility mode to duplicate the equatorial, ecliptical, and az/el of JPL
Horizons. The first two agree to 1 mas in my tests. Az/el agrees .0001 deg, which is the Horizons
output precision.
#endif
// Select the input time scale. Exactly one must be #defined, and the others commented out.
#define UTC
// #define UT1
// #define TT
// If #define is un-commented, convert a JPL ASCII ephemeris to binary.
// After the binary is created, this #define should be commented out.
// #define MAKE_EPHEMERIS
// If #define is un-commented, use the delta T you provide. Otherwise use the SofaJpl delta T model.
// Has no effect if time scale is UTC.
// #define MANUAL_DELTA_T
// If un-commented, use a precession nutation model compatible with JPL Horizons: IAU 1976 / 80
// precession / nutation, no frame bias, apply IERS pole offsets.
// NOTE: To obtain coordinates in the horizontal system or the ITRS, compatible with Horizons, it is
// also necessary to activate the MANUAL_DELTA_T option and supply the correct delta T of date.
// In addition, for compatible coordinates in the horizontal system set the deflection of the
// vertical parameters (xi and eta) to zero.
// #define HORIZONS
using System;
// Everything in SofaJpl is in the HirosePS.SofaJpl namespace.
using SofaJpl = HirosePS.SofaJpl;
static class SofaJplDemo {
// The leap second table. One is included in the SofaJpl distribution .zip file.
static string _leapsecTable = @"C:\Users\Stan\Documents\astro\SofaJpl_2_1_2\leapSecTable.txt";
#if MAKE_EPHEMERIS
// JPL ASCII ephemeris and header. Must be downloaded from JPL. Not needed after they are
// converted to a binary ephemeris.
static string _asciiEphemeris = @"C:\Users\Stan\Documents\astro\jpl\ascp2000.422";
static string _ephemerisHeader = @"C:\Users\Stan\Documents\astro\jpl\header.422";
#endif
// The binary JPL ephemeris to create and use.
static string _binaryEphemeris = @"binp2000_2099.422";
// mathematical constants
const double _degPerHour = 15.0;
// an enumeration of the three possible angle output formats
enum angleFormat { D, DM, DMS };
// output format desired by the user
static angleFormat _format;
// These control the output resolution of angle and time. The values are computed automatically
// from the angle accuracy set by the user.
static double _angleResolution; // resolution units per degree
static double _timeResolution; // resolution units per hour
// Number of decimal places in the floating point 'f' format to achieve unit vector rectangular
// coordinate precision comparable to the angles.
static int _vectorResolution;
static void Main(string[] args) {
// Set the display format and angle accuracy.
_format = angleFormat.DMS;
double angleAccuracy = SofaJpl.Angle.DmsToRad(0, 0, 0.1);
// Set the epoch.
int year = 2018;
int month = 6;
int day = 6;
int hour = 4;
int minute = 0;
double second = 0.0;
bool julian = false; // true if date is in Julian calendar
SofaJpl.Duration deltaT;
#if UTC
SofaJpl.Duration ut1MinusUtc = SofaJpl.Duration.FromSeconds(0.07522);
#else
#if MANUAL_DELTA_T
deltaT = SofaJpl.Duration.FromSeconds(69.109);
#else
// Get delta T from the SofaJpl model.
deltaT = SofaJpl.DeltaT.GetDeltaT(
new SofaJpl.JulianDate(year, month, day, hour, minute, second, julian));
#endif
#endif
// Observer position.
string obsName = "Kitt Peak Observatory"; // name of topocenter (optional)
double lon = SofaJpl.Angle.DmsToRad(-111, 36, 00.0); // east longitude
double lat = SofaJpl.Angle.DmsToRad(31, 57, 48.0); // north latitude
float height = 2100f; // meters above ellipsoid
// This instance of the Topocenter class will supply the observer position and velocity
// with respect to the ICRS.
SofaJpl.Topocenter obs = new SofaJpl.Topocenter(lat, lon, height);
// Construct an Atmosphere object to apply refraction. The parameters required by the
// constructor are type "float", thus the f suffix on the numbers.
float altimiterSetting = 29.9f * SofaJpl.Atmosphere.MillibarsPerInchHg; // millibars
float degC = 12f; // Celsius
float dewPointC = 0f; // Celsius.
SofaJpl.Atmosphere atm = new SofaJpl.Atmosphere(height, degC, altimiterSetting,
dewPointC, false);
// polar motion (radians)
double poleX = SofaJpl.Angle.DmsToRad(0, 0, 0.0);
double poleY = SofaJpl.Angle.DmsToRad(0, 0, 0.0);
// deflection of the vertical at the observer (radians)
double xi = SofaJpl.Angle.DmsToRad(0, 0, 0);
double eta = SofaJpl.Angle.DmsToRad(0, 0, 0);
// Create a binary ephemeris from JPL ASCII files.
#if MAKE_EPHEMERIS
Console.WriteLine("\nCREATING BINARY JPL EPHEMERIS.\n");
SofaJpl.JplEphemeris.AsciiToBinary(_ephemerisHeader, _asciiEphemeris, _binaryEphemeris);
#endif
// Open the binary JPL ephemeris. This loads the entire ephemeris into RAM. Even if the
// target body is a star, a solar system ephemeris is required for parallax, aberration,
// and light deflection due to solar gravitation.
SofaJpl.JplEphemeris eph = new SofaJpl.JplEphemeris(_binaryEphemeris);
// There are several ways to specify the target body.
// Use the SofaJpl star catalog, a subset of the Hipparcos catalog complete to mag. 3.
#if false
// Create a HipparcosCatalog object.
SofaJpl.HipparcosCatalog catalog = new SofaJpl.HipparcosCatalog();
// Look up the star name in the catalog we just constructed, and construct a Star object.
// The current SofaJpl implementation requires a full match (not substring) to a designation
// in the name dictionary. It may be necessary to examine the dictionary (star_names.vot)
// with a text editor to find the correct name of a star.
SofaJpl.Body body = catalog.GetStar("name Vega", eph);
#endif
// Star, from data you supply. The first parameter is the epoch of the data, J2000.0 in
// this case.
#if false
SofaJpl.Body body = new SofaJpl.Star(SofaJpl.JulianDate.J2000Base, 293.0899579 / 15.0,
+69.6611767, 173.77, 597.482, -1738.313, 26.78, eph, "sig Dra");
#endif
// Solar system body in the JPL ephemeris. The correction for light time iterates until the
// solution converges to the specified angle accuracy.
#if true
SofaJpl.Body body = new SofaJpl.JplBody(SofaJpl.JplEphemeris.Body.Venus, eph,
angleAccuracy);
#endif
// Finished setting parameters. Calculate and display results.
// Compute the angle, time, and rectangular coordinate resolution consistent with the angle
// accuracy (radians) specified by user. These variables control output formatting.
_angleResolution = 1.0 / SofaJpl.Angle.RadiansToDegrees(angleAccuracy);
_timeResolution = _angleResolution * _degPerHour;
_vectorResolution = (int)(-Math.Log10(angleAccuracy) + 1.0);
if (_vectorResolution < 1)
_vectorResolution = 1; // Always display at least 1 decimal place.
// Load the leap second table.
SofaJpl.Utc.LoadTableFromFile(_leapsecTable);
// Calculate UT1 and TT. UTC is not calculated (yet) unless it's the input time scale.
SofaJpl.JulianDate ut1, tt;
SofaJpl.Utc utc;
#if UTC
utc = new SofaJpl.Utc(year, month, day, hour, minute, second, julian);
tt = utc.TerrestrialTime;
ut1 = new SofaJpl.JulianDate(year, month, day, hour, minute, second, julian) + ut1MinusUtc;
deltaT = tt - ut1;
#endif
#if UT1
ut1 = new SofaJpl.JulianDate(year, month, day, hour, minute, second, julian);
tt = ut1 + deltaT;
#endif
#if TT
tt = new SofaJpl.JulianDate(year, month, day, hour, minute, second, julian);
ut1 = tt - deltaT;
#endif
// Display date and time in UTC.
// The TimeFields class breaks down a JulianDate into year, month, etc.
SofaJpl.TimeFields tf1;
#if UTC
tf1 = utc.ToTimeFields(_timeResolution, julian);
Console.WriteLine("{0} UTC", tf1);
#else
// The time input from the user was UT1 or TT, so it's not safe to assume conversion to UTC
// is possible. It must be compared to the UTC table boundaries, which are in terms of TAI.
SofaJpl.JulianDate tai = tt - SofaJpl.Duration.TTMinusTai;
if (tai >= SofaJpl.Utc.DefaultTable.FirstTai && tai <= SofaJpl.Utc.DefaultTable.LastTai) {
utc = new SofaJpl.Utc(tt);
tf1 = utc.ToTimeFields(_timeResolution, julian);
Console.WriteLine("{0} UTC", tf1);
}
#endif
// Display UT1.
tf1 = ut1.ToTimeFields(_timeResolution, julian);
Console.WriteLine("{0} UT1", tf1);
// Display TT as date and time, and also as 2-part Julian date.
tf1 = tt.ToTimeFields(_timeResolution, julian);
Console.WriteLine("{0} TT", tf1);
Console.WriteLine("Dates are {0} calendar.", julian ? "Julian" : "Gregorian");
Console.WriteLine("JD {0} TT", tt);
// Display delta T.
SofaJpl.Sexagesimal sex1 = new SofaJpl.Sexagesimal(deltaT.ToHours(), _timeResolution);
Console.WriteLine("{0:+fhms} delta T", sex1);
// Display polar motion.
Console.WriteLine("\npolar motion");
displayPolarMotion(poleX, poleY, "x, y");
// Display topocenter geodetic and rectangular coordinates.
Console.WriteLine();
Console.WriteLine(obsName); // name of observatory
displayGeodetic(lon, lat, "E lon, N lat");
Console.WriteLine("{0:f1} meters above ellipsoid", height);
SofaJpl.Vector obsVec = obs.ToVector();
displayXyz(obsVec, "ITRS unit vector", "modulus (km)");
// Display deflection of the vertical. The same format as polar motion is appropriate.
Console.WriteLine("\ndeflection of the vertical");
displayPolarMotion(xi, eta, "xi, eta");
// Display atmosphere conditions.
Console.WriteLine("\n{0:f0} C ({1:f0} F) at observer",
atm.StationTemperature,
atm.StationTemperature * 1.8 + 32.0);
Console.WriteLine("{0,6:f1} mb ({1,5:f2}\" Hg) altimeter setting",
atm.AltimeterSetting,
atm.AltimeterSetting / SofaJpl.Atmosphere.MillibarsPerInchHg);
Console.WriteLine("{0,6:f1} mb ({1,5:f2}\" Hg) station pressure",
atm.StationPressure,
atm.StationPressure / SofaJpl.Atmosphere.MillibarsPerInchHg);
if (atm.HumidityIsRelative)
Console.WriteLine("{0:f0}% relative humidity", atm.Humidity);
else
Console.WriteLine("{0:f0} C ({1:f0} F) dew point",
atm.Humidity,
atm.Humidity * 1.8 + 32.0);
// Create rotation matrices for the coordinate transformations.
// ITRS to horizontal (east, north, zenith) system, including deflection of the vertical
SofaJpl.RMatrix itrstoHor = SofaJpl.RMatrix.ItrsToHor(lat, lon, xi, eta);
// polar motion matrix (terrestrial intermediate system to ITRS)
SofaJpl.RMatrix tirsToItrs = SofaJpl.RMatrix.TirsToItrs(tt, poleX, poleY);
// GCRS to true equator & equinox (IAU 2006 precession and 2000B nutation).
SofaJpl.RMatrix gcrsToMean06 = SofaJpl.RMatrix.Precess06(tt);
double eps06 = SofaJpl.RMatrix.MeanObliq06(tt); // mean obliquity
double dPsi00, dEps00; // nutation in longitude and obliquity
SofaJpl.RMatrix.NutationAngles00b(tt, out dPsi00, out dEps00);
SofaJpl.RMatrix gcrsToTrue06 = SofaJpl.RMatrix.Nutate(eps06, dPsi00, dEps00) *
gcrsToMean06;
#if HORIZONS // Generate coordinates compatible with JPL Horizons (1976/80 precession/nutation).
// In SofaJpl, 1976 precession includes frame bias, which must be removed for
// Horizons compatibility.
SofaJpl.RMatrix gcrsToMean = SofaJpl.RMatrix.Precess76(tt) *
SofaJpl.RMatrix.IcrsToJ2000.Transpose();
// Calculate the GCRS to ecliptic (mean equinox) rotation matrix.
double eps = SofaJpl.RMatrix.MeanObliq80(tt); // mean obl.
SofaJpl.RMatrix gcrsToEclipticMean = SofaJpl.RMatrix.GcrsToEclip(gcrsToMean, eps);
// GCRS vector to the CIP, based on the IAU 2006 precession & 2000 nutation models.
SofaJpl.Vector cipGcrs06 = gcrsToTrue06.Row(3);
// Transform it to spherical coords in the 1976/80 ecliptic and mean equinox system.
SofaJpl.Vector cipEclip06 = gcrsToEclipticMean * cipGcrs06;
SofaJpl.Spherical cipSph = new SofaJpl.Spherical(cipEclip06);
// Derive and apply the nutation angles to obtain the GCRS to true equator/equinox matrix.
double dPsi, dEps;
dPsi = SofaJpl.Angle.HalfPi - cipSph.LonEast;
dEps = cipSph.NPD - eps;
SofaJpl.RMatrix gcrsToTrue = SofaJpl.RMatrix.Nutate(eps, dPsi, dEps) * gcrsToMean;
// Compute the GCRS to terrestrial intermediate matrix
double gast = SofaJpl.Angle.Gast94(ut1);
SofaJpl.RMatrix gcrsToTirs = SofaJpl.RMatrix.GcrsToTirs(gcrsToTrue, gast);
#else // Generate coordinates compatible with IAU 2006/00 precession/nutation.
// Mean obliquity and nutation in obliquity have already been computed.
double eps = eps06;
double dEps = dEps00;
// GCRS to true equator/equinox matrix has already been computed.
SofaJpl.RMatrix gcrsToTrue = gcrsToTrue06;
// Get X and Y of the celestial intermediate pole.
double cipX, cipY;
gcrsToTrue.CipXY(out cipX, out cipY);
// Form the GCRS to celestial intermediate matrix.
SofaJpl.RMatrix gcrsToCirs = SofaJpl.RMatrix.GcrsToCirs(cipX, cipY,
SofaJpl.RMatrix.S06(tt, cipX, cipY));
// Greenwich apparent sidereal time.
double gast = SofaJpl.Angle.Gast06b(ut1, tt);
// Earth rotation angle
double era = SofaJpl.Angle.Era00(ut1);
// Compute the GCRS to terrestrial intermediate matrix
SofaJpl.RMatrix gcrsToTirs = SofaJpl.RMatrix.GcrsToTirs(gcrsToCirs, era);
#endif
// End code compatible with IAU 2006/00 precession/nutation. The remaining rotation
// matrix computations are common to IAU 2006/00 and JPL Horizons modes.
// GCRS to ecliptic and true equinox
SofaJpl.RMatrix gcrsToEcliptic = SofaJpl.RMatrix.GcrsToEclip(gcrsToTrue, eps + dEps);
// GCRS to ITRS, including polar motion
SofaJpl.RMatrix gcrsToItrs = tirsToItrs * gcrsToTirs;
// GCRS to horizontal, including deflection of the vertical.
SofaJpl.RMatrix gcrsToHor = itrstoHor * gcrsToItrs;
// Display barycentric coordinates of the body.
Console.WriteLine("\n{0} barycentric position & velocity", body.Name);
SofaJpl.PVVector pv1 = body.Barycentric(tt);
displayRaDec(pv1.Position, "RA, dec (ICRS)");
displayXyz(pv1.Position, "unit vector (ICRS)", "distance (km)");
displayXyz(pv1.Velocity, "velocity unit vector", "km/day");
// heliocentric coordinates
Console.WriteLine("\n{0} heliocentric position & velocity", body.Name);
pv1 = body.Heliocentric(tt);
displayRaDec(pv1.Position, "RA, dec (ICRS)");
displayXyz(pv1.Position, "unit vector (ICRS)", "distance (km)");
displayXyz(pv1.Velocity, "velocity unit vector", "km/day");
// geocentric coordinates
Console.WriteLine("\n{0} geocentric geometric position & velocity", body.Name);
pv1 = body.GeocentricGeometric(tt);
displayRaDec(pv1.Position, "RA, dec (ICRS)");
displayXyz(pv1.Position, "unit vector (ICRS)", "geometric distance (km)");
displayXyz(pv1.Velocity, "velocity unit vector", "km/day");
Console.WriteLine("\n{0} geocentric astrometric place", body.Name);
SofaJpl.Vector vec1 = body.GeocentricAstrometric(tt);
displayRaDec(vec1, "RA, dec (ICRS)");
displayModulus(vec1, "astrometric distance (km)");
Console.WriteLine("\n{0} geocentric apparent place", body.Name);
vec1 = body.GeocentricApparent(tt);
displayRaDec(vec1, "RA, dec (ICRS)");
displayRaDec(gcrsToTrue * vec1, "equinox RA, dec");
#if ! HORIZONS
displayRaDec(gcrsToCirs * vec1, "intermediate RA, dec");
#endif
displayEcliptical(gcrsToEcliptic * vec1, "ecliptic true lon, lat");
// Geographic position of the body in the ITRS. If the user provided the parameters, this
// includes polar motion.
SofaJpl.Spherical sph1 = new SofaJpl.Spherical(gcrsToItrs * vec1);
displayGeodetic(sph1.LonEast, sph1.Lat, "E lon, N lat (ITRS)");
displayGeodetic(sph1.LonWest, sph1.Lat, "W lon, N lat (ITRS)");
displayHourAngle(SofaJpl.Angle.NormPlus(sph1.LonWest + lon), "LHA");
// Semidiameter. All bodies have a property to give radius in km. Stars default to zero.
// Solar system bodies are automatically initialized to their adopted IAU radii. This
// property can be modified by the user.
double sd = SofaJpl.Angle.Semidiameter(vec1.Modulus(), body.Radius);
displaySemidiameter(sd, "geocentric semidiameter");
// Greenwich apparent sidereal time
Console.WriteLine();
displaySiderealTime(gast, "Greenwich apparent sidereal time");
#if ! HORIZONS
// Earth rotation angle.
displayHourAngle(era, "Earth rotation angle");
#endif
// topocentric coordinates
// Get the GCRS position and velocity of the topocenter in a PVVector (position and velocity
// vector). The parameter passed to ToGcrsPV() must be the ITRS to GCRS rotation matrix, but
// what we calculated earlier does the opposite transformation. Thus it's transposed to
// reverse its sense.
SofaJpl.PVVector obsPV = obs.ToGcrsPV(gcrsToItrs.Transpose());
Console.WriteLine("\n{0} topocentric geometric place", body.Name);
vec1 = body.TopocentricGeometric(tt, obsPV).Position;
displayRaDec(vec1, "RA, dec (ICRS)");
displayModulus(vec1, "km");
Console.WriteLine("\n{0} topocentric astrometric place", body.Name);
vec1 = body.TopocentricAstrometric(tt, obsPV.Position);
displayRaDec(vec1, "RA, dec (ICRS)");
displayModulus(vec1, "km");
Console.WriteLine("\n{0} topocentric apparent place", body.Name);
vec1 = body.TopocentricApparent(tt, obsPV);
displayRaDec(vec1, "RA, dec (ICRS)");
// Transform apparent place from the GCRS to the horizontal system. If the user supplied
// the parameters, this includes polar motion and deflection of the vertical.
vec1 = gcrsToHor * vec1;
// Convert to spherical coordinates.
SofaJpl.Spherical sphUnref = new SofaJpl.Spherical(vec1);
// Use the Atmosphere object created earlier to generate refracted spherical coordinates.
// The application of refraction is iterative, so the desired accuracy must be passed
// to Refract().
SofaJpl.Spherical sphRefr = atm.Refract(sphUnref, angleAccuracy);
Console.WriteLine("\n{0} azimuth, zenith distance, elevation", body.Name);
displayAzZd(sphUnref, "az, unrefracted ZD");
displayAzEl(sphUnref, "az, unrefracted el");
displayAzEl(sphRefr, "az, refracted el");
// Semidiameter. All bodies have a property to give radius in km. Stars default to zero, but
// solar system bodies are automatically initialized to their adopted IAU radii.
sd = SofaJpl.Angle.Semidiameter(vec1.Modulus(), body.Radius);
displaySemidiameter(sd, "topocentric semidiameter");
// If applicable, display phase angle: the separation angle, at the body, between
// vectors directed to the Sun and the observer. These are the negatives of vectors to the
// body's heliocentric and topocentric positions in the ICRS. If we omit the negations,
// both vectors are off by 180 degrees and thus the angle between them is still correct.
if (body.IsSolarSystemBody && body.IsSun == false) {
SofaJpl.JplBody earth = new SofaJpl.JplBody(
SofaJpl.JplEphemeris.Body.Earth, eph, angleAccuracy);
SofaJpl.JplBody sun = new SofaJpl.JplBody(
SofaJpl.JplEphemeris.Body.Sun, eph, angleAccuracy);
SofaJpl.Vector earthVec = body.TopocentricGeometric(tt, obsPV).Position;
SofaJpl.Vector sunVec = body.Heliocentric(tt).Position;
double pa = earthVec.SeparationAngle(sunVec);
Console.WriteLine("{0:f0}° phase angle (0 = full, 180 = new)",
SofaJpl.Angle.RadiansToDegrees(pa));
}
} // end Main()
// helper methods to display data with the format and precision selected by user
/// <summary>
/// Display the modulus of a vector.
/// </summary>
/// <param name="vec">the vector</param>
/// <param name="label">string to display after the modulus</param>
static void displayModulus(SofaJpl.Vector vec, string label) {
System.Text.StringBuilder sb1 = new System.Text.StringBuilder("{0:e", 8);
// The _vectorResolution value is appropriate for the 'f' format. But the 'e' format gives
// one more significant digit, so subtract 1.
sb1.Append(_vectorResolution - 1);
sb1.Append("} ");
Console.WriteLine(sb1.ToString() + label, vec.Modulus());
}
/// <summary>
/// Display hour angle as a sexagesimal.
/// </summary>
/// <param name="angle">hour angle (radians) </param>
/// <param name="str">string to display after the angle</param>
static void displayHourAngle(double angle, string str) {
string formatString;
switch (_format) {
case angleFormat.D:
formatString = "{0:3c°'\"} {1}";
break;
case angleFormat.DM:
formatString = "{0:3b°'\"} {1}";
break;
default:
formatString = "{0:3a°'\"} {1}";
break;
}
SofaJpl.Sexagesimal sex = new SofaJpl.Sexagesimal(
SofaJpl.Angle.RadiansToDegrees(angle), _angleResolution, 360);
Console.WriteLine(formatString, sex, str);
}
/// <summary>
/// Display sidereal time as a sexagesimal.
/// </summary>
/// <param name="gast">sidereal time (radians)</param>
/// <param name="label">string to display after the angle</param>
/// <remarks>Display degrees (not hours) unless format is DMS.</remarks>
static void displaySiderealTime(double gast, string label) {
string formatString;
switch (_format) {
case angleFormat.D:
formatString = "{0:3c°'\"} {1}";
break;
case angleFormat.DM:
formatString = "{0:3b°'\"} {1}";
break;
default:
formatString = "{0:2ahms} {1}";
break;
}
SofaJpl.Sexagesimal sex;
if (_format == angleFormat.DMS)
sex = new SofaJpl.Sexagesimal(SofaJpl.Angle.RadiansToHours(gast),
_timeResolution, 24);
else
sex = new SofaJpl.Sexagesimal(SofaJpl.Angle.RadiansToDegrees(gast),
_angleResolution, 360);
Console.WriteLine(formatString, sex, label);
}
/// <summary>
/// Display semidiameter as a sexagesimal.
/// </summary>
/// <param name="sd">semidiameter (radians)</param>
/// <param name="label">string to display after the angle</param>
static void displaySemidiameter(double sd, string label) {
string formatString;
switch (_format) {
case angleFormat.D:
formatString = "{0:c°'\"} {1}";
break;
case angleFormat.DM:
formatString = "{0:e°'\"} {1}";
break;
default:
formatString = "{0:f°'\"} {1}";
break;
}
sd = SofaJpl.Angle.RadiansToDegrees(sd);
SofaJpl.Sexagesimal sex = new SofaJpl.Sexagesimal(sd, _angleResolution);
Console.WriteLine(formatString, sex, label);
}
/// <summary>
/// Display polar motion angles as sexagesimals.
/// </summary>
/// <param name="poleX">pole X (radians)</param>
/// <param name="poleY">pole Y (radians)</param>
/// <param name="label">string to display after the angles</param>
static void displayPolarMotion(double poleX, double poleY, string label) {
string formatString;
switch (_format) {
case angleFormat.D:
formatString = "{0:+c°'\"} {1:+c°'\"} {2}";
break;
case angleFormat.DM:
formatString = "{0:+e°'\"} {1:+e°'\"} {2}";
break;
default:
formatString = "{0:+f°'\"} {1:+f°'\"} {2}";
break;
}
double xDeg = SofaJpl.Angle.RadiansToDegrees(poleX);
double yDeg = SofaJpl.Angle.RadiansToDegrees(poleY);
SofaJpl.Sexagesimal sexX = new SofaJpl.Sexagesimal(xDeg, _angleResolution);
SofaJpl.Sexagesimal sexY = new SofaJpl.Sexagesimal(yDeg, _angleResolution);
Console.WriteLine(formatString, sexX, sexY, label);
}
/// <summary>
/// Display a vector as right ascension and declination.
/// </summary>
/// <param name="vec">vector to the body with respect to the equatorial system</param>
/// <param name="label">string to display after the angles</param>
/// <remarks>Unless the format is DMS, display RA in degrees not hours.</remarks>
static void displayRaDec(SofaJpl.Vector vec, string label) {
string formatString;
switch (_format) {
case angleFormat.D:
formatString = "{0:3c°'\"} {1:+2c°'\"} {2}";
break;
case angleFormat.DM:
formatString = "{0:3b°'\"} {1:+2b°'\"} {2}";
break;
default:
formatString = "{0:2ahms} {1:+2a°'\"} {2}";
break;
}
SofaJpl.Spherical sph = new SofaJpl.Spherical(vec);
SofaJpl.Sexagesimal raSex;
if (_format == angleFormat.DMS)
raSex = new SofaJpl.Sexagesimal(SofaJpl.Angle.RadiansToHours(sph.RA),
_timeResolution, 24);
else
raSex = new SofaJpl.Sexagesimal(SofaJpl.Angle.RadiansToDegrees(sph.RA),
_angleResolution, 360);
SofaJpl.Sexagesimal decSex = new SofaJpl.Sexagesimal(
SofaJpl.Angle.RadiansToDegrees(sph.Dec), _angleResolution);
Console.WriteLine(formatString, raSex, decSex, label);
}
/// <summary>
/// Display a vector as azimuth and elevation.
/// </summary>
/// <param name="vec">vector to the body with respect to the horizontal system</param>
/// <param name="label">string to display after the angles</param>
static void displayAzEl(SofaJpl.Vector vec, string label) {
string formatString;
switch (_format) {
case angleFormat.D:
formatString = "{0:3c°'\"} {1:+2c°'\"} {2}";
break;
case angleFormat.DM:
formatString = "{0:3b°'\"} {1:+2b°'\"} {2}";
break;
default:
formatString = "{0:3a°'\"} {1:+2a°'\"} {2}";
break;
}
SofaJpl.Spherical sph = new SofaJpl.Spherical(vec);
SofaJpl.Sexagesimal azSex = new SofaJpl.Sexagesimal(
SofaJpl.Angle.RadiansToDegrees(sph.Az), _angleResolution, 360);
SofaJpl.Sexagesimal elSex = new SofaJpl.Sexagesimal(
SofaJpl.Angle.RadiansToDegrees(sph.El), _angleResolution);
Console.WriteLine(formatString, azSex, elSex, label);
}
/// <summary>
/// Display a vector as ecliptic longitude and latitude.
/// </summary>
/// <param name="vec">vector to the body with respect to the ecliptic system</param>
/// <param name="label">string to display after the angles</param>
static void displayEcliptical(SofaJpl.Vector vec, string label) {
string formatString;
switch (_format) {
case angleFormat.D:
formatString = "{0:3c°'\"} {1:+2c°'\"} {2}";
break;
case angleFormat.DM:
formatString = "{0:3b°'\"} {1:+2b°'\"} {2}";
break;
default:
formatString = "{0:3a°'\"} {1:+2a°'\"} {2}";
break;
}
SofaJpl.Spherical sph = new SofaJpl.Spherical(vec);
SofaJpl.Sexagesimal lonSex = new SofaJpl.Sexagesimal(
SofaJpl.Angle.RadiansToDegrees(sph.LonEast), _angleResolution, 360);
SofaJpl.Sexagesimal latSex = new SofaJpl.Sexagesimal(
SofaJpl.Angle.RadiansToDegrees(sph.Lat), _angleResolution);
Console.WriteLine(formatString, lonSex, latSex, label);
}
/// <summary>
/// Display geodetic longitude and latitude.
/// </summary>
/// <param name="lon">east longitude (radians)</param>
/// <param name="lat">north latitude (radians)</param>
/// <param name="label">string to display after the angles</param>
static void displayGeodetic(double lon, double lat, string label) {
string formatString;
switch (_format) {
case angleFormat.D:
formatString = "{0:3c°'\"} {1:+2c°'\"} {2}";
break;
case angleFormat.DM:
formatString = "{0:3b°'\"} {1:+2b°'\"} {2}";
break;
default:
formatString = "{0:3a°'\"} {1:+2a°'\"} {2}";
break;
}
SofaJpl.Sexagesimal lonSex = new SofaJpl.Sexagesimal(
SofaJpl.Angle.RadiansToDegrees(lon), _angleResolution, 360);
SofaJpl.Sexagesimal latSex = new SofaJpl.Sexagesimal(
SofaJpl.Angle.RadiansToDegrees(lat), _angleResolution);
Console.WriteLine(formatString, lonSex, latSex, label);
}
/// <summary>
/// Display a vector as xyz components of a unit vector and modulus.
/// </summary>
/// <param name="v">vector</param>
/// <param name="xyzLabel">label for the unit vector</param>
/// <param name="modulusLabel">label for the modulus</param>
static void displayXyz(SofaJpl.Vector v, string xyzLabel, string modulusLabel) {
System.Text.StringBuilder sb1 = new System.Text.StringBuilder("{0:f", 23);
sb1.Append(_vectorResolution);
sb1.Append("} {1:f");
sb1.Append(_vectorResolution);
sb1.Append("} {2:f");
sb1.Append(_vectorResolution);
sb1.Append("} ");
SofaJpl.Vector uv = v.Unit();
Console.WriteLine(sb1.ToString() + xyzLabel, uv.X, uv.Y, uv.Z);
displayModulus(v, modulusLabel);
}
/// <summary>
/// Given a Spherical object, display azimuth and elevation as sexagesimals.
/// </summary>
/// <param name="sph">spherical coordinates in the horizontal system</param>
/// <param name="label">string to display after the angles</param>
static void displayAzEl(SofaJpl.Spherical sph, string label) {
string formatString;
switch (_format) {
case angleFormat.D:
formatString = "{0:3c°'\"} {1: 2c°'\"} {2}";
break;
case angleFormat.DM:
formatString = "{0:3b°'\"} {1: 2b°'\"} {2}";
break;
default:
formatString = "{0:3a°'\"} {1: 2a°'\"} {2}";
break;
}
double az = SofaJpl.Angle.RadiansToDegrees(sph.Az);
double el = SofaJpl.Angle.RadiansToDegrees(sph.El);
SofaJpl.Sexagesimal azSex = new SofaJpl.Sexagesimal(az, _angleResolution, 360);
SofaJpl.Sexagesimal elSex = new SofaJpl.Sexagesimal(el, _angleResolution);
Console.WriteLine(formatString, azSex, elSex, label);
}
/// <summary>
/// Given a Spherical object, display azimuth and zenith distance as sexagesimals.
/// </summary>
/// <param name="sph">spherical coordinates</param>
/// <param name="label">string to display after the angles</param>
static void displayAzZd(SofaJpl.Spherical sph, string label) {
string formatString;
switch (_format) {
case angleFormat.D:
// Zenith distance needs a different format than elevation. There's never a negative
// sign, and there can be up to 3 digits before the decimal point.
formatString = "{0:3c°'\"} {1:3c°'\"} {2}";
break;
case angleFormat.DM:
formatString = "{0:3b°'\"} {1:3b°'\"} {2}";
break;
default:
formatString = "{0:3a°'\"} {1:3a°'\"} {2}";
break;
}
double az = SofaJpl.Angle.RadiansToDegrees(sph.Az);
double zd = SofaJpl.Angle.RadiansToDegrees(sph.ZenithDistance);
SofaJpl.Sexagesimal azSex = new SofaJpl.Sexagesimal(az, _angleResolution, 360);
SofaJpl.Sexagesimal zdSex = new SofaJpl.Sexagesimal(zd, _angleResolution);
Console.WriteLine(formatString, azSex, zdSex, label);
}
} // end class
[back]
(last modified 2018-09-09)