cSlidingMass1.c

/***************************************************
 * Automatically generated by Maple.
 * Created On: Wed Mar 07 14:32:12 2012.
***************************************************/
#ifdef WMI_WINNT
#define EXP __declspec(dllexport)
#else
#define EXP
#endif
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/* EDIT: those are not needed
#include <mplshlib.h>
static MKernelVector kv;
EXP ALGEB M_DECL SetKernelVector(MKernelVector kv_in, ALGEB args) { kv=kv_in; return(kv->toMapleNULL()); }
*/

/***************************************************
* Variable Definition for System:

* State variable(s):
*    x[ 0] = `Main.'SlidingMass1::TF1'.v_rel`(t)
*    x[ 1] = `Main.Position.value`(t)
*
* Output variable(s):
*    y[ 0] = `Main.Position.value`(t)
*    y[ 1] = `Main.'SlidingMass1::Velocity'`(t)
*    y[ 2] = `Main.'SlidingMass1::M1'.a`(t)
*
* Input variable(s):
*    u[ 0] = `Main.'SlidingMass1::Force'`(t)
*
* Parameter(s):
*    p[ 0] = `Main.'SlidingMass1::M'` (default = 1.)
*
************************************************/

/* Fixed parameters */
#define NDIFF 2
#define NDFA 2
#define NEQ 8
#define NPAR 1
#define NINP 1
#define NDISC 0
#define NIX1 6
#define NOUT 3
#define NCON 0
#define NEVT 0
#ifdef EVTHYST
#define NZC 2*NEVT
#else
#define NZC NEVT
#endif

typedef struct {
  double h;    /* Integration step size */
  double *w;    /* Float workspace */
  long *iw;    /* Integer workspace */
  long err;    /* Error flag */
  char *buf;    /* Error message */
  /* EDIT: add a member which is the pointer to the input function */
  void (*inpfn)(double, double*);
} SolverStruct;

static void SolverError(SolverStruct *S, char *errmsg)
{
  sprintf(S->buf,"Error at t=%20.16e: %s\n",S->w[0],errmsg);
  /* EDIT: this is not needed
  if(S->err==-1) kv->error(S->buf);
  */
  S->err=1;
}

static double dsn_zero=0.0;
static unsigned char dsn_undefC[8] = { 0, 0, 0, 0, 0, 0, 0xF8, 0x7F };
static double *dsn_undef = (double *)&dsn_undefC;
static unsigned char dsn_posinfC[8] = { 0, 0, 0, 0, 0, 0, 0xF0, 0x7F };
static double *dsn_posinf = (double *)&dsn_posinfC;
static unsigned char dsn_neginfC[8] = { 0, 0, 0, 0, 0, 0, 0xF0, 0xFF };
static double *dsn_neginf = (double *)&dsn_neginfC;
#define trunc(v) ( (v>0.0) ? floor(v) : ceil(v) )


static void fp(long N, double T, double *Y, double *YP)
{
  double m, v;

  if( Y[0]>=0. ) {
    m = 1.;
    v = 1.;
  }
  else {
    m = 1.;
    v = -1.;
  }
  Y[4] = v/m;
  m = 1.;
  v = fabs(Y[0]);
  Y[5] = v/m;
  m = 1.;
  v = -Y[0];
  YP[1] = v/m;
  if( Y[5]<0.01 ) {
    m = 1.;
    v = 2452.41970901797979*Y[0];
  }
  else {
    m = 1.;
    v = 5.*Y[4]*(exp(-10.*Y[5])+4.)+0.001*Y[0];
  }
  Y[3] = v/m;
  m = Y[8];
  v = Y[7]+Y[3];
  Y[2] = v/m;
  m = 1.;
  v = -Y[2];
  YP[0] = v/m;
}

static void otp(double T, double *Y, double *YP)
{
  Y[6] = -Y[0];
}

/* EDIT: This function is not needed
static void inpfn(double T, double *U)
{
  if(T<1.0)
    U[0] = 50.;
  else
    U[0] = 0.;
}
*/

static void SolverUpdate(double *u, double *p, long first, long internal, SolverStruct *S)
{
  long i;
  /* EDIT: inpfn must be substituted with S->inpfn */
  (S->inpfn)(S->w[0],u);

  for(i=0; i<NINP; i++) S->w[i+NDIFF+NIX1-NINP+1]=u[i];
  if(p) for(i=0; i<NPAR; i++) S->w[i+NEQ+1]=p[i];
  fp(NEQ,S->w[0],&S->w[1],&S->w[NEQ+NPAR+1]);
  if(S->w[NEQ+NPAR+1]-S->w[NEQ+NPAR+1]!=0.0) {
    SolverError(S,"index-1 and derivative evaluation failure");
    return;
  }
  if(internal) return;
}

static void SolverOutputs(double *y, SolverStruct *S)
{
  otp(S->w[0],&S->w[1],&S->w[NEQ+NPAR+1]);
  y[ 0]=S->w[ 2];
  y[ 1]=S->w[ 7];
  y[ 2]=S->w[ 3];
}

static void EulerStep(double *u,SolverStruct *S)
{
  long i;
  S->w[0]+=S->h;
  for(i=1;i<=NDIFF;i++) S->w[i]+=S->h*S->w[NEQ+NPAR+i];
  SolverUpdate(u,NULL,0,0,S);
}

