Skip to content
This repository has been archived by the owner on May 29, 2022. It is now read-only.

Commit

Permalink
Reapplying JD_MJD changes submitted by Brian Tolman, along with some …
Browse files Browse the repository at this point in the history 
…UpdateStats tests. Changes were originally reverted after merging due to SGLTk issues. When re-applying changes, git merge-base encountered issues recognizing all changed files. Creating a patch from full diff was the best solution to resolving this issue.
  • Loading branch information
Bryan Parsons committed Jul 12, 2017
1 parent 53b809e commit beff4fb
Show file tree
Hide file tree
Showing 32 changed files with 972 additions and 491 deletions.
1 change: 0 additions & 1 deletion core/apps/time/timeconvert.cpp
View file
Expand Up @@ -199,7 +199,6 @@ void TimCvt::process()
string eight(8, ' '); // eight spaces

GPSWeekZcount wz(ct);
CivilTime civ(ct);

cout << endl
<< eight << leftJustify("Month/Day/Year H:M:S", 32)
Expand Down
2 changes: 1 addition & 1 deletion core/lib/GNSSCore/GlobalTropModel.cpp
View file
Expand Up @@ -601,7 +601,7 @@ namespace gpstk
// @param time CommonTime of interest
void GlobalTropModel::setTime(const CommonTime& time)
{
double mjd = static_cast<MJD>(time).mjd;
double mjd = static_cast<MJD>(time).asLongDouble();
setTime(mjd);
}

Expand Down
4 changes: 2 additions & 2 deletions core/lib/TimeHandling/Epoch.hpp
View file
Expand Up @@ -743,7 +743,7 @@ namespace gpstk
long double Epoch::JD() const
throw(Epoch::EpochException)
{
return get<JulianDate>().jd;
return get<JulianDate>().JD();
}

/// Get Modified Julian Date MJD
Expand All @@ -752,7 +752,7 @@ namespace gpstk
long double Epoch::MJD() const
throw(Epoch::EpochException)
{
return get<gpstk::MJD>().mjd; // gpstk to distinguish from Epoch::MJD
return get<gpstk::MJD>().asLongDouble();
}

/// Get year.
Expand Down
241 changes: 210 additions & 31 deletions core/lib/TimeHandling/JulianDate.cpp
View file
Expand Up @@ -39,12 +39,172 @@
#include <cmath>
#include "JulianDate.hpp"
#include "TimeConstants.hpp"
#include "StringUtils.hpp"

