climate/climate.cpp

Go to the documentation of this file.

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

#include "climate.h"
#include "solar_angles.h"

#define RAD(x) (x*PI)/180
                         //ET,CT,MT ,PT ,AT, HT
double std_meridians[] = {75,90,105,120,135,150};
double surface_angles[8] = {0,45,-45,90,-90,135,-135,180};

int tmy2_reader::open(const char *file){
    char *newname = gl_findfile((char *)file, NULL, 2);
    if (newname==NULL)
        return 0;
    fp = fopen(newname, "r");

    if(fp == NULL){
        return 0;
    }
    
    // read in the header (use the c code given in the manual)
    fgets(buf,500,fp);
    sscanf(buf,"%*s %s %s %d %*s %d %d %*s %d %d",data_city,data_state,&tz_offset,&lat_degrees,&lat_minutes,&long_degrees,&long_minutes);
    return 4;
}

int tmy2_reader::next(){
    // read the next line into the buffer using fgets.
    char *val;
    val = fgets(buf,500,fp);
    
    if(val != NULL)
        return 1;
    else
        return 0;
}

int tmy2_reader::header_info(char* city, char* state, int* degrees, int* minutes, int* long_deg, int* long_min){
    strcpy(city,data_city);
    strcpy(state,data_state);
    *degrees = lat_degrees;
    *minutes = lat_minutes;
    *long_deg = long_degrees;
    *long_min = long_minutes;
    return 1;
}

int tmy2_reader::read_data(double *dnr, double *dhr, double *tdb, double *rh, int* month, int* day, int* hour){
    int tmp_dnr, tmp_dhr, tmp_tdb, tmp_rh;
    //sscanf(buf, "%*2s%2d%2d%2d%*14s%4d%*2s%4d%*40s%4d%8*s%3d%*s",month,day,hour,&tmp_dnr,&tmp_dhr,&tmp_tdb,&tmp_rh);
    sscanf(buf, "%*2s%2d%2d%2d%*14s%4d%*2s%4d%*34s%4d%8*s%3d%*s",month,day,hour,&tmp_dnr,&tmp_dhr,&tmp_tdb,&tmp_rh);
                /* 3__5__7__9___23_27__29_33___67_71_79__82 */
    *dnr = tmp_dnr;
    *dhr = tmp_dhr;
    *tdb = ((double)tmp_tdb)/10.0;
    /* *tdb = ((double)tmp_tdb)/10.0 * 9.0 / 5.0 + 32.0; */
    *rh = ((double)tmp_rh)/100.0;
    return 1;
}

double tmy2_reader::calc_solar(COMPASS_PTS cpt, short doy, double lat, double sol_time, double dnr, double dhr,double vert_angle = 90)
{
    SolarAngles *sa = new SolarAngles();
    double surface_angle = surface_angles[cpt];
    double cos_incident = sa->cos_incident(lat,RAD(vert_angle),RAD(surface_angle),sol_time,doy);

    return (double)dnr * cos_incident + dhr;
}

void tmy2_reader::close(){
    fclose(fp);
}

CLASS *climate::oclass = NULL;
climate *climate::defaults = NULL;

climate::climate(MODULE *module)
{
    if (oclass==NULL)
    {
        oclass = gl_register_class(module,"climate",PC_BOTTOMUP); 
        if (gl_publish_variable(oclass,
            PT_char32, "city", PADDR(city),
            PT_char1024,"tmyfile",PADDR(tmyfile),
            PT_double,"temperature[degF]",PADDR(temperature),
            PT_double,"humidity[%]",PADDR(humidity),
            PT_double,"solar_flux[W/sf]",PADDR(solar_flux),
                PT_SIZE, 9,
            NULL)<1) GL_THROW("unable to publish properties in %s",__FILE__);
        strcpy(city,"");
        strcpy(tmyfile,"");
        temperature = 59.0;
        temperature_raw = 15.0;
        humidity = 0.75;
        //solar_flux = malloc(8 * sizeof(double));
        solar_flux[0] = solar_flux[1] = solar_flux[2] = solar_flux[3] = solar_flux[4] = solar_flux[5] = solar_flux[6] = solar_flux[7] = solar_flux[8] = 0.0; // W/sf normal
        //solar_flux_S = solar_flux_SE = solar_flux_SW = solar_flux_E = solar_flux_W = solar_flux_NE = solar_flux_NW = solar_flux_N = 0.0; // W/sf normal
        tmy = NULL;
        defaults = this;
    }
}

int climate::create(void) 
{
    memcpy(this,defaults,sizeof(*this));
    return 1;
}

