core/unit.c

Go to the documentation of this file.

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <errno.h>
#include <ctype.h>
#include <math.h>

#include "platform.h"
#include "output.h"
#include "unit.h"
#include "exception.h"
#include "globals.h"
#include "find.h"

/* fundamental physical/economic constants */
static double c = 2.997925e8;       
static double e = 1.602189246e-19;  
static double h = 6.62617636e-34;   
static double k = 1.38066244e-23;   
static double m = 0.910953447e-30;  
static double s = 1.233270e4;       
static int precision = 10;          
static UNIT *unit_list = NULL;
static char filepath[1024] = "unitfile.txt";
static int linenum=0;

/* locate an underived unit */
UNIT *unit_find_underived(char *unit)
{
    UNIT *p;

    /* scan list for existing entry */
    for (p=unit_list; p!=NULL; p=p->next)
    {
        if (strcmp(p->name,unit)==0)
            return p;
    }
    return NULL;
}

/* define a constant */
int unit_constant(char name, double value)
{
    char buffer[2] = {name,'\0'};
    if (unit_find_underived(buffer)!=NULL)
        throw_exception("%s(%d): constant definition of '%s' failed; unit already defined", filepath, linenum, buffer);
    switch (name) {
    case 'c':
        c=value;
        break;
    case 'e':
        e=value;
        break;
    case 'h':
        h=value;
        break;
    case 'k':
        k=value;
        break;
    case 'm':
        m=value;
        break;
    case 's':
        s=value;
        break;
    default:
        throw_exception("%s(%d): constant '%c' is not valid", filepath,linenum,name);
        return 0;
    }
    return 1;
}

/* define a primary unit */
int unit_primary(char *name,double c,double e,double h,double k,double m,double s,double a,double b,int prec)
{
    UNIT *p = (UNIT*)malloc(sizeof(UNIT));
    if (unit_find_underived(name)!=NULL)
        throw_exception("%s(%d): primary definition of '%s' failed; unit already defined", filepath, linenum, name);
    if (p==NULL)
    {
        throw_exception("%s(%d): memory allocation failed",filepath,linenum);
        return 0;
    }
    strncpy(p->name,name,sizeof(p->name)-1);
    p->c = c;
    p->e = e;
    p->h = h;
    p->k = k;
    p->m = m;
    p->s = s;
    p->a = a;
    p->b = b;
    p->prec = prec;
    p->next = unit_list;
    unit_list = p;
    return 1;
}

/* determine the precision of a term */
int unit_precision(char *term)
{
    char* p;
    char mant[256];

    /* get mantissa */
    if ((p=strchr(term,'e'))!=NULL || (p=strchr(term,'E'))!=NULL) {

        /* has exponent */
        strncpy(mant, term, p-term);
    }
    else {
        strcpy(mant, term);
    }

    return (int)(strlen(mant) - (strchr(mant,'.')!=NULL ? 1 : 0));
}