namespace gpstk
{
const unsigned int JulianDate::JDLEN(17); // < # dec digits uint64_t
const double JulianDate::JDFACT(1.0e-17); // 1.0e-JDLEN
//const uint64_t JulianDate::JDHALFDAY(500000000L); // JDLEN=9 digits
//const uint64_t JulianDate::JDHALFDAY(500000000000L); // JDLEN=12 digits
const uint64_t JulianDate::JDHALFDAY(50000000000000000L);// JDLEN=17 digits

// Constructor from long double JD
// Warning - precision lost on systems where long double == double (WIN)
JulianDate::JulianDate(long double jd)
{
if(jd < 0.0L)
GPSTK_THROW(InvalidParameter("Invalid input"));
jday = static_cast<long>(jd+0.5);
jd -= static_cast<long double>(jday);
if(jd >= 0.5) { dday = static_cast<uint64_t>((jd-0.5L)/JDFACT); jday += 1L; }
else dday = static_cast<uint64_t>((jd+0.5L)/JDFACT);
fday = static_cast<uint64_t>((jd/JDFACT-dday)/JDFACT);
}

// Constructor from Julian day (not JD) and seconds-of-day
JulianDate::JulianDate(long jd, int isod, double fsod, TimeSystem ts)
{
if(jd < 0)
GPSTK_THROW(InvalidParameter("Invalid jday input"));
if(isod < 0 || isod >= SEC_PER_DAY)
GPSTK_THROW(InvalidParameter("Invalid sec-of-day input"));
if(fsod < 0.0 || fsod >= 1.0)
GPSTK_THROW(InvalidParameter("Invalid frac-sec-of-day input"));

jday = jd;
dday = static_cast<uint64_t>(0);
fday = static_cast<uint64_t>(0);

double fracday;
if(isod > 0) {
fracday = static_cast<double>(isod)/SEC_PER_DAY; // NB SEC_PER_DAY=int
dday = static_cast<uint64_t>(fracday/JDFACT);
fracday -= static_cast<double>(dday)*JDFACT;
fday = static_cast<uint64_t>(fracday/(JDFACT*JDFACT));
}
if(fsod > 0.0) {
fracday = fsod/SEC_PER_DAY;
uint64_t fdday = static_cast<uint64_t>(fracday/JDFACT);
fracday -= static_cast<double>(fdday)*JDFACT;
uint64_t ffday = static_cast<uint64_t>(fracday/(JDFACT*JDFACT));

if(ffday > 0) { fday += ffday; }
if(fday >= 2*JDHALFDAY) { fday -= 2*JDHALFDAY; dday += 1L; }
if(fdday > 0) { dday += fdday; }
if(dday >= 2*JDHALFDAY) { dday -= 2*JDHALFDAY; jday += 1L; }
}
timeSystem = ts;
}

// Constructor from long int(JD) and frac(JD)
void JulianDate::fromIntFrac(long ijd, double fjd, TimeSystem ts)
{
if(ijd < 0 || fjd < 0.0 || fjd >= 1.0)
GPSTK_THROW(InvalidParameter("Invalid input"));

bool rnd(fjd >= 0.5);
jday = ijd + (rnd ? 1L : 0L);
if(rnd) dday = static_cast<uint64_t>((fjd-0.5)/JDFACT);
else dday = static_cast<uint64_t>((fjd+0.5)/JDFACT);
fday = static_cast<uint64_t>((fjd - dday*JDFACT)/(JDFACT*JDFACT));
timeSystem = ts;
}

// Constructor (except for system ) from string
void JulianDate::fromString(std::string instr)
{
// parse the string
int sign,exp,index;
std::string str = StringUtils::parseScientific(instr, sign, exp, index);

// cannot have negative
if(sign < 0) GPSTK_THROW(Exception("Negative JD"));

// mod the string to make exp=0
if(exp != 0) { index += exp; exp = 0; }
// pad the string to put index within the string
if(index < 0) // left leading zeros
str = std::string(-index,'0') + str;
else if(index >= str.size()) // trailing zeros
str = str + std::string(str.size()-index,'0');

// break into 3 strings int (.) 17-dig 17-dig
std::string istr("0");
if(index > 0) {
istr = str.substr(0,index);
StringUtils::stripLeading(str,istr);
}
// 64 bit long max value is 9223372036854775807, 19 digits
std::string dstr = (str.length() > 0 ? str.substr(0,JDLEN) : "0");
StringUtils::leftJustify(dstr,JDLEN,'0');
// truncate string after 17 digits, 17+17=34 digits past index
std::string fstr = (str.length() > JDLEN ? str.substr(JDLEN,JDLEN) : "0");
StringUtils::leftJustify(fstr,JDLEN,'0');

bool rnd(dstr[0] >= '5');
jday = std::strtol(istr.c_str(),0,10) + (rnd ? 1 : 0);
dday = strtoull(dstr.c_str(),0,10);
if(rnd) dday -= JDHALFDAY;
else dday += JDHALFDAY; // this accnts for 0.5d JD-jday
fday = strtoull(fstr.c_str(),0,10);
}

// write full JD to string.
std::string JulianDate::asString(const int prec) const
{
long j(jday);
uint64_t d(dday),f(fday);
static const char ten('9'+1);

if(dday < JDHALFDAY) { d += JDHALFDAY; j -= 1L; }
else { d -= JDHALFDAY; }

// write the part to right of '.' first
std::string str;
std::ostringstream oss;
oss << std::setfill('0') << std::setw(JDLEN) << d;
oss << std::setw(JDLEN) << f;
oss << std::setfill(' ');
str = oss.str();

// one very special case : JD == 0.0 to all precision
if(j == -1L && prec == -1) { j=0L; str = std::string(2*JDLEN,'0'); }

// trim to prec digits
if(prec > -1) {
if(prec < str.length()) {
if(str[prec] >= '5' && str[prec] <= '9') {
// round the string at prec
int k(prec-1);
bool rnd(str[k]=='9');
str[k] = (rnd ? '0' : str[k]+1);

while(rnd && --k >= 0) {
str[k] += 1;
rnd = (str[k]==ten);
if(rnd) str[k] = '0';
}
if(rnd) j += 1L;
}
// truncate at prec-1
str = str.substr(0,prec);
}
else
// 0-pad on the right
str = str + std::string(prec-str.length(),'0');
}

oss.str("");
oss << j << "." << str;
return oss.str();
}

JulianDate& JulianDate::operator=( const JulianDate& right )
{
jd = right.jd;
jday = right.jday;
dday = right.dday;
fday = right.fday;
timeSystem = right.timeSystem;
return *this;
}
Expand All @@ -53,16 +213,16 @@ namespace gpstk
{
try
{
long double temp_jd( jd + 0.5 );
long jday( static_cast<long>( temp_jd ) );
long double sod =
( temp_jd - static_cast<long double>( jday ) ) * SEC_PER_DAY;

// fraction of day
double frod = static_cast<double>(dday)*JDFACT;
frod += static_cast<double>(fday)*JDFACT*JDFACT;
long sod = static_cast<long>(frod*SEC_PER_DAY); // truncate
frod -= static_cast<double>(sod)/SEC_PER_DAY;
// fractional seconds of day
double frsod = frod*SEC_PER_DAY;

CommonTime ct;
return ct.set( jday,
static_cast<long>( sod ),
static_cast<double>( sod - static_cast<long>( sod ) ),
timeSystem );
return ct.set( jday, sod, frsod, timeSystem );
}
catch (InvalidParameter& ip)
{
Expand All @@ -73,14 +233,13 @@ namespace gpstk

void JulianDate::convertFromCommonTime( const CommonTime& ct )
{
long jday, sod;
long isod;
double fsod;
ct.get( jday, sod, fsod, timeSystem );
ct.get( jday, isod, fsod, timeSystem );

jd = static_cast<long double>( jday ) +
( static_cast<long double>( sod )
+ static_cast<long double>( fsod ) ) * DAY_PER_SEC
- 0.5;
double frac((static_cast<double>(isod)+fsod)/SEC_PER_DAY);
dday = static_cast<uint64_t>(frac/JDFACT);
fday = static_cast<uint64_t>((frac/JDFACT-static_cast<double>(dday))/JDFACT);
}

std::string JulianDate::printf( const std::string& fmt ) const
Expand All @@ -91,7 +250,7 @@ namespace gpstk
std::string rv( fmt );

rv = formattedPrint( rv, getFormatPrefixFloat() + "J",
"JLf", jd );
"JLf", JD() );
rv = formattedPrint( rv, getFormatPrefixInt() + "P",
"Ps", timeSystem.asString().c_str() );
return rv;
Expand Down Expand Up @@ -130,7 +289,7 @@ namespace gpstk
switch( i->first )
{
case 'J':
jd = asLongDouble( i->second );
this->fromString(i->second);
break;

case 'P':
Expand Down Expand Up @@ -159,22 +318,26 @@ namespace gpstk

void JulianDate::reset()
{
jd = 0.0;
jday = 0L;
dday = fday = static_cast<uint64_t>(0);
timeSystem = TimeSystem::Unknown;
}

bool JulianDate::operator==( const JulianDate& right ) const
{
/// Any (wildcard) type exception allowed, otherwise must be same time systems
if ((timeSystem != TimeSystem::Any &&
right.timeSystem != TimeSystem::Any) &&
timeSystem != right.timeSystem)
if((timeSystem != TimeSystem::Any && right.timeSystem != TimeSystem::Any)
&& timeSystem != right.timeSystem)
{
return false;

if( std::abs(jd - right.jd) < CommonTime::eps )
}
if(jday == right.jday
&& std::abs(((dday-right.dday)*JDFACT + (fday-right.fday)*JDFACT) * JDFACT)
< CommonTime::eps)
{
return true;
}

return false;
}

Expand All @@ -186,18 +349,15 @@ namespace gpstk
bool JulianDate::operator<( const JulianDate& right ) const
{
/// Any (wildcard) type exception allowed, otherwise must be same time systems
if ((timeSystem != TimeSystem::Any &&
right.timeSystem != TimeSystem::Any) &&
if (timeSystem != TimeSystem::Any && right.timeSystem != TimeSystem::Any &&
timeSystem != right.timeSystem)
{
gpstk::InvalidRequest ir("CommonTime objects not in same time system, cannot be compared");
gpstk::InvalidRequest ir("CommonTime objects not in same time system,"
" cannot be compared");
GPSTK_THROW(ir);
}

if( jd < right.jd )
{
return true;
}
if( jday < right.jday || dday < right.dday || fday < right.fday ) return true;
return false;
}

Expand All @@ -217,4 +377,23 @@ namespace gpstk
return ( !operator<( right ) );
}

// Compute long double JD
// Warning - precision lost when cast to double,
// and on systems where long double is implemented as double (eg. WIN)
// return long double Julian Date JD (not Julian day jday)
long double JulianDate::JD(void) const
{
long double jd = static_cast<long double>(jday);
// NB uint64_t can't do negative
if(dday < JDHALFDAY) {
jd += (static_cast<long double>(dday+JDHALFDAY)
+ static_cast<long double>(fday)*JDFACT)*JDFACT-1.0L;
}
else {
jd += (static_cast<long double>(dday-JDHALFDAY)
+ static_cast<long double>(fday)*JDFACT)*JDFACT;
}
return jd;
}

} // namespace

0 comments on commit beff4fb

Please sign in to comment.