int climate::init(OBJECT *parent)
{
    // ignore "" files
    if (strcmp(tmyfile,"")==0)
        return 1;
    
    // open access to the TMY file
    char *found_file = gl_findfile(tmyfile,getenv("GLPATH"),FF_READ);
    if (found_file==NULL || file.open(found_file)<3) // TODO: get proper values for solar
    {
        gl_error("weather file '%s' access failed", tmyfile);
        return 0;
    }

    // begin parsing the TMY file
    int line=0;
    tmy = (TMYDATA*)malloc(sizeof(TMYDATA)*8760);
    if (tmy==NULL)
    {
        gl_error("TMY buffer allocation failed");
        return 0;
    }

    int month,day, hour;
    double dnr,dhr;
    char cty[50];
    char st[3];
    int lat_deg,lat_min,long_deg,long_min;
    file.header_info(cty,st,&lat_deg,&lat_min,&long_deg,&long_min);
    double degrees = (double)lat_deg + ((double)lat_min) / 60;
    double long_degrees = (double)long_deg + ((double)long_min) / 60;
    SolarAngles *sa = new SolarAngles();
    double tz_meridian =  15 * abs(file.tz_offset);//std_meridians[-file.tz_offset-5];
    while (line<8760 && file.next())
    {
        
        file.read_data(&dnr,&dhr,&temperature,&humidity,&month,&day,&hour);
        int doy = sa->day_of_yr(month,day);
        int hoy = (doy - 1) * 24 + (hour-1); 
        if (hoy>=0 && hoy<8760){
            tmy[hoy].temp_raw = temperature;
            tmy[hoy].temp = temperature;
            if(0 == gl_convert("degC", "degF", &(tmy[hoy].temp))){
                gl_error("climate::init unable to gl_convert() 'degC' to 'degF'!");
                return 0;
            }
            tmy[hoy].rh = humidity;
            // calculate the solar radiation
            double sol_time = sa->solar_time((double)hour,doy,RAD(tz_meridian),RAD(long_degrees));
            double sol_rad = 0.0; 
            
            for(COMPASS_PTS c_point = CP_S; c_point < CP_LAST;c_point=COMPASS_PTS(c_point+1)){
                if(c_point == CP_H)
                    sol_rad = file.calc_solar(CP_E,doy,RAD(degrees),sol_time,dnr,dhr,0.0);//(double)dnr * cos_incident + dhr;
                else
                    sol_rad = file.calc_solar(c_point,doy,RAD(degrees),sol_time,dnr,dhr);//(double)dnr * cos_incident + dhr;
                /* TMY2 solar radiation data is in Watt-hours per square meter. */
                if(0 == gl_convert("W/m^2", "W/sf", &(sol_rad))){
                    gl_error("climate::init unable to gl_convert() 'W/m^2' to 'W/sf'!");
                    return 0;
                }
                tmy[hoy].solar[c_point] = sol_rad;
            }
            
        }
        else
            gl_error("%s(%d): day %d, hour %d is out of allowed range 0-8759 hours", tmyfile,line,day,hour);
        
        line++;
    }
    file.close();
    return 1;
}

TIMESTAMP climate::sync(TIMESTAMP t0) 
{
    if (tmy!=NULL)
    {
        DATETIME ts;
        gl_localtime(t0,&ts);
        int hoy = ts.yearday*24 + ts.hour;
        temperature = tmy[hoy].temp;
        temperature_raw = tmy[hoy].temp_raw;
        humidity = tmy[hoy].rh;
        memcpy(solar_flux,tmy[hoy].solar,CP_LAST*sizeof(double));
        return -(t0+(3600*TS_SECOND-t0%(3600 *TS_SECOND))); 
    }
    return TS_NEVER;
}

// IMPLEMENTATION OF CORE LINKAGE: climate

EXPORT int create_climate(OBJECT **obj, 
                          OBJECT *parent) 
{
    *obj = gl_create_object(climate::oclass,sizeof(climate));
    if (*obj!=NULL)
    {
        climate *my = OBJECTDATA(*obj,climate);
        gl_set_parent(*obj,parent);
        my->create();
        return 1;
    }
    return 0;
}

EXPORT int init_climate(OBJECT *obj, 
                        OBJECT *parent) 
{
    if (obj!=NULL)
    {
        climate *my = OBJECTDATA(obj,climate);
        return my->init(parent);
    }
    return 0;
}

EXPORT TIMESTAMP sync_climate(OBJECT *obj, 
                              TIMESTAMP t0) 
{
    TIMESTAMP t1 = OBJECTDATA(obj,climate)->sync(t0);
    obj->clock = t0;
    return t1;
}

csv::csv()
{
    fp=NULL;
    hdr=NULL;
    data=NULL;
    cols=0;
}

csv::~csv()
{
    if (fp) fclose(fp);
    if (hdr) delete [] hdr;
    if (data) delete [] data;
}

int csv::open(const char *file, const unsigned n, ...) // return the number of columns founds
{
    char *ff = gl_findfile((char*)file,NULL,FF_READ);
    if (ff==NULL) return 0;
    fp = fopen(ff,"r");
    if (fp==NULL) return 0;
    if (fgets(buffer,sizeof(buffer),fp)==NULL) return 0;
    char temp[sizeof(buffer)];
    strcpy(temp,buffer);
    hdr = new int[cols=n];
    if (hdr==NULL) return 0;

    // process each requested column
    va_list ptr;
    va_start(ptr,n);
    unsigned ncol;
    for (ncol=0; ncol<n; ncol++)
    {
        const char *arg=va_arg(ptr,const char*);
        char *p;
        int i=0;

        // scan the header list
        strcpy(buffer,temp);
        for (p=strtok(buffer,","); p!=NULL; p=strtok(NULL,","),i++)
        {
            if (strcmp(p,arg)==0)
            {
                hdr[ncol] = i;
                break;
            }
        }

        // didn't find it
        if (p==NULL) 
            return ncol;
    }
    data = new char*[ncol];
    return ncol;
}

int csv::next(void)
{
    if (fgets(buffer,sizeof(buffer),fp)==NULL) return 0;
    unsigned n, i=0;
    char *p;
    for (p=strtok(buffer,","); p!=NULL; p=strtok(NULL,","),i++)
    {
        for (n=0; n<cols; n++)
        {
            if (hdr[n]==i)
            {
                data[n]=p;
                break;
            }
        }
    }
    return cols;
}

char *csv::read(unsigned n)
{
    if (n>=0 && n<cols)
        return data[n];
    else
        return NULL;
}

---

GridLAB-D^TM Version 1.0

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