johnpye/nrel/asc_sunpos.c

/*  ASCEND modelling environment
    Copyright (C) 2008 Carnegie Mellon University

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, 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 General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*//** @file
    Wrapper for sunpos_nrel.c to allow access to the calculation from ASCEND.
*/

/* include the external function API from libascend... */
#include <ascend/compiler/extfunc.h>

/* include error reporting API as well, so we can send messages to user */
#include <ascend/utilities/error.h>

/* for accessing the DATA instance */
#include <ascend/compiler/child.h>
#include <ascend/general/list.h>
#include <ascend/compiler/module.h>
#include <ascend/compiler/childinfo.h>
#include <ascend/compiler/parentchild.h>
#include <ascend/compiler/slist.h>
#include <ascend/compiler/type_desc.h>
#include <ascend/compiler/packages.h>
#include <ascend/compiler/symtab.h>
#include <ascend/compiler/instquery.h>
#include <ascend/compiler/instmacro.h>
#include <ascend/compiler/instance_types.h>

/* the code that we're wrapping... */
#include "spa.h"

#ifndef PI
# define PI 3.141592653589793238462
#endif

static ExtBBoxInitFunc sunpos_nrel_prepare;
static ExtBBoxFunc sunpos_nrel_calc;
static ExtBBoxInitFunc julian_day_nrel_prepare;
static ExtBBoxFunc julian_day_nrel_calc;

static const char *sunpos_nrel_help = "\
Calculate sun position using NREL SPA code. Inputs are:\n\
  * time (relative to reference time)\n\
  * pressure (instantaneous atmospheric pressure)\n\
  * temperature (instantaneous absolute atmospheric temperature)\n\
  * reference time (Julian Day value expressed as seconds)\n\
The reference time allows this function to use the same time variable as the\
rest of your simulation; the reference time is expected to be pre-calculated\
from a year-month-day calculation (see 'julian_day_nrel' external relation).";

static const char *julian_day_nrel_help = "Calculate the Julian Day from "
    "year, month, day, hour, minute, second and timezone inputs. "
    "Intended for once-off use in ASCEND models to calculate the time offset "
    "eg for the start of a weather file. Acceptable dates are in the range "
    "of -2000 BC to AD 6000. All of the inputs should be as 'factor' type "
    "variables (to avoid needless time unit conversions), except for the "
    "timezone, which should be in time units eg '8{h}'.";

/*------------------------------------------------------------------------------
  REGISTRATION FUNCTION
*/

/**
    This is the function called from 'IMPORT "johnpye/nrel/sunpos_nrels";'

    It sets up the functions contained in this external library
*/
extern
ASC_EXPORT int sunpos_nrel_register(){
    int result = 0;

    ERROR_REPORTER_HERE(ASC_USER_WARNING,"SUNPOS_NREL is still EXPERIMENTAL. Use with caution.");

#define CALCFN(NAME,INPUTS,OUTPUTS) \
    result += CreateUserFunctionBlackBox(#NAME \
        , NAME##_prepare \
        , NAME##_calc /* value */ \
        , (ExtBBoxFunc*)NULL /* derivatives not provided yet*/ \
        , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
        , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
        , INPUTS,OUTPUTS /* inputs, outputs */ \
        , NAME##_help /* help text */ \
        , 0.0 \
    ) /* returns 0 on success */

    CALCFN(sunpos_nrel,4,2);
    CALCFN(julian_day_nrel,7,1);

#undef CALCFN

    if(result){
        ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
    }
    return result;
}

