#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_ */