/* define a derived unit */
int unit_derived(char *name,char *derivation)
{
    double c=0,e=0,h=0,k=0,m=0,s=0,a=0,b=0;
    int prec=0;
    double scalar=1.0;
    char lastOp = '\0', nextOp = '\0';
    UNIT *lastUnit = NULL;
    char *p=derivation;
    if (unit_find_underived(name)!=NULL)
        throw_exception("%s(%d): derived definition of '%s' failed; unit already defined", filepath, linenum, name);

    /* extract scalar */
    if (sscanf(p,"%lf",&scalar)==1)
    {
        /* advance pointer */
        if ((p=strchr(p,' '))!=NULL) 
            p++;

        /* reset pointer */
        else
            p=derivation;
    }

    /* scan for terms */
    while (*p!='\0')
    {
        char term[32];
        UNIT *pUnit;

        /* extract operation */
        if (sscanf(p,"%[^-*/^+]",term)!=1)
            throw_exception("%s(%d): unable to read unit '%s'", filepath,linenum,p);

        /* non-exponential ops */
        if (nextOp!='^')
        {
            /* find unit */
            pUnit = unit_find_underived(term);

            if (pUnit==NULL)
            {
                double val;
                if (sscanf(term,"%lf",&val)==1)
                {
                    static UNIT local;
                    if (nextOp=='+' || nextOp=='-') /* bias operation */
                    {
                        UNIT bias = {"",0,0,0,0,0,0,0,val,unit_precision(term)};
                        local = bias;
                    }
                    else /* scale operation */
                    {
                        UNIT offset = {"",0,0,0,0,0,0,val,0,unit_precision(term)};
                        local = offset;
                    }
                    pUnit = &local;
                }
                else
                    throw_exception("%s(%d): unable to find or derive unit '%s'", filepath,linenum,p);
            }
        }

        /* advance to next term */
        p += strlen(term);

        /* add this term to result */
        switch (nextOp) {
        case '\0':
            /* first term */
            if (a==0)
            {
                c = pUnit->c;
                e = pUnit->e;
                h = pUnit->h;
                k = pUnit->k;
                m = pUnit->m;
                s = pUnit->s;
                a = pUnit->a * scalar;
                b = pUnit->b;
                prec = pUnit->prec;
                nextOp = '*'; /* in case ^ follows first term */
            } /* else last term */
            break;
        case '*':
            c += pUnit->c;
            e += pUnit->e;
            h += pUnit->h;
            k += pUnit->k;
            m += pUnit->m;
            s += pUnit->s;
            a *= pUnit->a;
            prec = min(prec,pUnit->prec);
            break;
        case '/':
            c -= pUnit->c;
            e -= pUnit->e;
            h -= pUnit->h;
            k -= pUnit->k;
            m -= pUnit->m;
            s -= pUnit->s;
            a /= pUnit->a;
            prec = min(prec,pUnit->prec);
            break;
        case '+':
            b += pUnit->b;
            break;
        case '-':
            b -= pUnit->b;
            break;
        case '^':
            if (lastUnit!=NULL)
            {
                int repeat = 0;
                if (sscanf(term,"%d",&repeat)!=1 || repeat<2)
                    throw_exception("%s(%d): syntax error at '%s'", filepath,linenum,term);
                repeat--;
                switch(lastOp) {
                case '*':
                    c += lastUnit->c*repeat;
                    e += lastUnit->e*repeat;
                    h += lastUnit->h*repeat;
                    k += lastUnit->k*repeat;
                    m += lastUnit->m*repeat;
                    s += lastUnit->s*repeat;
                    a *= pow(lastUnit->a,repeat);
                    prec = min(prec,lastUnit->prec);
                    break;
                case '/':
                    c -= lastUnit->c*repeat;
                    e -= lastUnit->e*repeat;
                    h -= lastUnit->h*repeat;
                    k -= lastUnit->k*repeat;
                    m -= lastUnit->m*repeat;
                    s -= lastUnit->s*repeat;
                    a /= pow(lastUnit->a,repeat);
                    prec = min(prec,lastUnit->prec);
                    break;
                default:
                    throw_exception("%s(%d): ^ not allowed after '%c' at '%s'", filepath,linenum,lastOp,term);
                    return 0;
                }
            }
            break;
        default:
            throw_exception("%s(%d): '%c' is not recognized at '%s'", filepath,linenum,lastOp,term);
            return 0;
            break;
        }
        lastOp = nextOp;
        if ((nextOp=*p)!='\0')
            p++;
        lastUnit = pUnit;
    }

    return unit_primary(name,c,e,h,k,m,s,a,b,prec);
}