static void SolverSetup(double t0, double *ic, double *u,double *p, double *y, double h, SolverStruct *S)
{
  long i;

  S->h = h;
  S->iw=NULL;
  S->w[0] = t0;
  S->w[1] =  0.00000000000000000e+00;
  S->w[2] =  0.00000000000000000e+00;
  S->w[3] =  1.00000000000000000e+02;
  S->w[4] =  0.00000000000000000e+00;
  S->w[5] =  1.00000000000000000e+00;
  S->w[6] =  0.00000000000000000e+00;
  S->w[7] = -0.00000000000000000e+00;
  S->w[8] =  1.00000000000000000e+02;
  S->w[9] =  1.00000000000000000e+00;
  if(ic) for(i=0;i<NDIFF;i++) {
    S->w[i+1]=ic[i];
    S->w[i+NEQ+NPAR+1]=0.0;
  }
  SolverUpdate(u,p,1,0,S);
  SolverOutputs(y,S);
}

/*
  Parametrized simulation driver
*/
/* EDIT: M_DECL requires mplshlib.h, so has to be removed
   Also, in order to avoid use of inpfn() which blocks external interface
   (eg from ruby FFI), prototype has to be changed in order to accept 
   the pointer to the input function
*/
EXP long ParamDriverC(double t0, void (*inpfn)(double, double*), double dt, long npts, double *ic, double *p, double *out, char *errbuf, long internal)
{
  double u[NINP],y[NOUT],w[1+2*NEQ+NPAR+NDFA+NEVT];
  long i,j;
  SolverStruct S;

  /* Setup */
  for(i=0;i<npts*(NOUT+1);i++) out[i]=*dsn_undef;
  S.w=w;
  /* EDIT: Save the input function pointer into S */
  S.inpfn=inpfn;
  if(internal==0) S.err=0; else S.err=-1;
  S.buf=errbuf;
  SolverSetup(t0,ic,u,p,y,dt,&S);
  /* Output */
  out[0]=t0; for(j=0;j<NOUT;j++) out[j+1]=y[j];
  /* Integration loop */  
  for(i=1;i<npts;i++) {
    /* Take a step with states */    
    EulerStep(u,&S);
    if( S.err>0 ) break;
    /* Output */
    SolverOutputs(y,&S);
    out[i*(NOUT+1)]=S.w[0]; for(j=0;j<NOUT;j++) out[i*(NOUT+1)+j+1]=y[j];
  }

  return(i);
}

/* EDIT: this function is not necessary
EXP ALGEB M_DECL ParamDriver( MKernelVector kv_in, ALGEB *args )
*** REMOVED BODY
*/


/* 
EDIT: This part hast to be added
The original code generated by MapleSim uses a static function inpfn() to 
generate the inputs. This function gets called at each iteration step. Being 
hardcoded, it blocks the use of the rest as library or from FFI interfaces. 
So the above edits are meant to use a function pointer (stored into the S 
struct) that can customized and passed to the ParamDriverC() funtion.
*/
static void inputs(double T, double *U)
{
  if(T<1.0)
    U[0] = 50.;
  else
    U[0] = 0.;
}


int main(int argc, char *argv[])
{
  long i,j,iter;
  char errbuf[1000]; // buffer per ev messaggi di errore
  int npts;          // numero di iterazioni a passo fisso
  double t0, dt;     // Tempo iniziale e passo temporale
  double p[1];       // tebella parametri (qui uno sol)
  double ic[NDIFF] = {0.,0.}; // Initial Conditions
  double u[NINP];

  p[0] = 1.0;
  t0 = 0.0;
  dt = 0.01;
  if(argc<=1) {
    printf("Need one argument as number of points!\n");
    exit(1);
  }
  npts = atoi(argv[1]);
  double out[npts*(NOUT+1)+1]; // tabella dati in uscita
  u[0] = 50.;
  for(i=0;i<npts*(NOUT+1);i++) out[i]=*dsn_undef;
  iter = ParamDriverC(t0, &inputs, dt, npts, ic, p, out, errbuf, 0);
  printf("Result after %ld steps:\n",iter);
  printf("   i             T            s            v            a\n");
  for(i=0;i<npts;i++) {
    printf("%4ld ",i);
    for(j=0;j<NOUT+1;j++) printf("%12.4e ", out[i*(NOUT+1)+j]);
    printf("\n");
  }
  if(iter < npts) printf("Error: %s\n",errbuf);
  exit(0);
}

MSUpdater.rb

#!/usr/bin/env ruby

require "strscan"
TAB = "  "

class StringScanner
  def scan_and_put(match)
    puts scan_until(match).gsub(/\t/, TAB)
  end  
end

code = File.open(ARGV[0]) {|f| f.read}

scanner = StringScanner.new(code)

scanner.scan_and_put(/#include <math.h>/m)
scanner.scan_until(/toMapleNULL\(\)\); \}/)

scanner.scan_and_put(/char \*buf;\s*\/\* Error message \*\//m)
puts "  void (*inpfn)(double, double*); /* Input function pointer /*"

scanner.scan_and_put(/errmsg\);/)
scanner.scan_until(/error\(S->buf\);/)

scanner.scan_and_put(/static void otp/)
scanner.scan_and_put(/\}\n/)
scanner.scan_until(/static void inpfn/)
scanner.scan_until(/\}\n/)

scanner.scan_and_put(/long i;\n/)
scanner.scan_until(/for/)
puts "  (S->inpfn)(S->w[0],u);\n  for"

scanner.scan_and_put(/SolverOutputs\(y,S\);/)
puts "}\n\nEXP long ParamDriverC(double t0, void (*inpfn)(double, double*), double dt, long npts, double *ic, double *p, double *out, char *errbuf, long internal)\n{"

scanner.scan_and_put(/S\.w=w;/)
puts "  S.inpfn=inpfn;"
scanner.scan_and_put(/return\(i\);\n\}/)

The command:

ruby MSUpdater.rb <originalMapleSimCFile.c>

should give the shareable result.