MATLAB mex gateway function for NREL SOLPOS

solpos_mex.c

/*==========================================================
 * solpos.c - NREL solar position calculator
 *
 * calculates solar position from location and time
 *
 * The calling syntax is:
 *
 *		angles = solpos(location,datetime,weather)
 *
 * This is a MEX-file for MATLAB.
 *========================================================*/

#include "mex.h"

/* ===============================================================
 * INSERT THE NREL SOLPOS.C SOURCE CODE HERE
 * http://rredc.nrel.gov/solar/codesandalgorithms/solpos/solpos.c
 *
 * ALSO DOWNLOAD HEADER TO SAME FOLDER BEFORE COMPILING IN MATLAB
 * http://rredc.nrel.gov/solar/codesandalgorithms/solpos/solpos00.h
 *
 * BASED ON STEST00 FROM NREL
 * http://rredc.nrel.gov/solar/codesandalgorithms/solpos/stest00.c
 *
 * FOR MORE INFO SEE MAIN NREL WEBSITE
 * http://rredc.nrel.gov/solar/codesandalgorithms/solpos/
 * =============================================================== */

/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    
    /* variable declarations here */
    struct posdata pd, *pdat; /* declare a posdata struct and a pointer for
                               * it (if desired, the structure could be
                               * allocated dynamically with malloc) */
    long retval;              /* to capture S_solpos return codes */
    
    /* variable declarations here */
    /* inputs */
    double *location;       /* latitude, longitude, UTC timezone */
    double *datetime;       /* year, month, day, hour, minute, seconds */
    double *weather;        /* pressure, temperature */
    
    /* outputs */
    double *angles;         /* zenith, azimuth */
    
    /* check for proper number of arguments */
    if(nrhs!=3) {
        mexErrMsgIdAndTxt("DBM:solpos:nrhs","Three inputs required.");
    }
    if(nlhs!=1) {
        mexErrMsgIdAndTxt("DBM:solpos:nlhs","One output required.");
    }
    
    /* check that number of rows & columns in input arguments */
    if(mxGetM(prhs[0])!=1) {
        mexErrMsgIdAndTxt("DBM:solpos:notRowVector","Location should be row vector.");
    }
    if(mxGetN(prhs[0])!=3) {
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput","Location should be 3 column row vector with latitude, longitude, timezone.");
    }
    if(mxGetM(prhs[1])!=1) {
        mexErrMsgIdAndTxt("DBM:solpos:notRowVector","Datetime should be row vector.");
    }
    if(mxGetN(prhs[1])!=6) {
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput","Datetime should be 6 column row vector with year, month, day, hour, minute, second.");
    }
    if(mxGetM(prhs[2])!=1) {
        mexErrMsgIdAndTxt("DBM:solpos:notRowVector","Weather should be row vector.");
    }
    if(mxGetN(prhs[2])!=2) {
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput","Weather should be 6 column row vector with pressure and temperature.");
    }
    
    /* get pointer to the real data in the input matrix */
    location = mxGetPr(prhs[0]);
    datetime = mxGetPr(prhs[1]);
    weather = mxGetPr(prhs[2]);
    
    /* create the output matrix */
    plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL);
    
    /* get a pointer to the real data in the output matrix */
    angles = mxGetPr(plhs[0]);
    
    /* code here */

    /**************  Begin demo program **************/
    
    pdat = &pd; /* point to the structure for convenience */
    
    /* Initialize structure to default values. (Optional only if ALL input
     * parameters are initialized in the calling code, which they are not
     * in this example.) */
    S_init(pdat);
    
    pdat->function = ( (S_SOLAZM  | S_REFRAC) & ~S_DOY );
    
    pdat->timezone  =  location[2]; // -8; //UTC
    pdat->longitude = location[1]; //-119.2022; //degrees
    pdat->latitude  =  location[0]; //35.56836; //degrees
    
    pdat->press     = weather[0]; //1015.62055; //mbar
    pdat->temp      = weather[1]; //40.0; //deg C
    
    pdat->tilt      = 0;
    pdat->aspect    = 180;
    
    pdat->year      = datetime[0];
    //using ~S_DOY so it will calculate the day of year for me
    pdat->month = datetime[1];
    pdat->day   = datetime[2];
    
    pdat->hour      = datetime[3];
    pdat->minute    = datetime[4];
    pdat->second    = datetime[5];
    
    /* call the computational routine */
    retval = S_solpos(pdat);  /* S_solpos function call */
    //S_decode(retval, pdat);    /* ALWAYS look at the return code! */
    
    /*============================================================================
     * decodes the error codes from S_solpos return value
     *----------------------------------------------------------------------------*/
    
    if ( retval & (1L << S_YEAR_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the year: %d [1950-2050]\n",
                pdat->year);
    if ( retval & (1L << S_MONTH_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the month: %d\n",
                pdat->month);
    if ( retval & (1L << S_DAY_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the day-of-month: %d\n",
                pdat->day);
    if ( retval & (1L << S_DOY_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the day-of-year: %d\n",
                pdat->daynum);
    if ( retval & (1L << S_HOUR_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the hour: %d\n",
                pdat->hour);
    if ( retval & (1L << S_MINUTE_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the minute: %d\n",
                pdat->minute);
    if ( retval & (1L << S_SECOND_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the second: %d\n",
                pdat->second);
    if ( retval & (1L << S_TZONE_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the time zone: %f\n",
                pdat->timezone);
    if ( retval & (1L << S_INTRVL_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the interval: %d\n",
                pdat->interval);
    if ( retval & (1L << S_LAT_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the latitude: %f\n",
                pdat->latitude);
    if ( retval & (1L << S_LON_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the longitude: %f\n",
                pdat->longitude);
    if ( retval & (1L << S_TEMP_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the temperature: %f\n",
                pdat->temp);
    if ( retval & (1L << S_PRESS_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the pressure: %f\n",
                pdat->press);
    if ( retval & (1L << S_TILT_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the tilt: %f\n",
                pdat->tilt);
    if ( retval & (1L << S_ASPECT_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the aspect: %f\n",
                pdat->aspect);
    if ( retval & (1L << S_SBWID_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the shadowband width: %f\n",
                pdat->sbwid);
    if ( retval & (1L << S_SBRAD_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the shadowband radius: %f\n",
                pdat->sbrad);
    if ( retval & (1L << S_SBSKY_ERROR) )
        mexErrMsgIdAndTxt("DBM:solpos:incorrectInput", "S_decode ==> Please fix the shadowband sky factor: %f\n",
                pdat->sbsky);
    
    /* output off MEX function */
    angles[0] = pdat->zenref;
    angles[1] = pdat->azim;
    
}