void unit_init(void)
{
    static int tried=0;
    char *glpath = getenv("GLPATH");
    FILE *fp = NULL;
    char *tpath = find_file(filepath, NULL, 4);
    /* try only once */
    if (tried)
        return;
    else
        tried = 1;
    
    if(tpath != NULL)
        fp = fopen(tpath, "r");

    /* locate unit file on GLPATH if not found locally */
    if (fp==NULL && glpath!=NULL)
    {
        char envbuf[1024];
        char *dir;
        strcpy(envbuf,glpath);
        dir = strtok(envbuf,";");
        while (dir!=NULL)
        {
            strcpy(filepath,dir);
            strcat(filepath,"/unitfile.txt");
            fp = fopen(filepath,"r");
            if (fp!=NULL)
                break;
            dir = strtok(NULL,";");
        }
    }

    /* unit file not found */
    if (fp==NULL)
    {
        output_error("%s: %s (GLPATH=%s)", filepath, strerror(errno),glpath);
        return;
    }

    /* scan file */
    while (!feof(fp) && !ferror(fp))
    {
        char buffer[1024], *p;
        if (fgets(buffer,sizeof(buffer),fp)==NULL)
            break;
        linenum++;

        /* wipe comments */
        p = strchr(buffer,';');
        if (p!=NULL)
            *p = '\0';

        /* remove trailing whitespace */
        p = buffer+strlen(buffer)-1;
        while (iswspace(*p) && p>buffer)
            *p--='\0';

        /* ignore blank lines or lines starting with white space*/
        if (buffer[0]=='\0' || iswspace(buffer[0]))
            continue;

        /* constant definition */
        if (buffer[0]=='#') {
            double value=0.0;
            char name;
            if (sscanf(buffer+1,"%c=%lf",&name,&value)!=2) 
                output_error("%s(%d): constant definition '%s' is not valid", filepath,linenum,buffer);
            else
                unit_constant(name,value);
            continue;
        }

        /* scalar definition */
        if (buffer[0]=='*')
        {
            output_warning("%s(%d): scalar '%s' ignored (not supported yet)", filepath,linenum,buffer+1);
            continue;
        }

        /* primary definition */
        {
            char name[256];
            double c, e, h, k, m, s, a, b;
            int prec;
            if (sscanf(buffer,"%[^=]=%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%d",name,&c,&e,&h,&k,&m,&s,&a,&b,&prec)==10) 
            {
                unit_primary(name,c,e,h,k,m,s,a,b,prec);
                continue;
            }
        }

        /* derived definition */
        {
            char name[256];
            char derivation[256];
            if (sscanf(buffer,"%[^=]=%[-+/*^0-9. %A-Za-z]",name,derivation)==2) 
            {
                unit_derived(name,derivation);
                continue;
            }
        }

        /* not recognized */
        throw_exception("%s(%d): unit '%s' is not recognized", filepath,linenum,buffer);
    }

    /* check status */
    if (ferror(fp))
        throw_exception("unitfile: %s", strerror(errno));
    else
        output_verbose("%s loaded ok", filepath);

    /* done */
    fclose(fp);
}

int unit_convert(char *from, char *to, double *pValue)
{
    UNIT *pFrom = unit_find(from);
    UNIT *pTo = unit_find(to);
    if (pFrom!=NULL && pTo!=NULL)
        return unit_convert_ex(pFrom,pTo,pValue);
    return 0;
}

int unit_convert_ex(UNIT *pFrom, UNIT *pTo, double *pValue)
{
    if (pTo->c==pFrom->c && pTo->e==pFrom->e && pTo->h==pFrom->h && pTo->k==pFrom->k && pTo->m==pFrom->m && pTo->s==pFrom->s)
    {
        *pValue = (*pValue-pFrom->b)*pFrom->a/pTo->a + pTo->b;
        return 1;
    }
    else
        return 0;
}

UNIT *unit_find(char *unit) 
{
    UNIT *p;

    /* first time */
    if (unit_list==NULL) unit_init();

    /* scan list for existing entry */
    for (p=unit_list; p!=NULL; p=p->next)
    {
        if (strcmp(p->name,unit)==0)
            return p;
    }

    /* derive entry if possible */
    if (unit_derived(unit,unit))
        return unit_list;
    else
        return NULL;
}

