/*psn_mex/psn1d_mex.c*/
/***********************************************************************
    Finite Volume Poisson PDE Solver: C-Library & Matlab Toolbox
    Implements numerical solution of Poisson PDE
    in 2D  Cartesian and Cylindrical coordinates

    Copyright (C) 2004 Igor Kaufman
    Copyright (C) 2004 Lancaster University, UK

    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 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 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


    Author's email: i.kaufman@lancaster.ac.uk

    APPLICATION    : POISSON1D MATLAB MEXDLL
    VERSION        : 1.0
************************************************************************/
/*17.01.04*/

#include <string.h>
#include "mex.h"
#include "psn.h"

/*
MEX-DLL for Matlab.
Provides a solution for 1D Poisson equation in matter
*/

enum ERR_CODE {
  ERR_NOERROR,
  ERR_NOARG,
  ERR_NOOUT,
  ERR_CALC
};


void GetErrMsg(int ErrCode, char* buf, char* buf1);
void _export mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])
{


   int ErrCode=ERR_NOERROR, n, N,i;
   double *par=NULL, *v=NULL,*e=NULL, *u;
   char buf[20], buf1[100];
   psn_1d_struct pst;

__try {
   /*Read Input argument*/
   if (nrhs>=2) {
      for (i=0;i<nrhs;i++) {
        switch (i) {
          case 0:
            par=mxGetPr(prhs[i]);
            break;
          case 1:
            v=mxGetPr(prhs[i]);
            n = (mxGetM(prhs[i])*mxGetN(prhs[i]));
            break;
          case 2:
            e=mxGetPr(prhs[2]);
            break;
        }

      }
   }  else {
      ErrCode=ERR_NOARG;
      return;
   }

   if (nlhs)  {
     plhs[0]=mxCreateDoubleMatrix(n,1, mxREAL);
     u=mxGetPr(plhs[0]);
   } else {
      ErrCode=ERR_NOOUT;
     return;
   }

   //mexPrintf("n=%u\n",n);
   N=n-2;

   pst.N=N;
   pst.grid=par[0];
   pst.x[0]=par[1];
   pst.x[1]=par[2];
   pst.alpha[0]=par[3];
   pst.alpha[1]=par[4];
   pst.beta[0]=par[5];
   pst.beta[1]=par[6];
   pst.gamma[0]=par[7];
   pst.gamma[1]=par[8];
   pst.v=v;
   pst.u=u;
   pst.e=e;
   pst.a=(double*)mxCalloc(n,sizeof(double));
   pst.b=(double*)mxCalloc(n,sizeof(double));
   pst.c=(double*)mxCalloc(n,sizeof(double));
   if (psn_1d_solver(&pst)) {
      ErrCode=ERR_CALC;
      return;
   }
   //mexPrintf("here1\n");

} __finally {
     if (ErrCode) {
        GetErrMsg(ErrCode, buf, buf1);
        mexErrMsgTxt(buf1);
     }
}

}

void GetErrMsg(int ErrCode, char* buf, char* buf1)
{
   if (buf) {
     sprintf(buf,"%hd",(short int)ErrCode);
   }

   if (buf1)
   switch (ErrCode) {
     case ERR_NOARG:
       sprintf (buf1,"no input arguments");
       break;
     case ERR_NOOUT:
       sprintf(buf1,"no output arguments");
       break;
     case ERR_CALC:
       sprintf(buf1,"error in numerical routine");
       break;
    default:
       sprintf (buf1,"unknown error");
   }
}