#define GET_CHILD_VAL(NAME) \
    inst = ChildByChar(data,AddSymbol(#NAME)); \
    if(!inst){ \
        ERROR_REPORTER_HERE(ASC_USER_ERROR \
            ,"Couldn't locate '" #NAME "' in DATA, please check usage of SUNPOS."\
        );\
        return 1;\
    }\
    if(InstanceKind(inst)!=REAL_CONSTANT_INST){\
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member '" #NAME "' must be a real_constant");\
        return 1;\
    }\
    NAME = RC_INST(inst)->value;

/**
    This function is called when the black-box relation is being instantiated.

    This just gets the data member and checks that it's valid, and stores
    it in the blackbox data field.
*/
static int sunpos_nrel_prepare(struct BBoxInterp *bbox,
       struct Instance *data,
       struct gl_list_t *arglist
){
    struct Instance *inst;
    double latitude, longitude, elevation;

    /* get the latitude */
    GET_CHILD_VAL(latitude);
    CONSOLE_DEBUG("Latitude: %0.3f",latitude);
    if(latitude > PI/2 || latitude < -PI/2){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"'latitude' is out of allowable range -PI/2 to PI/2.");
        return 1;
    }

    /* get the longitude */
    GET_CHILD_VAL(longitude);
    CONSOLE_DEBUG("Longitude: %0.3f",longitude);
    if(longitude > PI || longitude < -PI){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"'latitude' is out of allowable range -PI to PI.");
        return 1;
    }

    /* get the elevation */
    GET_CHILD_VAL(elevation);
    CONSOLE_DEBUG("Elevation: %0.3f m",elevation);
    if(elevation < -6500000){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"'elevation' is out of allowable range (must be > -6,500 km)");
        return 1;
    }

    spa_data *S = ASC_NEW(spa_data);
    S->latitude = latitude * 180/PI;
    S->longitude = longitude * 180/PI;
    S->elevation = elevation;
    S->function = SPA_ZA_JD;

    ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared position for sun position.\n");
    bbox->user_data = (void *)S;
    return 0;
}

#define CALCPREPARE(NIN,NOUT) \
    /* a few checks about the input requirements */ \
    if(ninputs != NIN)return -1; \
    if(noutputs != NOUT)return -2; \
    if(inputs==NULL)return -3; \
    if(outputs==NULL)return -4; \
    if(bbox==NULL)return -5; \
    \
    /* the 'user_data' in the black box object will contain the */\
    /* coefficients required for this fluid; cast it to the required form: */\
    const spa_data *sunpos1 = (const spa_data *)bbox->user_data


/**
    Evaluation function for 'sunpos'
    @return 0 on success
*/
static int sunpos_nrel_calc(struct BBoxInterp *bbox,
        int ninputs, int noutputs,
        double *inputs, double *outputs,
        double *jacobian
){
    CALCPREPARE(4,2);

    double t, p, T, t_offset;

    t = inputs[0]; /* convert from JD seconds to JD days */
    p = inputs[1] / 100. /* convert Pa to mbar */;
    T = inputs[2] - 273.15 /* convert ��C to K */;
    t_offset = inputs[3];

    spa_data S = *sunpos1;
    S.pressure = p;
    S.temperature = T;
    S.jd = (t + t_offset) / 3600 / 24; /* convert to days */

    int res = spa_calculate(&S);
    //CONSOLE_DEBUG("Sun position: t = %f JD, p  %f mbar, T = %f C: res = %d, az = %f, zen = %f",S.jd, p, T, res, S.azimuth, S.zenith);

    /* returned values are in degrees, need to convert back to base SI: radians */
    outputs[0] = S.zenith * PI/180.;
    outputs[1] = S.azimuth180 * PI/180.;

    switch(res){
    case 0: break;
    case 16: CONSOLE_DEBUG("Calculated julian day (t + offset) = %f is out of permitted range",S.jd); break;
    default: CONSOLE_DEBUG("Error code %d returned from spa_calculate",res);
    }

    /* 0 on success, non-zero is error code from spa_calculate (would prob be input parameters out-of-range) */
    return res;
}

/*---------- SUNPOS_JULIAN_DAY ------------*/

static int julian_day_nrel_prepare(struct BBoxInterp *bbox,
       struct Instance *data,
       struct gl_list_t *arglist
){
    bbox->user_data = NULL;
    return 0;
}