int unit_test(void)
{
    typedef struct {double value; char *unit;} VALUE;
    typedef struct {VALUE from; VALUE to; double precision;} TEST;
#ifndef PI
#define PI 3.141592635
#endif
    TEST test[] = {
        {1, "ratio",    1, "unit"},
        {1, "%",        0.01, "unit"},
        {1, "pu",       1, "1/unit"},
        {1, "/%",       1, "1/%"},

        {PI, "unit",    3.141592635, "unit"},
        {2*PI, "rad",   1, "unit"},
        {360, "deg",    1, "unit"},
        {400, "grad",   1, "unit"},
        {4, "quad",     1, "unit"},
        {4*PI, "sr",    1, "unit"},

        {1, "R",        5.0/9.0, "K", 0.001},
        {0, "degC",     273.14, "K", 0.001},
        {0, "degF",     459.65, "R", 0.001},
        {1, "kg",       1000, "g"},
        {1, "N",        1, "m*kg/s^2"},
        {1, "J",        1, "N*m"},
        
        {1, "min",      60,"s"},
        {1, "h",        60,"min"},
        {1, "day",      24,"h"},
        {1, "yr",       365,"day"},
        {1, "syr",      365.24,"day"},

        {1, "cm",       0.01, "m"},
        {1, "mm",       0.1, "cm"},
        {1, "km",       1000,"m"},
        {1, "in",       2.54, "cm"},
        {1, "ft",       12, "in"},
        {1, "yd",       3, "ft"},
        {1, "mile",     5280, "ft"},

        {1, "sf",       1,"ft^2"},
        {1, "sy",       1,"yd^2"},

        {1, "cf",       1,"ft^3"},
        {1, "cy",       1,"yd^3"},
        {1,"l",         0.001,"m^3"},
        {1,"gal",       3.785412, "l"},

        {2.2046, "lb",  1, "kg"},
        {1, "tonne",    1000, "kg"},
        
        {1, "mph",      1, "mile/h"},
        {1, "fps",      1, "ft/s"},
        {1, "fpm",      1, "ft/min"},
        {1, "gps",      1, "gal/s"},
        {1, "gpm",      1, "gal/min"},
        {1, "gph",      1, "gal/h"},
        {1, "cfm",      1, "ft^3/min"},
        {1, "ach",      1, "1/h"},
        
        {1, "Hz",       1, "1/s"},
        
        {1, "W",        1, "J/s"},
        {1, "kW",       1000, "W"},
        {1, "MW",       1000, "kW"},
        {1, "Wh",       1, "W*h"},
        {1, "kWh",      1000, "Wh"},
        {1, "MWh",      1000, "kWh"},
        {1, "Btu",      0.293, "Wh"},
        {1, "kBtu",     1000, "Btu"},
        {1, "MBtu",     1000, "kBtu"},
        {1, "ton",      12000, "Btu/h"},
        {1, "tons",     1, "ton*s"},
        {1, "tonh",     1, "ton*h"},
        {1, "hp",       746, "W"},
        {1, "W",        1, "V*A"},
        {1, "C",        1, "A*s"},
        {1, "F",        1, "C/V"},
        {1, "ohm",      1,"V/A"},
        {1, "H",        1,"ohm*s"},
    };
    int n, failed=0, succeeded=0;
    output_verbose("performing units tests");
    output_test("BEGIN: units tests");
    for (n=0; n<sizeof(test)/sizeof(test[0]); n++)
    {
        double v = test[n].from.value;
        if (test[n].precision==0)
            test[n].precision = 1e-4;

        /* forward test */
        if (!unit_convert(test[n].from.unit,test[n].to.unit,&v))
        {
            output_test("FAILED: conversion from %s to %s not possible", test[n].from.unit,test[n].to.unit);
            failed++;
        }
        else if (fabs(v-test[n].to.value)>test[n].precision)
        {
            output_test("FAILED: incorrect unit conversion %g %s -> %g %s (got %g %s instead)", 
                test[n].from.value,test[n].from.unit,
                test[n].to.value,test[n].to.unit,
                v,test[n].to.unit);
            failed++;
        }

        /* reverse test */
        else if (!unit_convert(test[n].to.unit,test[n].from.unit,&v))
        {
            output_test("FAILED: conversion from %s to %s not possible", test[n].to.unit,test[n].from.unit);
            failed++;
        }
        else if (fabs(v-test[n].from.value)>test[n].precision)
        {
            output_test("FAILED: incorrect unit conversion %g %s -> %g %s (got %g %s instead)", 
                test[n].to.value,test[n].to.unit,
                test[n].from.value,test[n].from.unit,
                v,test[n].from.unit);
            failed++;
        }

        else
        {
            output_test("SUCCESS: %g %s = %g %s",               
                    test[n].from.value,test[n].from.unit,
                    test[n].to.value,test[n].to.unit);
            succeeded++;
        }
    }
    output_test("END: %d units tested", n);
    output_debug("units tested: %d ok, %d failed (see '%s' for details).", succeeded, failed, global_testoutputfile);
    return failed;
}

---

GridLAB-D^TM Version 1.0

An open-source project initiated by the US Department of Energy