Computational routine
eng

##
fromws_c

### File content

    #include "scicos_block4.h"
    /*    Masoud Najafi, Alan Layec September 2007 */
    /*    Copyright INRIA
     *    Scicos block simulator
     *    From workspace block
     */
    #include "../stack-c.h"
    #include <stdio.h>
    #include <string.h>
    #include "../machine.h"
    #include <math.h>

    #if WIN32
    #define NULL    0
    #endif

    #define Fnlength  block->ipar[0]
    #define FName     block->ipar[1]
    #define Method     block->ipar[1+Fnlength]
    #define ZC        block->ipar[2+Fnlength]
    #define OutEnd    block->ipar[3+Fnlength]
    #define T0        ptr->workt[0]
    #define TNm1      ptr->workt[nPoints-1]
    #define TP        (TNm1-0)

    extern int C2F(cvstr) __PARAMS((integer *,integer *,char *,integer *,unsigned long int));
    extern int C2F(mgetnc)();
    extern void C2F(mopen)();
    extern int C2F(cluni0) __PARAMS((char *name, char *nams, integer *ln, long int name_len,
                                    long int nams_len));
    extern void C2F(mclose) __PARAMS((integer *fd, double *res));
    extern void sciprint __PARAMS((char *fmt,...));
    int Mytridiagldltsolve(double *d__, double * l, double * b, int n);
    int Myevalhermite2(double *t, double *xa, double *xb, double *ya, double *yb, double *da, double *db, double *h, double *dh, double *ddh, double *dddh, int *i);

    //int Myevalhermite(double *t, double *xa, double *xb, double *ya, double *yb, double *da, double *db, double *h, double *dh, double *ddh, double *dddh, int *i);

    static char fmtd[3]={'d','l','\000'};
    static char fmti[3]={'i','l','\000'};
    static char fmtl[3]={'l','l','\000'};
    static char fmts[3]={'s','l','\000'};
    static char fmtc[3]={'c','l','\000'};
    static char fmtul[3]={'u','l','\000'};
    static char fmtus[3]={'u','s','\000'};
    static char fmtuc[3]={'u','c','\000'};

    #ifdef hppa
    #undef FILENAME_MAX
    #define FILENAME_MAX 4096
    #endif
    /* work struct for that block */
    typedef struct {
      int nPoints;
      int Yt;
      int Yst;
      int cnt1;
      int cnt2;
      int EVindex;
      int PerEVcnt;
      int firstevent;
      double *D;
      void *work;
      double *workt;
    } fromwork_struct ;

    void fromws_c(scicos_block *block,int flag)
    {
      double t,y1,y2,t1,t2,r;
      int i,inow;
      double *spline, *A_d, *A_sd, *qdy;
      /* double  a,b,c,*y;*/
      double d1,d2,h, dh, ddh, dddh;
      int fd;
      char *status;
      int swap = 1;
      double res;
      int out_n;
      long int lout;
      char filename[FILENAME_MAX];
      char str[100];
      int ierr;
      int Ytype, nPoints, mY, YsubType, my, ytype,j,jfirst;
      int Ydim[10];
      int cnt1, cnt2, EVindex, PerEVcnt;
      /* generic pointer */
      SCSREAL_COP *y_d,*y_cd,*ptr_d, *ptr_T, *ptr_D;
      SCSINT8_COP *y_c,*ptr_c;
      SCSUINT8_COP *y_uc, *ptr_uc;
      SCSINT16_COP *y_s,*ptr_s;
      SCSUINT16_COP *y_us,*ptr_us;
      SCSINT32_COP *y_l,*ptr_l;
      SCSUINT32_COP *y_ul,*ptr_ul;
     /* the struct ptr of that block */
      fromwork_struct *ptr;

     /* for path of TMPDIR/workspace */
      char env[256];
      char sep[2];
    #ifdef _MSC_VER
      sep[0]='\\';
    #else
      sep[0]='/';
    #endif
      sep[1]='\0';

     my=GetOutPortRows(block,1); /* number of rows of Outputs*/

     if (flag==4){
       C2F(cvstr)(&(Fnlength),&(FName),str,(j=1,&j),(unsigned long)strlen(str));
       str[Fnlength] = '\0';

       /* retrieve path of TMPDIR/workspace */
       strcpy(env,getenv("TMPDIR"));
       strcat(env,sep);
       strcat(env,"Workspace");
       strcat(env,sep);
       strcat(env,str);

       /* open tmp file */
       status = "r";
       lout=FILENAME_MAX;
       C2F(cluni0)(env, filename, &out_n,1,lout);
       C2F(mopen)(&fd,env,status,&swap,&res,&ierr);

       if (ierr!=0) {
         sciprint("The '%s' variable does not exist.\n",str);
         set_block_error(-3);
         return;
       }

       /* read x */
       C2F(mgetnc) (&fd, &Ydim[0], (j=nsiz,&j), fmti, &ierr);  /* read sci id */
       C2F(mgetnc) (&fd, &Ydim[6], (j=1,&j), fmti, &ierr);  /* read sci type */
       C2F(mgetnc) (&fd, &Ydim[7], (j=3,&j), fmti, &ierr);  /* read sci header */

       Ytype=Ydim[6];
       nPoints=Ydim[7];
       mY=Ydim[8];
       YsubType=Ydim[9];
       ytype=GetOutType(block,1);     /* output type */

       if (mY!=my) {
         sciprint("Data dimensions are inconsistent:\n\r Variable size=%d \n\r Block output size=%d.\n",mY,my);
         set_block_error(-3);
         return;
       }

       if (Ytype==1) { /*real/complex cases*/
         switch (YsubType)
         {
         case 0: if (ytype!=10) {
                   sciprint("Output should be of Real type.\n");
                   set_block_error(-3);
                   C2F(mclose)(&fd,&res);
                   return;
                 }
                 break;

         case 1: if (ytype!=11) {
                   sciprint("Output should be of complex type.\n");
                   set_block_error(-3);
                   C2F(mclose)(&fd,&res);
                   return;
                 }
                 break;
         }
       }
       else if(Ytype==8) { /*integer cases*/
         switch (YsubType)
         {
         case 1: if (ytype!=81) {
                   sciprint("Output should be of int8 type.\n");
                   set_block_error(-3);
                   C2F(mclose)(&fd,&res);
                   return;
                 }
                 break;

         case 2: if (ytype!=82) {
                   sciprint("Output should be of int16 type.\n");
                   set_block_error(-3);
                   C2F(mclose)(&fd,&res);
                   return;
                 }
                 break;

         case 4: if (ytype!=84) {
                   sciprint("Output should be of int32 type.\n");
                   set_block_error(-3);
                   C2F(mclose)(&fd,&res);
                   return;
                 }
                 break;

         case 11:if (ytype!=811) {
                   sciprint("Output should be of uint8 type.\n");
                   set_block_error(-3);
                   C2F(mclose)(&fd,&res);
                   return;
                 }
                 break;

         case 12:if (ytype!=812) {
                   sciprint("Output should be of uint16 type.\n");
                   set_block_error(-3);
                   C2F(mclose)(&fd,&res);
                   return;
                 }
                 break;

         case 14:if (ytype!=814) {
                   sciprint("Output should be of uint32 type.\n");
                   set_block_error(-3);
                   C2F(mclose)(&fd,&res);
                   return;
                 }
                 break;
         }
       }
       if((*(block->work)=(fromwork_struct*) scicos_malloc(sizeof(fromwork_struct)))==NULL) {
         set_block_error(-16);
         C2F(mclose)(&fd,&res);
         return;
       }
       ptr = *(block->work);
       ptr->D=NULL;
       ptr->workt=NULL;
       ptr->work=NULL;

       if (Ytype==1) { /*real/complex case*/
         switch (YsubType) {
         case 0 : // Real
           if((ptr->work=(void *) scicos_malloc(nPoints*mY*sizeof(double)))==NULL) {
             set_block_error(-16);
             scicos_free(ptr);
             *(block->work)=NULL;
             C2F(mclose)(&fd,&res);
             return;
           }
           ptr_d = (SCSREAL_COP *) ptr->work;
           C2F(mgetnc) (&fd, ptr_d, (j=nPoints*mY,&j), fmtd, &ierr);  /* read double data */
           break;
         case 1:  // complex
           if((ptr->work=(void *) scicos_malloc(2*nPoints*mY*sizeof(double)))==NULL) {
             set_block_error(-16);
             scicos_free(ptr);
             *(block->work)=NULL;
             C2F(mclose)(&fd,&res);
             return;
           }
           ptr_d = (SCSREAL_COP *) ptr->work;
           C2F(mgetnc) (&fd, ptr_d, (j=2*nPoints*mY,&j), fmtd, &ierr);  /* read double data */
           break;
         }
       }else if(Ytype==8) { /*integer case*/
         switch (YsubType) {
         case 1 ://int8
           if((ptr->work=(void *) scicos_malloc(nPoints*mY*sizeof(char)))==NULL) {
             set_block_error(-16);
             scicos_free(ptr);
             *(block->work)=NULL;
             C2F(mclose)(&fd,&res);
             return;
           }
           ptr_c = (SCSINT8_COP *) ptr->work;
           C2F(mgetnc) (&fd, ptr_c, (j=nPoints*mY,&j), fmtc, &ierr);  /* read char data */
           break;
         case 2 :  // int16
           if((ptr->work=(void *) scicos_malloc(nPoints*mY*sizeof(short)))==NULL) {
             set_block_error(-16);
             scicos_free(ptr);
             *(block->work)=NULL;
             C2F(mclose)(&fd,&res);
             return;
           }
           ptr_s = (SCSINT16_COP *) ptr->work;
           C2F(mgetnc) (&fd, ptr_s, (j=nPoints*mY,&j), fmts, &ierr);  /* read short data */
           break;
         case 4 :   // int32
           if((ptr->work=(void *) scicos_malloc(nPoints*mY*sizeof(long)))==NULL) {
             set_block_error(-16);
             scicos_free(ptr);
             *(block->work)=NULL;
             C2F(mclose)(&fd,&res);
             return;
           }
           ptr_l = (SCSINT32_COP *) ptr->work;
           C2F(mgetnc) (&fd, ptr_l, (j=nPoints*mY,&j), fmtl, &ierr);  /* read short data */
           break;
         case 11 :   // uint8
           if((ptr->work=(void *) scicos_malloc(nPoints*mY*sizeof(unsigned char)))==NULL) {
             set_block_error(-16);
             scicos_free(ptr);
             *(block->work)=NULL;
             C2F(mclose)(&fd,&res);
             return;
           }
           ptr_uc = (SCSUINT8_COP *) ptr->work;
           C2F(mgetnc) (&fd, ptr_uc, (j=nPoints*mY,&j), fmtuc, &ierr);  /* read short data */
           break;
         case 12 : // uint16
           if((ptr->work=(void *) scicos_malloc(nPoints*mY*sizeof(unsigned short)))==NULL) {
             set_block_error(-16);
             scicos_free(ptr);
             *(block->work)=NULL;
             C2F(mclose)(&fd,&res);
             return;
           }
           ptr_us = (SCSUINT16_COP *) ptr->work;
           C2F(mgetnc) (&fd, ptr_us, (j=nPoints*mY,&j), fmtus, &ierr);  /* read short data */
           break;
         case 14 :  // uint32
           if((ptr->work=(void *) scicos_malloc(nPoints*mY*sizeof(unsigned long)))==NULL) {
             set_block_error(-16);
             scicos_free(ptr);
             *(block->work)=NULL;
             C2F(mclose)(&fd,&res);
             return;
           }
           ptr_ul = (SCSUINT32_COP *) ptr->work;
           C2F(mgetnc) (&fd, ptr_ul, (j=nPoints*mY,&j), fmtul, &ierr);  /* read short data */
           break;
         }
       }

       /* read t */
       C2F(mgetnc) (&fd, &Ydim[0], (j=nsiz,&j), fmti, &ierr);  /* read sci id */
       C2F(mgetnc) (&fd, &Ydim[6], (j=1,&j), fmti, &ierr);  /* read sci type */
       C2F(mgetnc) (&fd, &Ydim[7], (j=3,&j), fmti, &ierr);  /* read sci header */

       if (nPoints!=Ydim[7]) {
         sciprint("The size of the Time(%d) and Data(%d) vectors are inconsistent.\n",Ydim[7],nPoints);
         set_block_error(-3);
         *(block->work)=NULL;
         scicos_free(ptr->work);
         scicos_free(ptr);
         C2F(mclose)(&fd,&res);
         return;
       }

       if ((Ydim[6]!=1) | (Ydim[9]!=0)) {
         sciprint("The Time vector type is not ""double"".\n");
         set_block_error(-3);
         *(block->work)=NULL;
         scicos_free(ptr->work);
         scicos_free(ptr);
         C2F(mclose)(&fd,&res);
         return;
       }

       if((ptr->workt=(double *) scicos_malloc(nPoints*sizeof(double)))==NULL) {
         set_block_error(-16);
         *(block->work)=NULL;
         scicos_free(ptr->work);
         scicos_free(ptr);
         C2F(mclose)(&fd,&res);
         return;
       }
       ptr_T = (SCSREAL_COP *) ptr->workt;
       C2F(mgetnc) (&fd, ptr_T, (j=nPoints,&j), fmtd, &ierr);  /* read data of t */
       /* close the file*/
       C2F(mclose)(&fd,&res);

       //================================
       // check for an increasing time data
       for(j = 0; j < nPoints-1; j++) {
         if(ptr_T[j] > ptr_T[j+1]) {
           sciprint("The time vector should be an increasing vector.\n");
           set_block_error(-3);
           *(block->work)=NULL;
           scicos_free(ptr->workt);
           scicos_free(ptr->work);
           scicos_free(ptr);
           return;
         }
       }
       //=================================
       if ((Method>1)&&((Ytype==1)&&( YsubType==0))) {/* double */

         if((ptr->D=(double *) scicos_malloc(nPoints*mY*sizeof(double)))==NULL) {
           set_block_error(-16);
           *(block->work)=NULL;
           scicos_free(ptr->workt);
           scicos_free(ptr->work);
           scicos_free(ptr);
           return;
         }

         if((spline=(double *) scicos_malloc((3*nPoints-2)*sizeof(double)))==NULL){
           sciprint("Allocation problem in spline.\n");
           set_block_error(-16);
           *(block->work)=NULL;
           scicos_free(ptr->D);
           scicos_free(ptr->workt);
           scicos_free(ptr->work);
           scicos_free(ptr);
           return;
         }
         A_d=spline;
         A_sd=A_d+nPoints;
         qdy=A_sd+nPoints-1;

         for (j=0;j<mY;j++){

           for (i=0;i<=nPoints-2;i++){
    	 A_sd[i] = 1.0 / (ptr_T[i+1] - ptr_T[i]);
    	 qdy[i] = (ptr_d[i+1+j*nPoints] - ptr_d[i+j*nPoints]) * A_sd[i]*A_sd[i];
           }

           for (i=1;i<=nPoints-2;i++){
    	 A_d[i] = 2.0*(A_sd[i-1] +A_sd[i]);
    	 ptr->D[i+j*nPoints] = 3.0*(qdy[i-1]+qdy[i]);
           }

           if (Method==2){
    	 A_d[0] =  2.0*A_sd[0];
    	 ptr->D[0+j*nPoints] = 3.0 * qdy[0];
             A_d[nPoints-1] =  2.0*A_sd[nPoints-2];
             ptr->D[nPoints-1+j*nPoints] =  3.0 * qdy[nPoints-2];
    	 Mytridiagldltsolve(A_d, A_sd, &ptr->D[j*nPoints], nPoints);
           }

           if (Method==3){
    	 /*  s'''(x(2)-) = s'''(x(2)+) */
    	 r = A_sd[1]/A_sd[0];
    	 A_d[0]= A_sd[0]/(1.0 + r);
    	 ptr->D[j*nPoints]=((3.0*r+2.0)*qdy[0]+r*qdy[1])/((1.0+r)*(1.0+r));
    	 /*  s'''(x(n-1)-) = s'''(x(n-1)+) */
    	 r = A_sd[nPoints-3]/A_sd[nPoints-2];
    	 A_d[nPoints-1] = A_sd[nPoints-2]/(1.0 + r);
    	 ptr->D[nPoints-1+j*nPoints] = ((3.0*r+2.0)*qdy[nPoints-2]+r*qdy[nPoints-3])/((1.0+r)*(1.0+r));
    	 Mytridiagldltsolve(A_d, A_sd, &ptr->D[j*nPoints], nPoints);
           }
         }
         scicos_free(spline);
         /*-----------------------------------------------------------*/
       }else if ((Method>1)&&((Ytype==1)&&( YsubType==1))) {/* Complex */

        if((ptr->D=(double *) scicos_malloc(2*nPoints*mY*sizeof(double)))==NULL) {
           set_block_error(-16);
           *(block->work)=NULL;
           scicos_free(ptr->workt);
           scicos_free(ptr->work);
           scicos_free(ptr);
           return;
         }

         if((spline=(double *) scicos_malloc((3*nPoints-2)*sizeof(double)))==NULL){
           sciprint("Allocation problem in spline");
           set_block_error(-16);
           *(block->work)=NULL;
           scicos_free(ptr->D);
           scicos_free(ptr->workt);
           scicos_free(ptr->work);
           scicos_free(ptr);
           return;
         }
         A_d=spline;
         A_sd=A_d+nPoints;
         qdy=A_sd+nPoints-1;
         for (j=0;j<mY;j++){
           for (i=0;i<=nPoints-2;i++){
    	 A_sd[i] = 1.0 / (ptr_T[i+1] - ptr_T[i]);
    	 qdy[i] = (ptr_d[i+1+j*nPoints] - ptr_d[i+j*nPoints]) * A_sd[i]*A_sd[i];
           }

           for (i=1;i<=nPoints-2;i++){
    	 A_d[i] = 2.0*(A_sd[i-1] +A_sd[i]);
    	 ptr->D[i+j*nPoints] = 3.0*(qdy[i-1]+qdy[i]);
           }

           if (Method==2){
    	 A_d[0] =  2.0*A_sd[0];
    	 ptr->D[0+j*nPoints] = 3.0 * qdy[0];
             A_d[nPoints-1] =  2.0*A_sd[nPoints-2];
             ptr->D[nPoints-1+j*nPoints] =  3.0 * qdy[nPoints-2];
    	 Mytridiagldltsolve(A_d, A_sd, &ptr->D[j*nPoints], nPoints);
           }

           if (Method==3){
    	 /*  s'''(x(2)-) = s'''(x(2)+) */
    	 r = A_sd[1]/A_sd[0];
    	 A_d[0]= A_sd[0]/(1.0 + r);
    	 ptr->D[+j*nPoints]=((3.0*r+2.0)*qdy[0]+r*qdy[1])/((1.0+r)*(1.0+r));
    	 /*  s'''(x(n-1)-) = s'''(x(n-1)+) */
    	 r = A_sd[nPoints-3]/A_sd[nPoints-2];
    	 A_d[nPoints-1] = A_sd[nPoints-2]/(1.0 + r);
    	 ptr->D[nPoints-1+j*nPoints] = ((3.0*r+2.0)*qdy[nPoints-2]+r*qdy[nPoints-3])/((1.0+r)*(1.0+r));
    	 Mytridiagldltsolve(A_d, A_sd, &ptr->D[j*nPoints], nPoints);
           }
           /* ********* */
           for (i=0;i<=nPoints-2;i++){
    	 A_sd[i] = 1.0 / (ptr_T[i+1] - ptr_T[i]);
    	 qdy[i] = (ptr_d[nPoints+i+1+j*nPoints] - ptr_d[nPoints+i+j*nPoints]) * A_sd[i]*A_sd[i];
           }

           for (i=1;i<=nPoints-2;i++){
    	 A_d[i] = 2.0*(A_sd[i-1] +A_sd[i]);
    	 ptr->D[i+j*nPoints+nPoints] = 3.0*(qdy[i-1]+qdy[i]);
           }

           if (Method==2){
    	 A_d[0] =  2.0*A_sd[0];
    	 ptr->D[nPoints+0+j*nPoints] = 3.0 * qdy[0];
             A_d[nPoints-1] =  2.0*A_sd[nPoints-2];
             ptr->D[nPoints+nPoints-1+j*nPoints] =  3.0 * qdy[nPoints-2];
    	 Mytridiagldltsolve(A_d, A_sd, &ptr->D[nPoints+j*nPoints], nPoints);
           }

           if (Method==3){
    	 /*  s'''(x(2)-) = s'''(x(2)+) */
    	 r = A_sd[1]/A_sd[0];
    	 A_d[0]= A_sd[0]/(1.0 + r);
    	 ptr->D[nPoints+j*nPoints]=((3.0*r+2.0)*qdy[0]+r*qdy[1])/((1.0+r)*(1.0+r));
    	 /*  s'''(x(n-1)-) = s'''(x(n-1)+) */
    	 r = A_sd[nPoints-3]/A_sd[nPoints-2];
    	 A_d[nPoints-1] = A_sd[nPoints-2]/(1.0 + r);
    	 ptr->D[nPoints+nPoints-1+j*nPoints] = ((3.0*r+2.0)*qdy[nPoints-2]+r*qdy[nPoints-3])/((1.0+r)*(1.0+r));
    	 Mytridiagldltsolve(A_d, A_sd, &ptr->D[nPoints+j*nPoints], nPoints);
           }
         }
         /* ********* */
         scicos_free(spline);
       }

       //===================================
       cnt1=nPoints-1;
       cnt2=nPoints;
       for (i=0;i<nPoints;i++){ // finding the first positive time instant
         if (ptr->workt[i]>=0 ) {
           cnt1=i-1;
           cnt2=i;
           break;
         }
       }
       ptr->nPoints=nPoints;
       ptr->Yt=Ytype;
       ptr->Yst=YsubType;
       ptr->cnt1=cnt1;
       ptr->cnt2=cnt2;
       ptr->EVindex=0;
       ptr->PerEVcnt=0;
       ptr->firstevent=1;
       return;
       /*******************************************************/
       /*******************************************************/
     }else if (flag==1){   /* event date computation */
       ptr = *(block->work);
       nPoints=ptr->nPoints;
       t=get_scicos_time();
       cnt1=ptr->cnt1;
       cnt2=ptr->cnt2;
       EVindex= ptr->EVindex;
       PerEVcnt=ptr->PerEVcnt;
       t1=t;

       if (ZC==1){
         if (OutEnd==2) t-=(PerEVcnt)*TP;
         inow=nPoints-1;
         for (i=cnt1;i<nPoints;i++){
           if (i==-1) continue;
           if (t<ptr->workt[i]) {
    	 inow=i-1;
    	 if (inow<cnt2){
    	   cnt2=inow;
    	 }else{
    	   cnt1=cnt2;
    	   cnt2=inow;
    	 }
    	 break;
           }
         }
       }else{
         if (OutEnd==2){ if (TP!=0) r=floor((t/TP));else r=0;t-=((int)r)*TP;};
         inow=nPoints-1;
         for (i=0;i<nPoints;i++){
           if (t<ptr->workt[i]){inow=i-1; break;}
         }
       }
       ptr->cnt1=cnt1;
       ptr->cnt2=cnt2;
       ptr->EVindex=EVindex;
       ptr->PerEVcnt=PerEVcnt;
       for (j=0;j<my;j++){

         if (ptr->Yt==1){
           if ((ptr->Yst==0)||(ptr->Yst==1)){ /*  if Real or complex*/
    	 y_d = GetRealOutPortPtrs(block,1);
    	 ptr_d=(double*) ptr->work;
    	 ptr_D=(double*) ptr->D;

    	 if (inow>=nPoints-1){
    	   if (OutEnd==0){
    	     y_d[j]=0.0;// outputs set to zero
    	   }else if (OutEnd==1){
    	     y_d[j]=ptr_d[nPoints-1+(j)*nPoints]; // hold outputs at the end
    	   }
    	 }else if (Method==0){
    	   if (inow<0){
    	     y_d[j]=0.0;
    	   }else{
    	     y_d[j]=ptr_d[inow+(j)*nPoints];
    	   }
    	 }else if (Method==1){
    	   if (inow<0){inow=0;}
    	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=ptr_d[inow+j*nPoints];
    	   y2=ptr_d[inow+1+j*nPoints];
    	   y_d[j]=(y2-y1)*(t-t1)/(t2-t1)+y1;
    	 }else if (Method>=2){
    	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=ptr_d[inow+j*nPoints];
    	   y2=ptr_d[inow+1+j*nPoints];
    	   d1=ptr_D[inow+j*nPoints];
    	   d2=ptr_D[inow+1+j*nPoints];
    	   Myevalhermite2(&t, &t1,&t2, &y1,&y2, &d1,&d2, &h, &dh, &ddh, &dddh, &inow);
    	   y_d[j]=h;
    	 }
           }
           if (ptr->Yst==1){/*  --------------Complex----------------------*/
    	 y_cd = GetImagOutPortPtrs(block,1);
    	 if (inow>=nPoints-1){
    	   if (OutEnd==0){
    	     y_cd[j]=0.0;// outputs set to zero
    	   }else if (OutEnd==1){
    	     y_cd[j]=ptr_d[nPoints*my+nPoints-1+(j)*nPoints]; // hold outputs at the end
    	   }
    	 }else if (Method==0){
    	   if (inow<0){
    	     y_cd[j]=0.0;// outputs set to zero
    	   }else{
    	     y_cd[j]=ptr_d[nPoints*my+inow+(j)*nPoints];
    	   }
    	 }else if (Method==1){
    	   if (inow<0){inow=0;} // extrapolation for 0<t<X(0)
    	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=ptr_d[nPoints*my+inow+j*nPoints];
    	   y2=ptr_d[nPoints*my+inow+1+j*nPoints];
    	   y_cd[j]=(y2-y1)*(t-t1)/(t2-t1)+y1;
    	 }else if (Method>=2){
    	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=ptr_d[inow+j*nPoints+nPoints];
    	   y2=ptr_d[inow+1+j*nPoints+nPoints];
    	   d1=ptr_D[inow+j*nPoints+nPoints];
    	   d2=ptr_D[inow+1+j*nPoints+nPoints];
    	   Myevalhermite2(&t, &t1,&t2, &y1,&y2, &d1,&d2, &h, &dh, &ddh, &dddh, &inow);
    	   y_cd[j]=h;
    	 }
           }
         }else if (ptr->Yt==8){
           switch (ptr->Yst){
           case 1: // ---------------------int8 char  ----------------------------
    	 y_c = Getint8OutPortPtrs(block,1);
    	 ptr_c=(char*) ptr->work;
    	 //y_c[j]=ptr_c[inow+(j)*nPoints];
    	 if (inow>=nPoints-1){
    	   if (OutEnd==0){
    	     y_c[j]=0;// outputs set to zero
    	   }else if (OutEnd==1){
    	     y_c[j]=ptr_c[nPoints-1+(j)*nPoints]; // hold outputs at the end
    	   }
    	 }else if (Method==0){
    	   if (inow<0){
    	     y_c[j]=0;
    	   }else{
    	     y_c[j]=ptr_c[inow+(j)*nPoints];
    	   }
    	 }else if (Method>=1){
        	   if (inow<0){inow=0;}
    	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=(double)ptr_c[inow+j*nPoints];
    	   y2=(double)ptr_c[inow+1+j*nPoints];
    	   y_c[j] =(char)((y2-y1)*(t-t1)/(t2-t1)+y1);
    	 }
    	 break;
           case 2: // ---------------------int16 short---------------------
    	 y_s = Getint16OutPortPtrs(block,1);
    	 ptr_s=(short*) ptr->work;
    	 // y_s[j]=ptr_s[inow+(j)*nPoints];
    	 if (inow>=nPoints-1){
    	   if (OutEnd==0){
    	     y_s[j]=0;// outputs set to zero
    	   }else if (OutEnd==1){
    	     y_s[j]=ptr_s[nPoints-1+(j)*nPoints]; // hold outputs at the end
    	   }
    	 }else if (Method==0){
    	   if (inow<0){
    	     y_s[j]=0;
    	   }else{
    	     y_s[j]=ptr_s[inow+(j)*nPoints];
    	   }
    	 }else if (Method>=1){
    	   if (inow<0){inow=0;}
    	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=(double)ptr_s[inow+j*nPoints];
    	   y2=(double)ptr_s[inow+1+j*nPoints];
    	   y_s[j] =(short)((y2-y1)*(t-t1)/(t2-t1)+y1);
    	 }
    	 break;
           case 4: // ---------------------int32 long---------------------
    	 y_l = Getint32OutPortPtrs(block,1);
    	 ptr_l=(long*) ptr->work;
    	 //y_l[j]=ptr_l[inow+(j)*nPoints];
    	 if (inow>=nPoints-1){
    	   if (OutEnd==0){
    	     y_l[j]=0;// outputs set to zero
    	   }else if (OutEnd==1){
    	     y_l[j]=ptr_l[nPoints-1+(j)*nPoints]; // hold outputs at the end
    	   }
    	 }else if (Method==0){
    	   if (inow<0){
    	     y_l[j]=0;
    	   }else{
    	     y_l[j]=ptr_l[inow+(j)*nPoints];
    	   }
    	 }else if (Method>=1){
    	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=(double)ptr_l[inow+j*nPoints];
    	   y2=(double)ptr_l[inow+1+j*nPoints];
    	   y_l[j] =(long)((y2-y1)*(t-t1)/(t2-t1)+y1);

    	 }
    	 break;
           case 11: //--------------------- uint8 uchar---------------------
    	 y_uc = Getuint8OutPortPtrs(block,1);
    	 ptr_uc=(unsigned char*) ptr->work;
    	 //y_uc[j]=ptr_uc[inow+(j)*nPoints];
    	 if (inow>=nPoints-1){
    	   if (OutEnd==0){
    	     y_uc[j]=0;// outputs set to zero
    	   }else if (OutEnd==1){
    	     y_uc[j]=ptr_uc[nPoints-1+(j)*nPoints]; // hold outputs at the end
    	   }
    	 }else if (Method==0){
    	   if (inow<0){
    	     y_uc[j]=0;
    	   }else{
    	     y_uc[j]=ptr_uc[inow+(j)*nPoints];
    	   }
    	 }else if (Method>=1){
      	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=(double)ptr_uc[inow+j*nPoints];
    	   y2=(double)ptr_uc[inow+1+j*nPoints];
    	   y_uc[j] =(unsigned char)((y2-y1)*(t-t1)/(t2-t1)+y1);
    	 }
           break;
           case 12: // ---------------------uint16 ushort---------------------
    	 y_us = Getuint16OutPortPtrs(block,1);
    	 ptr_us=(unsigned short*) ptr->work;
    	 // y_us[j]=ptr_us[inow+(j)*nPoints];
    	 if (inow>=nPoints-1){
    	   if (OutEnd==0){
    	     y_us[j]=0;// outputs set to zero
    	   }else if (OutEnd==1){
    	     y_us[j]=ptr_us[nPoints-1+(j)*nPoints]; // hold outputs at the end
    	   }
    	 }else if (Method==0){
    	   if (inow<0){
    	     y_us[j]=0;
    	   }else{
    	     y_us[j]=ptr_us[inow+(j)*nPoints];
    	   }
    	 }else if (Method>=1){
           	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=(double)ptr_us[inow+j*nPoints];
    	   y2=(double)ptr_us[inow+1+j*nPoints];
    	   y_us[j] =(unsigned short)((y2-y1)*(t-t1)/(t2-t1)+y1);
    	 }
    	 break;
           case 14: // ---------------------uint32 ulong---------------------
    	 y_ul = Getuint32OutPortPtrs(block,1);
    	 ptr_ul=(unsigned long*) ptr->work;
    	 // y_ul[j]=ptr_ul[inow+(j)*nPoints];
    	 if (inow>=nPoints-1){
    	   if (OutEnd==0){
    	     y_ul[j]=0;// outputs set to zero
    	   }else if (OutEnd==1){
    	     y_ul[j]=ptr_ul[nPoints-1+(j)*nPoints]; // hold outputs at the end
    	   }
    	 }else if (Method==0){
    	   if (inow<0){
    	     y_ul[j]=0;
    	   }else{
    	     y_ul[j]=ptr_ul[inow+(j)*nPoints];
    	   }
    	 }else if (Method>=1){
    	   t1=ptr->workt[inow];
    	   t2=ptr->workt[inow+1];
    	   y1=(double)ptr_ul[inow+j*nPoints];
    	   y2=(double)ptr_ul[inow+1+j*nPoints];
    	   y_ul[j] =(unsigned long)((y2-y1)*(t-t1)/(t2-t1)+y1);
    	 }
    	 break;
           }
         }
       }// for j loop
       /********************************************************************/
     }else if(flag==3){   /* event date computation */
      t=get_scicos_time();
       ptr = *(block->work);
       nPoints=ptr->nPoints;
       cnt1=ptr->cnt1;
       cnt2=ptr->cnt2;
       EVindex= ptr->EVindex;
       PerEVcnt=ptr->PerEVcnt;
       if (ZC==1) {/* generate Events only if ZC is active */
         if ((Method==1)||(Method==0)){
           //-------------------------
           if (ptr->firstevent==1){
    	 jfirst=nPoints-1; // finding first positive time instant
    	 for (j=0;j<nPoints;j++){
    	   if (ptr->workt[j]>0) {
    	     jfirst=j;
    	     break;
    	   }
    	 }
    	 block->evout[0]=ptr->workt[jfirst];
    	 EVindex=jfirst;
    	 ptr->EVindex=EVindex;
    	 ptr->firstevent=0;
    	 return;
           }
           //------------------------
           i=EVindex;
           //------------------------
           if (i<nPoints-1) {
    	 block->evout[0]=ptr->workt[i+1]-ptr->workt[i];
    	 EVindex=i+1;
           }
           //------------------------
           if (i==nPoints-1){
    	 if (OutEnd==2) {/*  Periodic*/
    	   cnt1=-1;
    	   cnt2=0;
    	   PerEVcnt++;/* When OutEnd==2 (perodic output)*/
    	   jfirst=nPoints-1; // finding first positive time instant
    	   for (j=0;j<nPoints;j++){
    	     if (ptr->workt[j]>0) {
    	       jfirst=j;
    	       break;
    	     }
    	   }
    	   block->evout[0]=ptr->workt[jfirst];
    	   EVindex=jfirst;
    	 }
           }
           //--------------------------
         }else if (Method<=3){
           if (ptr->firstevent==1){
    	   block->evout[0]=TP;
    	   ptr->firstevent=0;
           }else{
    	 if (OutEnd==2){
    	   block->evout[0]=TP;
    	 }
    	 PerEVcnt++;
           }
           cnt1=-1;
           cnt2=0;
         }
         ptr->cnt1=cnt1;
         ptr->cnt2=cnt2;
         ptr->EVindex=EVindex;
         ptr->PerEVcnt=PerEVcnt;
       }
       /***********************************************************************/
     }else if (flag==5){   /* finish */
       ptr = *(block->work);
       if (ptr!=NULL) {
         if (ptr->D!=NULL) {
           scicos_free(ptr->D);
         }
         if (ptr->work!=NULL) {
           scicos_free(ptr->work);
         }
         if (ptr->workt!=NULL) {
           scicos_free(ptr->workt);
         }
         scicos_free(ptr);
       }
     }
     /*************************************************************************/
    }
    int Mytridiagldltsolve(double *d__, double * l, double * b, int n)
    {
        int i__1;
        double temp;
        int i__;
        --b;
        --l;
        --d__;

        i__1 = n;
        for (i__ = 2; i__ <= i__1; ++i__) {
    	temp = l[i__ - 1];
    	l[i__ - 1] /= d__[i__ - 1];
    	d__[i__] -= temp * l[i__ - 1];
    	b[i__] -= l[i__ - 1] * b[i__ - 1];
        }
        b[n] /= d__[n];
        for (i__ = n - 1; i__ >= 1; --i__) {
    	b[i__] = b[i__] / d__[i__] - l[i__] * b[i__ + 1];
        }
        return 0;
    }

    int Myevalhermite2(double *t, double *xa, double *xb, double *ya, double *yb, double *da, double *db, double *h, double *dh, double *ddh, double *dddh, int *i)
    {
      double tmxa, p, c2, c3, dx;

      /*    if (old_i != *i) {*/
    /*        compute the following Newton form : */
    /*           h(t) = ya + da*(t-xa) + c2*(t-xa)^2 + c3*(t-xa)^2*(t-xb) */
    	dx = 1. / (*xb - *xa);
    	p = (*yb - *ya) * dx;
    	c2 = (p - *da) * dx;
    	c3 = (*db - p + (*da - p)) * (dx * dx);
    	/*	}	 old_i = *i;*/

    /*     eval h(t), h'(t), h"(t) and h"'(t), by a generalised Horner 's scheme */
        tmxa = *t - *xa;
        *h = c2 + c3 * (*t - *xb);
        *dh = *h + c3 * tmxa;
        *ddh = (*dh + c3 * tmxa) * 2.;
        *dddh = c3 * 6.;
        *h = *da + *h * tmxa;
        *dh = *h + *dh * tmxa;
        *h = *ya + *h * tmxa;
        return 0;
    } /* evalhermite_ */