src/sct_accuracy.c

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2007 by Aleksander Morgado Juez                         *
 *   scitime@aleksander_morgado.mm.st                                      *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU Library General Public License as       *
 *   published by the Free Software Foundation; either version 2 of the    *
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this program; if not, write to the                 *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "scitime.h"

static int _sct_get_exponent_from_double(const double number, int *pExponent);

SCT_RET sct_set_accuracy (sct_time_ptr pElem, SCT_ACCURACY accuracy)
{
    pElem->accuracy = accuracy;
    return SCT_RET_OK;
}


SCT_ACCURACY sct_get_accuracy (sct_time_ptr pElem)
{
    return pElem->accuracy;
}






SCT_ACCURACY sct_get_accuracy_from_jd(const double jd)
{
    int accuracyValue = 100;
    int exponent = 0;

    if( _sct_get_exponent_from_double(jd, &exponent) != SCT_RET_OK)
    {
        #ifdef SCT_DEBUG
            fprintf(stderr, "[scitime:%s] >> Could not get exponent...\n",
            __FUNCTION__);
        #endif
        return SCT_ACCURACY_DAYS;
    }

    // From http://en.wikipedia.org/wiki/IEEE_754
    //Double precision doubles structure:
    // - 1 bit for the sign (63)
    // - 11 bits for the exponent (52 to 62)
    // - 52 bit for the fraction

    //The number of correct digits in the whole double is given by:
    // 2^52 is = 4503599627370496 (16 digits) 
    // We will say it's 15 de maximum number of correct digits (skipping the
    //    last one)
    // So, the correct number of digits after the comma will be:
    //  c = 15 - exponent
    // For example, if the number is 10.5, exponent will be 1, so the correct
    //  number of digits after the comma will be 14
    // In this example, if the given value is a JD date, it says that its
    //  accuracy is correct until the 14th decimal value--> 10^-14 days!

    // Now accuracy must be converted to secs directly...
    //  1 sec = 1.15 E-5 days ... so... 5 full digits are due to getting
    //  the seconds of the number


    // So, finally... the correct number of digits AFTER the secs is 10
    accuracyValue = exponent -10;

    if(accuracyValue > SCT_ACCURACY_DAYS)
    {
        #ifdef SCT_DEBUG
            fprintf(stderr, "[scitime:%s] >> Accuracy is WORSE than days "
                "(1E%d secs)!\n", __FUNCTION__, accuracyValue);
        #endif
        return SCT_ACCURACY_DAYS;
    }
    if(accuracyValue < SCT_ACCURACY_FEMTOSECS)
    {
        #ifdef SCT_DEBUG
            fprintf(stderr, "[scitime:%s] >> Accuracy is BETTER than femtosecs "
                "(1E%d secs)! (Resetting to 1E-15 secs)\n",
            __FUNCTION__, accuracyValue);
        #endif
        return SCT_ACCURACY_FEMTOSECS;
    }

    return (SCT_ACCURACY)(accuracyValue);
}

SCT_ACCURACY sct_get_accuracy_from_secs(const double secs)
{
    int accuracyValue = 100;
    int exponent = 0;

    if( _sct_get_exponent_from_double(secs, &exponent) != SCT_RET_OK)
    {
        #ifdef SCT_DEBUG
            fprintf(stderr, "[scitime:%s] >> Could not get exponent...\n",
            __FUNCTION__);
        #endif
        return SCT_ACCURACY_DAYS;
    }

    // From http://en.wikipedia.org/wiki/IEEE_754
    //Double precision doubles structure:
    // - 1 bit for the sign (63)
    // - 11 bits for the exponent (52 to 62)
    // - 52 bit for the fraction

    //The number of correct digits in the whole double is given by:
    // 2^52 is = 4503599627370496 (16 digits) 
    // We will say it's 15 de maximum number of correct digits (skipping the
    //    last one)
    // So, the correct number of digits after the comma will be:
    //  c = 15 - exponent
    // For example, if the number is 10.5, exponent will be 1, so the correct
    //  number of digits after the comma will be 14

    // Now accuracy must be converted to secs directly... but it's already in
    //  seconds!! 

    // So... the correct number of digits is 15
    accuracyValue = exponent -15;

    if(accuracyValue > SCT_ACCURACY_DAYS)
    {
        #ifdef SCT_DEBUG
            fprintf(stderr, "[scitime:%s] >> Accuracy is WORSE than days "
                "(1E%d secs)!\n", __FUNCTION__, accuracyValue);
        #endif
        return SCT_ACCURACY_DAYS;
    }
    if(accuracyValue < SCT_ACCURACY_FEMTOSECS)
    {
        #ifdef SCT_DEBUG
            fprintf(stderr, "[scitime:%s] >> Accuracy is BETTER than femtosecs "
                "(1E%d secs)! (Resetting to 1E-15 secs)\n",
            __FUNCTION__, accuracyValue);
        #endif
        return SCT_ACCURACY_FEMTOSECS;
    }

    return (SCT_ACCURACY)(accuracyValue);
}


static int _sct_get_exponent_from_double(const double number, int *pExponent)
{
    char aux[50];
    char *expS = NULL;

    //Get value of the EXPONENT... can be done faster probably (TODO)
    memset(&aux[0],0,sizeof(char)*50);
    sprintf(aux,"%E",number);
    expS = strstr(aux,"E");
    if(expS == NULL)
    {
        #ifdef SCT_DEBUG
        fprintf(stderr, "[scitime:%s] >> Invalid string (%s) to get the "
            "exponent!\n", __FUNCTION__, aux);
        #endif
        return SCT_RET_INVALID_NUMBER;
    }
    else
    {
        *pExponent = atoi(&expS[1]);
    }
    return SCT_RET_OK;
}

---