static int julian_day_nrel_calc(struct BBoxInterp *bbox,
        int ninputs, int noutputs,
        double *inputs, double *outputs,
        double *jacobian
){
    CALCPREPARE(7,1);
    (void)sunpos1;

    int y,mon,d,h,m,s;
    double tz;

    y = inputs[0]; /* year */
    mon = inputs[1]; /* month */
    d = inputs[2]; /* day */
    h = inputs[3]; /* hour */
    m = inputs[4]; /* minute */
    s = inputs[5]; /* second */
    tz = inputs[6] / 3600.; /* timezone (in seconds, converted here to hours) */

    double t = julian_day(y,mon,d, h,m,s, tz) * 3600 * 24;

    outputs[0] = t;
    return 0;
}


1	jpye	2584	/* ASCEND modelling environment
2			Copyright (C) 2008 Carnegie Mellon University
3
4			This program is free software; you can redistribute it and/or modify
5			it under the terms of the GNU General Public License as published by
6			the Free Software Foundation; either version 2, or (at your option)
7			any later version.
8
9			This program is distributed in the hope that it will be useful,
10			but WITHOUT ANY WARRANTY; without even the implied warranty of
11			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12			GNU General Public License for more details.
13
14			You should have received a copy of the GNU General Public License
15	jpye	2649	along with this program. If not, see <http://www.gnu.org/licenses/>.
16	jpye	2584	//* @file
17			Wrapper for sunpos_nrel.c to allow access to the calculation from ASCEND.
18			*/
19
20			/* include the external function API from libascend... */
21			#include <ascend/compiler/extfunc.h>
22
23			/* include error reporting API as well, so we can send messages to user */
24			#include <ascend/utilities/error.h>
25
26			/* for accessing the DATA instance */
27			#include <ascend/compiler/child.h>
28			#include <ascend/general/list.h>
29			#include <ascend/compiler/module.h>
30			#include <ascend/compiler/childinfo.h>
31			#include <ascend/compiler/parentchild.h>
32			#include <ascend/compiler/slist.h>
33			#include <ascend/compiler/type_desc.h>
34			#include <ascend/compiler/packages.h>
35			#include <ascend/compiler/symtab.h>
36			#include <ascend/compiler/instquery.h>
37			#include <ascend/compiler/instmacro.h>
38			#include <ascend/compiler/instance_types.h>
39
40			/* the code that we're wrapping... */
41			#include "spa.h"
42
43			#ifndef PI
44	jpye	2599	# define PI 3.141592653589793238462
45	jpye	2584	#endif
46
47	jpye	2591	static ExtBBoxInitFunc sunpos_nrel_prepare;
48			static ExtBBoxFunc sunpos_nrel_calc;
49			static ExtBBoxInitFunc julian_day_nrel_prepare;
50			static ExtBBoxFunc julian_day_nrel_calc;
51	jpye	2584
52	jpye	2592	static const char *sunpos_nrel_help = "\
53			Calculate sun position using NREL SPA code. Inputs are:\n\
54			* time (relative to reference time)\n\
55			* pressure (instantaneous atmospheric pressure)\n\
56			* temperature (instantaneous absolute atmospheric temperature)\n\
57			* reference time (Julian Day value expressed as seconds)\n\
58			The reference time allows this function to use the same time variable as the\
59			rest of your simulation; the reference time is expected to be pre-calculated\
60			from a year-month-day calculation (see 'julian_day_nrel' external relation).";
61	jpye	2584
62	jpye	2591	static const char *julian_day_nrel_help = "Calculate the Julian Day from "
63			"year, month, day, hour, minute, second and timezone inputs. "
64			"Intended for once-off use in ASCEND models to calculate the time offset "
65			"eg for the start of a weather file. Acceptable dates are in the range "
66			"of -2000 BC to AD 6000. All of the inputs should be as 'factor' type "
67			"variables (to avoid needless time unit conversions), except for the "
68			"timezone, which should be in time units eg '8{h}'.";
69
70	jpye	2584	/*------------------------------------------------------------------------------
71			REGISTRATION FUNCTION
72			*/
73
74			/**
75	jpye	2591	This is the function called from 'IMPORT "johnpye/nrel/sunpos_nrels";'
76	jpye	2584
77			It sets up the functions contained in this external library
78			*/
79			extern
80			ASC_EXPORT int sunpos_nrel_register(){
81			int result = 0;
82
83	jpye	2857	ERROR_REPORTER_HERE(ASC_USER_WARNING,"SUNPOS_NREL is still EXPERIMENTAL. Use with caution.");
84	jpye	2584
85			#define CALCFN(NAME,INPUTS,OUTPUTS) \
86			result += CreateUserFunctionBlackBox(#NAME \
87	jpye	2591	, NAME##_prepare \
88	jpye	2584	, NAME##_calc /* value */ \
89			, (ExtBBoxFunc)NULL / derivatives not provided yet*/ \
90			, (ExtBBoxFunc)NULL / hessian not provided yet */ \
91			, (ExtBBoxFinalFunc)NULL / finalisation not implemented */ \
92			, INPUTS,OUTPUTS /* inputs, outputs */ \
93			, NAME##_help /* help text */ \
94			, 0.0 \
95			) /* returns 0 on success */
96
97	jpye	2592	CALCFN(sunpos_nrel,4,2);
98	jpye	2591	CALCFN(julian_day_nrel,7,1);
99	jpye	2584
100			#undef CALCFN
101
102			if(result){
103			ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
104			}
105			return result;
106			}
107
108			#define GET_CHILD_VAL(NAME) \
109			inst = ChildByChar(data,AddSymbol(#NAME)); \
110			if(!inst){ \
111			ERROR_REPORTER_HERE(ASC_USER_ERROR \
112			,"Couldn't locate '" #NAME "' in DATA, please check usage of SUNPOS."\
113			);\
114			return 1;\
115			}\
116			if(InstanceKind(inst)!=REAL_CONSTANT_INST){\
117			ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member '" #NAME "' must be a real_constant");\
118			return 1;\
119			}\
120			NAME = RC_INST(inst)->value;
121
122	jpye	2591	/**
123			This function is called when the black-box relation is being instantiated.
124
125			This just gets the data member and checks that it's valid, and stores
126			it in the blackbox data field.
127			*/
128			static int sunpos_nrel_prepare(struct BBoxInterp *bbox,
129			struct Instance *data,
130			struct gl_list_t *arglist
131			){
132			struct Instance *inst;
133			double latitude, longitude, elevation;
134
135	jpye	2584	/* get the latitude */
136			GET_CHILD_VAL(latitude);
137			CONSOLE_DEBUG("Latitude: %0.3f",latitude);
138			if(latitude > PI/2 \|\| latitude < -PI/2){
139			ERROR_REPORTER_HERE(ASC_USER_ERROR,"'latitude' is out of allowable range -PI/2 to PI/2.");
140			return 1;
141			}
142
143			/* get the longitude */
144			GET_CHILD_VAL(longitude);
145			CONSOLE_DEBUG("Longitude: %0.3f",longitude);
146			if(longitude > PI \|\| longitude < -PI){
147			ERROR_REPORTER_HERE(ASC_USER_ERROR,"'latitude' is out of allowable range -PI to PI.");
148			return 1;
149			}
150
151			/* get the elevation */
152			GET_CHILD_VAL(elevation);
153			CONSOLE_DEBUG("Elevation: %0.3f m",elevation);
154			if(elevation < -6500000){
155			ERROR_REPORTER_HERE(ASC_USER_ERROR,"'elevation' is out of allowable range (must be > -6,500 km)");
156			return 1;
157			}
158
159			spa_data *S = ASC_NEW(spa_data);
160			S->latitude = latitude * 180/PI;
161			S->longitude = longitude * 180/PI;
162			S->elevation = elevation;
163			S->function = SPA_ZA_JD;
164
165			ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared position for sun position.\n");
166			bbox->user_data = (void *)S;
167			return 0;
168			}
169
170			#define CALCPREPARE(NIN,NOUT) \
171			/* a few checks about the input requirements */ \
172			if(ninputs != NIN)return -1; \
173			if(noutputs != NOUT)return -2; \
174			if(inputs==NULL)return -3; \
175			if(outputs==NULL)return -4; \
176			if(bbox==NULL)return -5; \
177			\
178			/* the 'user_data' in the black box object will contain the */\
179			/* coefficients required for this fluid; cast it to the required form: */\
180			const spa_data sunpos1 = (const spa_data )bbox->user_data
181
182
183			/**
184			Evaluation function for 'sunpos'
185			@return 0 on success
186			*/
187	jpye	2591	static int sunpos_nrel_calc(struct BBoxInterp *bbox,
188	jpye	2584	int ninputs, int noutputs,
189			double inputs, double outputs,
190			double *jacobian
191			){
192	jpye	2592	CALCPREPARE(4,2);
193	jpye	2584
194	jpye	2592	double t, p, T, t_offset;
195	jpye	2584
196	jpye	2592	t = inputs[0]; /* convert from JD seconds to JD days */
197	jpye	2584	p = inputs[1] / 100. /* convert Pa to mbar */;
198			T = inputs[2] - 273.15 /* convert ��C to K */;
199	jpye	2592	t_offset = inputs[3];
200	jpye	2584
201			spa_data S = *sunpos1;
202			S.pressure = p;
203			S.temperature = T;
204	jpye	2592	S.jd = (t + t_offset) / 3600 / 24; /* convert to days */
205	jpye	2584
206			int res = spa_calculate(&S);
207	jpye	2599	//CONSOLE_DEBUG("Sun position: t = %f JD, p %f mbar, T = %f C: res = %d, az = %f, zen = %f",S.jd, p, T, res, S.azimuth, S.zenith);
208	jpye	2584
209			/* returned values are in degrees, need to convert back to base SI: radians */
210			outputs[0] = S.zenith * PI/180.;
211	jpye	2594	outputs[1] = S.azimuth180 * PI/180.;
212	jpye	2584
213	jpye	2599	switch(res){
214			case 0: break;
215			case 16: CONSOLE_DEBUG("Calculated julian day (t + offset) = %f is out of permitted range",S.jd); break;
216			default: CONSOLE_DEBUG("Error code %d returned from spa_calculate",res);
217			}
218
219	jpye	2584	/* 0 on success, non-zero is error code from spa_calculate (would prob be input parameters out-of-range) */
220			return res;
221			}
222
223	jpye	2591	/---------- SUNPOS_JULIAN_DAY ------------/
224	jpye	2584
225	jpye	2591	static int julian_day_nrel_prepare(struct BBoxInterp *bbox,
226			struct Instance *data,
227			struct gl_list_t *arglist
228			){
229			bbox->user_data = NULL;
230			return 0;
231			}
232	jpye	2584
233	jpye	2591	static int julian_day_nrel_calc(struct BBoxInterp *bbox,
234			int ninputs, int noutputs,
235			double inputs, double outputs,
236			double *jacobian
237			){
238			CALCPREPARE(7,1);
239			(void)sunpos1;
240
241			int y,mon,d,h,m,s;
242			double tz;
243
244	jpye	2599	y = inputs[0]; /* year */
245	jpye	2649	mon = inputs[1]; /* month */
246	jpye	2599	d = inputs[2]; /* day */
247			h = inputs[3]; /* hour */
248			m = inputs[4]; /* minute */
249			s = inputs[5]; /* second */
250			tz = inputs[6] / 3600.; /* timezone (in seconds, converted here to hours) */
251	jpye	2591
252			double t = julian_day(y,mon,d, h,m,s, tz) * 3600 * 24;
253	jpye	2649
254	jpye	2591	outputs[0] = t;
255			return 0;
256			}
257
258