de1e0e8a792831efa2fe871f5239eb896d1f88ac
[scilab.git] / scilab / modules / scicos_blocks / src / c / fromws_c.c
1 /*  Scicos
2 *
3 *  Copyright (C) INRIA - METALAU Project <scicos@inria.fr>
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18 *
19 * See the file ./license.txt
20 */
21 #include "scicos_block4.h"
22 /*    Masoud Najafi, Alan Layec September 2007 */
23 /*    Copyright INRIA
24  *    Scicos block simulator
25  *    From workspace block
26  */
27 #include "stack-c.h"
28 #include <stdio.h>
29 #include <string.h>
30 #include "machine.h"
31 #include <math.h>
32
33 #if WIN32
34 #define NULL    0
35 #endif
36
37 #define Fnlength  block->ipar[0]
38 #define FName     block->ipar[1]
39 #define Method     block->ipar[1+Fnlength]
40 #define ZC        block->ipar[2+Fnlength]
41 #define OutEnd    block->ipar[3+Fnlength]
42 #define T0        ptr->workt[0]
43 #define TNm1      ptr->workt[nPoints-1]
44 #define TP        (TNm1-0)
45
46 extern int C2F(cvstr) __PARAMS((integer *,integer *,char *,integer *,unsigned long int));
47 extern int C2F(mgetnc)();
48 extern void C2F(mopen)();
49 extern int C2F(cluni0) __PARAMS((char *name, char *nams, integer *ln, long int name_len,
50                                 long int nams_len));
51 extern void C2F(mclose) __PARAMS((integer *fd, double *res));
52 extern void sciprint __PARAMS((char *fmt,...));
53 int Mytridiagldltsolve(double *d, double * l, double * b, int n);
54 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);
55 /*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);*/
56
57 /* function to check and extract data coming from an hypermat */
58 int Ishm(int *fd,int *Ytype,int *nPoints,int *my,int *ny,int *YsubType);
59
60 static char fmtd[3]={'d','l','\000'};
61 static char fmti[3]={'i','l','\000'};
62 static char fmtl[3]={'l','l','\000'};
63 static char fmts[3]={'s','l','\000'}; 
64 static char fmtc[3]={'c','l','\000'};
65 static char fmtul[3]={'u','l','\000'};
66 static char fmtus[3]={'u','s','\000'};
67 static char fmtuc[3]={'u','c','\000'};
68
69 #ifdef hppa
70 #undef FILENAME_MAX
71 #define FILENAME_MAX 4096
72 #endif
73 /* work struct for that block */
74 typedef struct {
75   int nPoints;
76   int Hmat;
77   int Yt;
78   int Yst;
79   int cnt1;
80   int cnt2;
81   int EVindex;
82   int PerEVcnt;
83   int firstevent;
84   double *D;
85   void *work;
86   double *workt;
87 } fromwork_struct ;
88
89
90 void fromws_c(scicos_block *block,int flag)
91 {
92   double t,y1,y2,t1,t2,r;
93   double *spline, *A_d, *A_sd, *qdy;
94   /* double  a,b,c,*y;*/
95   double d1,d2,h, dh, ddh, dddh;
96   /* counter and indexes variables */
97   int i,inow;
98   int j,jfirst;
99   int cnt1, cnt2, EVindex, PerEVcnt;
100
101   /* variables to handle files of TMPDIR/Workspace */
102   int fd;
103   char *status;
104   int swap = 1;
105   double res;
106   int out_n;
107   long int lout;
108   char filename[FILENAME_MAX];
109   char str[100];
110   int ierr;
111
112   /* variables for type and dims of data coming from scilab */
113   int Ytype, YsubType, mY, nY;
114   int nPoints;
115   int Ydim[10];
116
117   /* variables for type and dims of data of the output port block */
118   int ytype, my, ny;
119
120   /* generic pointer */
121   SCSREAL_COP *y_d,*y_cd,*ptr_d, *ptr_T, *ptr_D;
122   SCSINT8_COP *y_c,*ptr_c;
123   SCSUINT8_COP *y_uc, *ptr_uc;
124   SCSINT16_COP *y_s,*ptr_s;
125   SCSUINT16_COP *y_us,*ptr_us;
126   SCSINT32_COP *y_l,*ptr_l;
127   SCSUINT32_COP *y_ul,*ptr_ul;
128
129   /* the struct ptr of that block */
130   fromwork_struct *ptr;
131
132   /* for path of TMPDIR/workspace */
133   char env[256];
134   char sep[2];
135 #ifdef _MSC_VER
136   sep[0]='\\';
137 #else
138   sep[0]='/';
139 #endif
140   sep[1]='\0';
141
142  /*retrieve dimensions of output port*/
143  my       = GetOutPortRows(block,1); /* number of rows of Outputs*/
144  ny       = GetOutPortCols(block,1); /* number of cols of Outputs*/
145  ytype    = GetOutType(block,1);     /* output type */
146
147  /* init */
148  if (flag==4) {
149    /* convert scilab code of the variable name to C string */
150    C2F(cvstr)(&(Fnlength),&(FName),str,(j=1,&j),(unsigned long)strlen(str));
151    str[Fnlength] = '\0';
152
153    /* retrieve path of TMPDIR/workspace */
154    strcpy(env,getenv("TMPDIR"));
155    strcat(env,sep);
156    strcat(env,"Workspace");
157    strcat(env,sep);
158    strcat(env,str);
159
160    /* open tmp file */
161    status = "rb"; /* "r" : read */
162                   /* "b" : binary format (required for Windows) */
163
164    lout=FILENAME_MAX;
165    C2F(cluni0)(env, filename, &out_n,1,lout);
166    C2F(mopen)(&fd,env,status,&swap,&res,&ierr);
167
168    if (ierr!=0) {
169      sciprint("The '%s' variable does not exist.\n",str);
170      set_block_error(-3);
171      return;
172    }
173
174    /* read x */
175    C2F(mgetnc) (&fd, &Ydim[0], (j=nsiz,&j), fmti, &ierr);  /* read sci id */
176    C2F(mgetnc) (&fd, &Ydim[6], (j=1,&j), fmti, &ierr);     /* read sci type */
177    if (Ydim[6]==17) {
178      if (!Ishm(&fd,&Ytype,&nPoints,&mY,&nY,&YsubType)) {
179        sciprint("Invalid variable type.\n");
180        set_block_error(-3); 
181        C2F(mclose)(&fd,&res);
182        return;
183      }
184      if (!((Ytype==1) || (Ytype==8))) {
185        sciprint("Invalid variable type.\n");
186        set_block_error(-3);
187        C2F(mclose)(&fd,&res);
188        return;
189      }
190    }
191    else if ((Ydim[6]==1)||(Ydim[6]==8)) {
192      C2F(mgetnc) (&fd, &Ydim[7], (j=3,&j), fmti, &ierr); /* read sci header */
193      Ytype    = Ydim[6]; /* data type        */
194      nPoints  = Ydim[7]; /* number of data   */
195      mY       = Ydim[8]; /* first dimension  */
196      nY       = 1;       /* second dimension */
197      YsubType = Ydim[9]; /* subtype          */
198    }
199    else {
200     sciprint("Invalid variable type.\n");
201     set_block_error(-3);
202     C2F(mclose)(&fd,&res);
203     return;
204    }
205
206    /* check dimension for output port and variable */
207    if ((mY!=my)||(nY!=ny)) {
208      sciprint("Data dimensions are inconsistent:\n\r Variable size=[%d,%d] \n\r"
209               "Block output size=[%d,%d].\n",mY,nY,my,ny);
210      set_block_error(-3);
211      C2F(mclose)(&fd,&res);
212      return;
213    }
214
215    /* check variable data type and output block data type */
216    if (Ytype==1) { /*real/complex cases*/
217      switch (YsubType)
218      {
219      case 0: if (ytype!=10) {
220                sciprint("Output should be of Real type.\n");
221                set_block_error(-3);
222                C2F(mclose)(&fd,&res);
223                return;
224              }
225              break;
226
227      case 1: if (ytype!=11) {
228                sciprint("Output should be of complex type.\n");
229                set_block_error(-3);
230                C2F(mclose)(&fd,&res);
231                return;
232              }
233              break;
234      }
235    }
236    else if(Ytype==8) { /*integer cases*/
237      switch (YsubType)
238      {
239      case 1: if (ytype!=81) {
240                sciprint("Output should be of int8 type.\n");
241                set_block_error(-3);
242                C2F(mclose)(&fd,&res);
243                return;
244              }
245              break;
246
247      case 2: if (ytype!=82) {
248                sciprint("Output should be of int16 type.\n");
249                set_block_error(-3);
250                C2F(mclose)(&fd,&res);
251                return;
252              }
253              break;
254
255      case 4: if (ytype!=84) {
256                sciprint("Output should be of int32 type.\n");
257                set_block_error(-3);
258                C2F(mclose)(&fd,&res);
259                return;
260              }
261              break;
262
263      case 11:if (ytype!=811) {
264                sciprint("Output should be of uint8 type.\n");
265                set_block_error(-3);
266                C2F(mclose)(&fd,&res);
267                return;
268              }
269              break;
270
271      case 12:if (ytype!=812) {
272                sciprint("Output should be of uint16 type.\n");
273                set_block_error(-3);
274                C2F(mclose)(&fd,&res);
275                return;
276              }
277              break;
278
279      case 14:if (ytype!=814) {
280                sciprint("Output should be of uint32 type.\n");
281                set_block_error(-3);
282                C2F(mclose)(&fd,&res);
283                return;
284              }
285              break;
286      }
287    }
288
289    /* allocation of the work structure of that block */
290    if((*(block->work)=(fromwork_struct*) scicos_malloc(sizeof(fromwork_struct)))==NULL) {
291      set_block_error(-16);
292      C2F(mclose)(&fd,&res);
293      return;
294    }
295    ptr = *(block->work);
296    ptr->D=NULL;
297    ptr->workt=NULL;
298    ptr->work=NULL;
299
300    if (Ytype==1) { /*real/complex case*/
301      switch (YsubType) {
302      case 0 : /* Real */
303        if((ptr->work=(void *) scicos_malloc(nPoints*mY*nY*sizeof(double)))==NULL) {
304          set_block_error(-16);
305          scicos_free(ptr);
306          *(block->work)=NULL;
307          C2F(mclose)(&fd,&res);
308          return;
309        }
310        ptr_d = (SCSREAL_COP *) ptr->work;
311        C2F(mgetnc) (&fd, ptr_d, (j=nPoints*mY*nY,&j), fmtd, &ierr);  /* read double data */
312        break;
313      case 1:  /* complex */
314        if((ptr->work=(void *) scicos_malloc(2*nPoints*mY*nY*sizeof(double)))==NULL) {
315          set_block_error(-16);
316          scicos_free(ptr);
317          *(block->work)=NULL;
318          C2F(mclose)(&fd,&res);
319          return;
320        }
321        ptr_d = (SCSREAL_COP *) ptr->work;
322        C2F(mgetnc) (&fd, ptr_d, (j=2*nPoints*mY*nY,&j), fmtd, &ierr);  /* read double data */
323        break;
324      }
325    }
326    else if(Ytype==8) { /*integer case*/
327      switch (YsubType) {
328      case 1 :/* int8 */
329        if((ptr->work=(void *) scicos_malloc(nPoints*mY*nY*sizeof(char)))==NULL) {
330          set_block_error(-16);
331          scicos_free(ptr);
332          *(block->work)=NULL;
333          C2F(mclose)(&fd,&res);
334          return;
335        }
336        ptr_c = (SCSINT8_COP *) ptr->work;
337        C2F(mgetnc) (&fd, ptr_c, (j=nPoints*mY*nY,&j), fmtc, &ierr);  /* read char data */
338        break;
339      case 2 :  /* int16 */
340        if((ptr->work=(void *) scicos_malloc(nPoints*mY*nY*sizeof(short)))==NULL) {
341          set_block_error(-16);
342          scicos_free(ptr);
343          *(block->work)=NULL;
344          C2F(mclose)(&fd,&res);
345          return;
346        }
347        ptr_s = (SCSINT16_COP *) ptr->work;
348        C2F(mgetnc) (&fd, ptr_s, (j=nPoints*mY*nY,&j), fmts, &ierr);  /* read short data */
349        break;
350      case 4 :   /* int32 */
351        if((ptr->work=(void *) scicos_malloc(nPoints*mY*nY*sizeof(long)))==NULL) {
352          set_block_error(-16);
353          scicos_free(ptr);
354          *(block->work)=NULL;
355          C2F(mclose)(&fd,&res);
356          return;
357        }
358        ptr_l = (SCSINT32_COP *) ptr->work;
359        C2F(mgetnc) (&fd, ptr_l, (j=nPoints*mY*nY,&j), fmtl, &ierr);  /* read short data */
360        break;
361      case 11 :   /* uint8 */
362        if((ptr->work=(void *) scicos_malloc(nPoints*mY*nY*sizeof(unsigned char)))==NULL) {
363          set_block_error(-16);
364          scicos_free(ptr);
365          *(block->work)=NULL;
366          C2F(mclose)(&fd,&res);
367          return;
368        }
369        ptr_uc = (SCSUINT8_COP *) ptr->work;
370        C2F(mgetnc) (&fd, ptr_uc, (j=nPoints*mY*nY,&j), fmtuc, &ierr);  /* read short data */
371        break;
372      case 12 : /* uint16 */
373        if((ptr->work=(void *) scicos_malloc(nPoints*mY*nY*sizeof(unsigned short)))==NULL) {
374          set_block_error(-16);
375          scicos_free(ptr);
376          *(block->work)=NULL;
377          C2F(mclose)(&fd,&res);
378          return;
379        }
380        ptr_us = (SCSUINT16_COP *) ptr->work;
381        C2F(mgetnc) (&fd, ptr_us, (j=nPoints*mY*nY,&j), fmtus, &ierr);  /* read short data */
382        break;
383      case 14 :  /* uint32 */
384        if((ptr->work=(void *) scicos_malloc(nPoints*mY*nY*sizeof(unsigned long)))==NULL) {
385          set_block_error(-16);
386          scicos_free(ptr);
387          *(block->work)=NULL;
388          C2F(mclose)(&fd,&res);
389          return;
390        }
391        ptr_ul = (SCSUINT32_COP *) ptr->work;
392        C2F(mgetnc) (&fd, ptr_ul, (j=nPoints*mY*nY,&j), fmtul, &ierr);  /* read short data */
393        break;
394      }
395    }
396
397    /* check Hmat */
398    if (Ydim[6]==17) {
399      ptr->Hmat=1;
400    }
401    else {
402      ptr->Hmat=0;
403    }
404
405    /* read t */
406    C2F(mgetnc) (&fd, &Ydim[0], (j=nsiz,&j), fmti, &ierr);  /* read sci id */
407    C2F(mgetnc) (&fd, &Ydim[6], (j=1,&j), fmti, &ierr);  /* read sci type */
408    C2F(mgetnc) (&fd, &Ydim[7], (j=3,&j), fmti, &ierr);  /* read sci header */
409
410    if (nPoints!=Ydim[7]) {
411      sciprint("The size of the Time(%d) and Data(%d) vectors are inconsistent.\n",Ydim[7],nPoints);
412      set_block_error(-3);
413      *(block->work)=NULL;
414      scicos_free(ptr->work);
415      scicos_free(ptr);
416      C2F(mclose)(&fd,&res);
417      return;
418    }
419
420    if ((Ydim[6]!=1) | (Ydim[9]!=0)) {
421      sciprint("The Time vector type is not ""double"".\n"); 
422      set_block_error(-3);
423      *(block->work)=NULL;
424      scicos_free(ptr->work);
425      scicos_free(ptr);
426      C2F(mclose)(&fd,&res);
427      return;
428    }
429
430    if((ptr->workt=(double *) scicos_malloc(nPoints*sizeof(double)))==NULL) {
431      set_block_error(-16);
432      *(block->work)=NULL;
433      scicos_free(ptr->work);
434      scicos_free(ptr);
435      C2F(mclose)(&fd,&res);
436      return;
437    }
438    ptr_T = (SCSREAL_COP *) ptr->workt;
439    C2F(mgetnc) (&fd, ptr_T, (j=nPoints,&j), fmtd, &ierr);  /* read data of t */
440
441    /* close the file*/
442    C2F(mclose)(&fd,&res);
443
444    /*================================*/
445    /* check for an increasing time data */
446    for(j = 0; j < nPoints-1; j++) {
447      if(ptr_T[j] > ptr_T[j+1]) {
448        sciprint("The time vector should be an increasing vector.\n");
449        set_block_error(-3);
450        *(block->work)=NULL;
451        scicos_free(ptr->workt);
452        scicos_free(ptr->work);
453        scicos_free(ptr);
454        return;
455      }
456    }
457    /*=================================*/
458    if ((Method>1)&&(Ytype==1)&&(!ptr->Hmat)) { /* double or complex */
459
460      if (YsubType==0) { /*real*/
461        if((ptr->D=(double *) scicos_malloc(nPoints*mY*sizeof(double)))==NULL) {
462          set_block_error(-16);
463          *(block->work)=NULL;
464          scicos_free(ptr->workt);
465          scicos_free(ptr->work);
466          scicos_free(ptr);
467          return;
468        }
469      }
470      else { /*complex*/
471        if((ptr->D=(double *) scicos_malloc(2*nPoints*mY*sizeof(double)))==NULL) {
472          set_block_error(-16);
473          *(block->work)=NULL;
474          scicos_free(ptr->workt);
475          scicos_free(ptr->work);
476          scicos_free(ptr);
477          return;
478        }
479      }
480
481      if((spline=(double *) scicos_malloc((3*nPoints-2)*sizeof(double)))==NULL) {
482        sciprint("Allocation problem in spline.\n");
483        set_block_error(-16);
484        *(block->work)=NULL;
485        scicos_free(ptr->D);
486        scicos_free(ptr->workt);
487        scicos_free(ptr->work);
488        scicos_free(ptr);
489        return;
490      }
491
492      A_d  = spline;
493      A_sd = A_d  + nPoints;
494      qdy  = A_sd + nPoints-1;
495
496      for (j=0;j<mY;j++) { /* real part */
497        for (i=0;i<=nPoints-2;i++) {
498          A_sd[i] = 1.0 / (ptr_T[i+1] - ptr_T[i]);
499          qdy[i]  = (ptr_d[i+1+j*nPoints] - ptr_d[i+j*nPoints]) * A_sd[i]*A_sd[i];
500        }
501
502        for (i=1;i<=nPoints-2;i++) {
503          A_d[i] = 2.0*(A_sd[i-1] +A_sd[i]);
504          ptr->D[i+j*nPoints] = 3.0*(qdy[i-1]+qdy[i]);
505        }
506
507        if (Method==2) {
508          A_d[0] =  2.0*A_sd[0];
509          ptr->D[0+j*nPoints] = 3.0 * qdy[0];
510          A_d[nPoints-1] =  2.0*A_sd[nPoints-2];
511          ptr->D[nPoints-1+j*nPoints] =  3.0 * qdy[nPoints-2];
512          Mytridiagldltsolve(A_d, A_sd, &ptr->D[j*nPoints], nPoints);
513        }
514
515        if (Method==3) {
516          /*  s'''(x(2)-) = s'''(x(2)+) */
517          r = A_sd[1]/A_sd[0];
518          A_d[0]= A_sd[0]/(1.0 + r);
519          ptr->D[j*nPoints]=((3.0*r+2.0)*qdy[0]+r*qdy[1])/((1.0+r)*(1.0+r));
520          /*  s'''(x(n-1)-) = s'''(x(n-1)+) */
521          r = A_sd[nPoints-3]/A_sd[nPoints-2];
522          A_d[nPoints-1] = A_sd[nPoints-2]/(1.0 + r);
523          ptr->D[nPoints-1+j*nPoints] = \
524                   ((3.0*r+2.0)*qdy[nPoints-2]+r*qdy[nPoints-3])/((1.0+r)*(1.0+r));
525          Mytridiagldltsolve(A_d, A_sd, &ptr->D[j*nPoints], nPoints);
526        }
527      }
528
529      if (YsubType==1) { /* imag part */
530        for (j=0;j<mY;j++) {
531          for (i=0;i<=nPoints-2;i++) {
532            A_sd[i] = 1.0 / (ptr_T[i+1] - ptr_T[i]);
533            qdy[i]  = (ptr_d[nPoints+i+1+j*nPoints] - ptr_d[nPoints+i+j*nPoints]) * A_sd[i]*A_sd[i];
534          }
535
536          for (i=1;i<=nPoints-2;i++) {
537            A_d[i] = 2.0*(A_sd[i-1] +A_sd[i]);
538            ptr->D[i+j*nPoints+nPoints] = 3.0*(qdy[i-1]+qdy[i]);
539          }
540
541          if (Method==2) {
542            A_d[0] =  2.0*A_sd[0];
543            ptr->D[nPoints+0+j*nPoints] = 3.0 * qdy[0];
544            A_d[nPoints-1] =  2.0*A_sd[nPoints-2];
545            ptr->D[nPoints+nPoints-1+j*nPoints] =  3.0 * qdy[nPoints-2];
546            Mytridiagldltsolve(A_d, A_sd, &ptr->D[nPoints+j*nPoints], nPoints);
547          }
548
549          if (Method==3) {
550            /*  s'''(x(2)-) = s'''(x(2)+) */
551            r = A_sd[1]/A_sd[0];
552            A_d[0]= A_sd[0]/(1.0 + r);
553            ptr->D[nPoints+j*nPoints]=((3.0*r+2.0)*qdy[0]+r*qdy[1])/((1.0+r)*(1.0+r));
554            /*  s'''(x(n-1)-) = s'''(x(n-1)+) */
555            r = A_sd[nPoints-3]/A_sd[nPoints-2];
556            A_d[nPoints-1] = A_sd[nPoints-2]/(1.0 + r);
557            ptr->D[nPoints+nPoints-1+j*nPoints] = \
558                     ((3.0*r+2.0)*qdy[nPoints-2]+r*qdy[nPoints-3])/((1.0+r)*(1.0+r));
559            Mytridiagldltsolve(A_d, A_sd, &ptr->D[nPoints+j*nPoints], nPoints);
560          }
561        }
562      }
563
564      scicos_free(spline);
565    }
566    /*===================================*/
567    cnt1=nPoints-1;
568    cnt2=nPoints;
569    for (i=0;i<nPoints;i++) { /* finding the first positive time instant */
570      if (ptr->workt[i]>=0 ) {
571        cnt1=i-1;
572        cnt2=i;
573        break;
574      }
575    }
576    ptr->nPoints=nPoints;
577    ptr->Yt=Ytype;
578    ptr->Yst=YsubType;
579    ptr->cnt1=cnt1;
580    ptr->cnt2=cnt2;
581    ptr->EVindex=0;
582    ptr->PerEVcnt=0;
583    ptr->firstevent=1;
584    return;
585    /*******************************************************/
586    /*******************************************************/
587  }
588  else if (flag==1){   /* output computation */
589
590    /* retrieve ptr of the structure of that block */
591    ptr = *(block->work);
592    nPoints=ptr->nPoints;
593    cnt1=ptr->cnt1;
594    cnt2=ptr->cnt2;
595    EVindex= ptr->EVindex;
596    PerEVcnt=ptr->PerEVcnt;
597
598    /* get current simulation time */
599    t=get_scicos_time();
600    t1=t;
601
602    if (ZC==1){ /*zero crossing enable*/
603      if (OutEnd==2) {
604        t-=(PerEVcnt)*TP;
605      }
606      inow=nPoints-1;
607      for (i=cnt1;i<nPoints;i++) {
608        if (i==-1) {
609          continue;
610        }
611        if (t<ptr->workt[i]) {
612          inow=i-1;
613          if (inow<cnt2) {
614            cnt2=inow;
615          }
616          else {
617           cnt1=cnt2;
618           cnt2=inow;
619          }
620          break;
621        }
622      }
623    }
624    else { /*zero crossing disable*/
625      if (OutEnd==2) {
626        if (TP!=0) {
627          r=floor((t/TP));
628        }
629        else {
630          r=0;
631        }
632        t-=((int)r)*TP;
633      }
634      inow=nPoints-1;
635      for (i=0;i<nPoints;i++) {
636        if (t<ptr->workt[i]) {
637          inow=i-1;
638          break;
639        }
640      }
641    }
642
643    ptr->cnt1=cnt1;
644    ptr->cnt2=cnt2;
645    ptr->EVindex=EVindex;
646    ptr->PerEVcnt=PerEVcnt;
647
648    /***************************/
649    /* hypermatrix case */
650    if (ptr->Hmat) {
651
652      for (j=0;j<my*ny;j++) {
653        if (ptr->Yt==1) {
654          if (ptr->Yst==0) { /* real case */
655            y_d = GetRealOutPortPtrs(block,1);
656            ptr_d=(double*) ptr->work;
657
658            if (inow>=nPoints-1) {
659              if (OutEnd==0){
660                y_d[j]=0.0; /* outputs set to zero */
661              }
662              else if (OutEnd==1) {
663                y_d[j]=ptr_d[(nPoints-1)*ny*my+j]; /* hold outputs at the end */
664              }
665            }
666            else {
667              if (inow<0) {
668                y_d[j]=0.0;
669              }
670              else {
671                y_d[j]=ptr_d[inow*ny*my+j];
672              }
673            }
674          }
675          else { /* complexe case */
676            y_d = GetRealOutPortPtrs(block,1);
677            y_cd = GetImagOutPortPtrs(block,1);
678            ptr_d=(double*) ptr->work;
679
680            if (inow>=nPoints-1) {
681              if (OutEnd==0) {
682                y_d[j]=0.0;  /* outputs set to zero */
683                y_cd[j]=0.0; /* outputs set to zero */
684              }
685              else if (OutEnd==1) {
686                y_d[j]=ptr_d[(nPoints-1)*ny*my+j]; /* hold outputs at the end */
687                y_cd[j]=ptr_d[nPoints*my*ny+(nPoints-1)*ny*my+j];    /* hold outputs at the end */
688              }
689            }
690            else {
691              if (inow<0) {
692               y_d[j]=0.0;  /* outputs set to zero */
693               y_cd[j]=0.0; /* outputs set to zero */
694              }
695              else {
696                y_d[j]=ptr_d[inow*ny*my+j];
697                y_cd[j]=ptr_d[nPoints*my*ny+inow*ny*my+j];
698              }
699            }
700          }
701        }
702       else if (ptr->Yt==8) {
703         switch (ptr->Yst) {
704         case 1: /* ---------------------int8 char  ---------------------------- */
705           y_c = Getint8OutPortPtrs(block,1);
706           ptr_c=(char*) ptr->work;
707           if (inow>=nPoints-1) {
708             if (OutEnd==0) {
709               y_c[j]=0; /* outputs set to zero */
710             }
711             else if (OutEnd==1) {
712               y_c[j]=ptr_c[(nPoints-1)*ny*my+j]; /* hold outputs at the end */
713             }
714           }
715           else {
716             if (inow<0) {
717               y_c[j]=0;
718             }
719             else {
720               y_c[j]=ptr_c[inow*ny*my+j];
721             }
722           }
723           break;
724
725         case 2: /* ---------------------int16 short--------------------- */
726           y_s = Getint16OutPortPtrs(block,1);
727           ptr_s=(short*) ptr->work;
728           if (inow>=nPoints-1) {
729             if (OutEnd==0) {
730               y_s[j]=0; /* outputs set to zero */
731             }
732             else if (OutEnd==1) {
733               y_s[j]=ptr_s[(nPoints-1)*ny*my+j]; /* hold outputs at the end */
734             }
735           }
736           else {
737             if (inow<0) {
738               y_s[j]=0;
739             }
740             else {
741               y_s[j]=ptr_s[inow*ny*my+j];
742             }
743           }
744           break;
745
746         case 4: /* ---------------------int32 long--------------------- */
747           y_l = Getint32OutPortPtrs(block,1);
748           ptr_l=(long*) ptr->work;
749           if (inow>=nPoints-1) {
750             if (OutEnd==0) {
751               y_l[j]=0;/* outputs set to zero */
752             }
753             else if (OutEnd==1) {
754               y_l[j]=ptr_l[(nPoints-1)*ny*my+j]; /* hold outputs at the end */
755             }
756           }
757           else {
758             if (inow<0) {
759               y_l[j]=0;
760             }
761             else {
762               y_l[j]=ptr_l[inow*ny*my+j];
763             }
764           }
765           break;
766
767         case 11: /*--------------------- uint8 uchar---------------------*/
768           y_uc = Getuint8OutPortPtrs(block,1);
769           ptr_uc=(unsigned char*) ptr->work;
770           if (inow>=nPoints-1) {
771             if (OutEnd==0) {
772               y_uc[j]=0;/* outputs set to zero */
773             }
774             else if (OutEnd==1) {
775               y_uc[j]=ptr_uc[(nPoints-1)*ny*my+j]; /* hold outputs at the end */
776             }
777           }
778           else {
779             if (inow<0) {
780               y_uc[j]=0;
781             }
782             else {
783               y_uc[j]=ptr_uc[inow*ny*my+j];
784             }
785           }
786         break;
787
788         case 12: /* ---------------------uint16 ushort--------------------- */
789           y_us = Getuint16OutPortPtrs(block,1);
790           ptr_us=(unsigned short*) ptr->work;
791           if (inow>=nPoints-1) {
792             if (OutEnd==0) {
793               y_us[j]=0;/* outputs set to zero */
794             }
795             else if (OutEnd==1) {
796               y_us[j]=ptr_us[(nPoints-1)*ny*my+j]; /* hold outputs at the end */
797             }
798           }
799           else {
800             if (inow<0) {
801               y_us[j]=0;
802             }
803             else {
804               y_us[j]=ptr_us[inow*ny*my+j];
805             }
806           }
807           break;
808
809         case 14: /* ---------------------uint32 ulong--------------------- */
810           y_ul = Getuint32OutPortPtrs(block,1);
811           ptr_ul=(unsigned long*) ptr->work;
812           if (inow>=nPoints-1) {
813             if (OutEnd==0) {
814               y_ul[j]=0;/* outputs set to zero */
815             }
816             else if (OutEnd==1) {
817               y_ul[j]=ptr_ul[(nPoints-1)*ny*my+j]; /* hold outputs at the end */
818             }
819           }
820           else {
821             if (inow<0) {
822               y_ul[j]=0;
823             }
824             else {
825               y_ul[j]=ptr_ul[inow*ny*my+j];
826             }
827           }
828           break;
829         }
830       }
831      } /* for j loop */
832    }
833    /****************************/
834    /* scalar of vectorial case */
835    else {
836      for (j=0;j<my;j++) {
837        if (ptr->Yt==1) {
838          if ((ptr->Yst==0)||(ptr->Yst==1)) { /*  if Real or complex*/
839            y_d = GetRealOutPortPtrs(block,1);
840            ptr_d=(double*) ptr->work;
841            ptr_D=(double*) ptr->D;
842
843            if (inow>=nPoints-1) {
844              if (OutEnd==0){
845                y_d[j]=0.0; /* outputs set to zero */
846              }
847              else if (OutEnd==1) {
848                y_d[j]=ptr_d[nPoints-1+(j)*nPoints]; /* hold outputs at the end */
849              }
850            }
851            else if (Method==0) {
852              if (inow<0) {
853                y_d[j]=0.0;
854              }
855              else {
856                y_d[j]=ptr_d[inow+(j)*nPoints];
857              }
858            }
859            else if (Method==1) {
860              if (inow<0) {
861                inow=0;
862              }
863              t1=ptr->workt[inow];
864              t2=ptr->workt[inow+1];
865              y1=ptr_d[inow+j*nPoints];
866              y2=ptr_d[inow+1+j*nPoints];
867              y_d[j]=(y2-y1)*(t-t1)/(t2-t1)+y1;
868            }
869            else if (Method>=2) {
870              t1=ptr->workt[inow];
871              t2=ptr->workt[inow+1];
872              y1=ptr_d[inow+j*nPoints];
873              y2=ptr_d[inow+1+j*nPoints];
874              d1=ptr_D[inow+j*nPoints];
875              d2=ptr_D[inow+1+j*nPoints];
876              Myevalhermite2(&t, &t1,&t2, &y1,&y2, &d1,&d2, &h, &dh, &ddh, &dddh, &inow);
877              y_d[j]=h;
878            }
879          }
880          if (ptr->Yst==1) { /*  --------------complex----------------------*/
881            y_cd = GetImagOutPortPtrs(block,1);
882            if (inow>=nPoints-1) {
883              if (OutEnd==0) {
884                y_cd[j]=0.0;/* outputs set to zero*/
885              }
886              else if (OutEnd==1) {
887                y_cd[j]=ptr_d[nPoints*my+nPoints-1+(j)*nPoints]; // hold outputs at the end
888              }
889            }
890            else if (Method==0){
891              if (inow<0){
892               y_cd[j]=0.0; /* outputs set to zero */
893              }
894              else {
895                y_cd[j]=ptr_d[nPoints*my+inow+(j)*nPoints];
896              }
897            }
898            else if (Method==1) {
899              if (inow<0) {
900                inow=0;
901              } /* extrapolation for 0<t<X(0) */
902              t1=ptr->workt[inow];
903              t2=ptr->workt[inow+1];
904              y1=ptr_d[nPoints*my+inow+j*nPoints];
905              y2=ptr_d[nPoints*my+inow+1+j*nPoints];
906              y_cd[j]=(y2-y1)*(t-t1)/(t2-t1)+y1;
907            }
908            else if (Method>=2) {
909              t1=ptr->workt[inow];
910              t2=ptr->workt[inow+1];
911              y1=ptr_d[inow+j*nPoints+nPoints];
912              y2=ptr_d[inow+1+j*nPoints+nPoints];
913              d1=ptr_D[inow+j*nPoints+nPoints];
914              d2=ptr_D[inow+1+j*nPoints+nPoints];
915              Myevalhermite2(&t, &t1,&t2, &y1,&y2, &d1,&d2, &h, &dh, &ddh, &dddh, &inow);
916              y_cd[j]=h;
917            }
918          }
919        }
920        else if (ptr->Yt==8) {
921          switch (ptr->Yst) {
922          case 1: /* ---------------------int8 char  ---------------------------- */
923            y_c = Getint8OutPortPtrs(block,1);
924            ptr_c=(char*) ptr->work;
925            /*y_c[j]=ptr_c[inow+(j)*nPoints];*/
926            if (inow>=nPoints-1) {
927              if (OutEnd==0) {
928                y_c[j]=0; /* outputs set to zero */
929              }
930              else if (OutEnd==1) {
931                y_c[j]=ptr_c[nPoints-1+(j)*nPoints]; /* hold outputs at the end */
932              }
933            }
934            else if (Method==0) {
935              if (inow<0) {
936                y_c[j]=0;
937              }
938              else {
939                y_c[j]=ptr_c[inow+(j)*nPoints];
940              }
941            }
942            else if (Method>=1){
943              if (inow<0) {
944                inow=0;
945              }
946              t1=ptr->workt[inow];
947              t2=ptr->workt[inow+1];
948              y1=(double)ptr_c[inow+j*nPoints];
949              y2=(double)ptr_c[inow+1+j*nPoints];
950              y_c[j] =(char)((y2-y1)*(t-t1)/(t2-t1)+y1);
951            }
952            break;
953          case 2: /* ---------------------int16 short--------------------- */
954            y_s = Getint16OutPortPtrs(block,1);
955            ptr_s=(short*) ptr->work;
956            /* y_s[j]=ptr_s[inow+(j)*nPoints]; */
957            if (inow>=nPoints-1) {
958              if (OutEnd==0) {
959                y_s[j]=0; /* outputs set to zero */
960              }
961              else if (OutEnd==1) {
962                 y_s[j]=ptr_s[nPoints-1+(j)*nPoints]; // hold outputs at the end
963              }
964            }
965            else if (Method==0) {
966              if (inow<0) {
967                y_s[j]=0;
968              }
969              else {
970                y_s[j]=ptr_s[inow+(j)*nPoints];
971              }
972            }
973            else if (Method>=1) {
974              if (inow<0) {
975                inow=0;
976              }
977              t1=ptr->workt[inow];
978              t2=ptr->workt[inow+1];
979              y1=(double)ptr_s[inow+j*nPoints];
980              y2=(double)ptr_s[inow+1+j*nPoints];
981              y_s[j] =(short)((y2-y1)*(t-t1)/(t2-t1)+y1);
982            }
983            break;
984          case 4: /* ---------------------int32 long--------------------- */
985            y_l = Getint32OutPortPtrs(block,1);
986            ptr_l=(long*) ptr->work;
987            /*y_l[j]=ptr_l[inow+(j)*nPoints];*/
988            if (inow>=nPoints-1) {
989              if (OutEnd==0) {
990                y_l[j]=0;/* outputs set to zero */
991              }
992              else if (OutEnd==1) {
993                y_l[j]=ptr_l[nPoints-1+(j)*nPoints]; /* hold outputs at the end */
994              }
995            }
996            else if (Method==0) {
997              if (inow<0) {
998                y_l[j]=0;
999              }
1000              else {
1001                y_l[j]=ptr_l[inow+(j)*nPoints];
1002              }
1003            }
1004            else if (Method>=1) {
1005              t1=ptr->workt[inow];
1006              t2=ptr->workt[inow+1];
1007              y1=(double)ptr_l[inow+j*nPoints];
1008              y2=(double)ptr_l[inow+1+j*nPoints];
1009              y_l[j] =(long)((y2-y1)*(t-t1)/(t2-t1)+y1);
1010            }
1011            break;
1012          case 11: /*--------------------- uint8 uchar---------------------*/
1013            y_uc = Getuint8OutPortPtrs(block,1);
1014            ptr_uc=(unsigned char*) ptr->work;
1015            /*y_uc[j]=ptr_uc[inow+(j)*nPoints];*/
1016            if (inow>=nPoints-1) {
1017              if (OutEnd==0) {
1018                y_uc[j]=0;/* outputs set to zero */
1019              }
1020              else if (OutEnd==1) {
1021                y_uc[j]=ptr_uc[nPoints-1+(j)*nPoints]; /* hold outputs at the end */
1022              }
1023            }
1024            else if (Method==0) {
1025              if (inow<0) {
1026                y_uc[j]=0;
1027              }
1028              else {
1029                y_uc[j]=ptr_uc[inow+(j)*nPoints];
1030              }
1031            }
1032            else if (Method>=1) {
1033              t1=ptr->workt[inow];
1034              t2=ptr->workt[inow+1];
1035              y1=(double)ptr_uc[inow+j*nPoints];
1036              y2=(double)ptr_uc[inow+1+j*nPoints];
1037              y_uc[j] =(unsigned char)((y2-y1)*(t-t1)/(t2-t1)+y1);
1038            }
1039          break;
1040          case 12: /* ---------------------uint16 ushort--------------------- */
1041            y_us = Getuint16OutPortPtrs(block,1);
1042            ptr_us=(unsigned short*) ptr->work;
1043            /* y_us[j]=ptr_us[inow+(j)*nPoints]; */
1044            if (inow>=nPoints-1) {
1045              if (OutEnd==0) {
1046                y_us[j]=0;/* outputs set to zero */
1047              }
1048              else if (OutEnd==1) {
1049                y_us[j]=ptr_us[nPoints-1+(j)*nPoints]; /* hold outputs at the end */
1050              }
1051            }
1052            else if (Method==0) {
1053              if (inow<0) {
1054                y_us[j]=0;
1055              }
1056              else {
1057                y_us[j]=ptr_us[inow+(j)*nPoints];
1058              }
1059            }
1060            else if (Method>=1) {
1061              t1=ptr->workt[inow];
1062              t2=ptr->workt[inow+1];
1063              y1=(double)ptr_us[inow+j*nPoints];
1064              y2=(double)ptr_us[inow+1+j*nPoints];
1065              y_us[j] =(unsigned short)((y2-y1)*(t-t1)/(t2-t1)+y1);
1066            }
1067            break;
1068          case 14: /* ---------------------uint32 ulong--------------------- */
1069            y_ul = Getuint32OutPortPtrs(block,1);
1070            ptr_ul=(unsigned long*) ptr->work;
1071            /* y_ul[j]=ptr_ul[inow+(j)*nPoints]; */
1072            if (inow>=nPoints-1) {
1073              if (OutEnd==0) {
1074                y_ul[j]=0;/* outputs set to zero */
1075              }
1076              else if (OutEnd==1) {
1077                y_ul[j]=ptr_ul[nPoints-1+(j)*nPoints]; /* hold outputs at the end */
1078              }
1079            }
1080            else if (Method==0) {
1081              if (inow<0) {
1082                y_ul[j]=0;
1083              }
1084              else {
1085                y_ul[j]=ptr_ul[inow+(j)*nPoints];
1086              }
1087            }
1088            else if (Method>=1) {
1089              t1=ptr->workt[inow];
1090              t2=ptr->workt[inow+1];
1091              y1=(double)ptr_ul[inow+j*nPoints];
1092              y2=(double)ptr_ul[inow+1+j*nPoints];
1093              y_ul[j] =(unsigned long)((y2-y1)*(t-t1)/(t2-t1)+y1);
1094            }
1095            break;
1096          }
1097        }
1098      } /* for j loop */
1099    }
1100    /********************************************************************/
1101  }
1102  else if(flag==3) {   /* event date computation */
1103    /* retrieve ptr of the structure of that block */
1104    ptr = *(block->work);
1105    nPoints=ptr->nPoints;
1106    cnt1=ptr->cnt1;
1107    cnt2=ptr->cnt2;
1108    EVindex= ptr->EVindex;
1109    PerEVcnt=ptr->PerEVcnt;
1110
1111    /* get current simulation time */
1112    t=get_scicos_time();
1113
1114    if (ZC==1) { /* generate Events only if ZC is active */
1115      if ((Method==1)||(Method==0)) {
1116        /*-------------------------*/
1117        if (ptr->firstevent==1) {
1118          jfirst=nPoints-1; /* finding first positive time instant */
1119          for (j=0;j<nPoints;j++) {
1120            if (ptr->workt[j]>0) {
1121              jfirst=j;
1122              break;
1123            }
1124          }
1125          block->evout[0]=ptr->workt[jfirst];
1126          EVindex=jfirst;
1127          ptr->EVindex=EVindex;
1128          ptr->firstevent=0;
1129          return;
1130        }
1131        /*------------------------*/
1132        i=EVindex;
1133        /*------------------------*/
1134        if (i<nPoints-1) {
1135          block->evout[0]=ptr->workt[i+1]-ptr->workt[i];
1136          EVindex=i+1;
1137        }
1138        /*------------------------*/
1139        if (i==nPoints-1) {
1140          if (OutEnd==2) {/*  Periodic*/
1141            cnt1=-1;
1142            cnt2=0;
1143            PerEVcnt++;/* When OutEnd==2 (perodic output)*/
1144            jfirst=nPoints-1; /* finding first positive time instant */
1145            for (j=0;j<nPoints;j++) {
1146              if (ptr->workt[j]>0) {
1147                jfirst=j;
1148                break;
1149              }
1150            }
1151            block->evout[0]=ptr->workt[jfirst];
1152            EVindex=jfirst;
1153          }
1154        }
1155        /*-------------------------- */
1156      }
1157      else if (Method<=3) {
1158        if (ptr->firstevent==1) {
1159          block->evout[0]=TP;
1160          ptr->firstevent=0;
1161        }
1162        else {
1163          if (OutEnd==2) {
1164            block->evout[0]=TP;
1165          }
1166          PerEVcnt++;
1167        }
1168        cnt1=-1;
1169        cnt2=0;
1170      }
1171      ptr->cnt1=cnt1;
1172      ptr->cnt2=cnt2;
1173      ptr->EVindex=EVindex;
1174      ptr->PerEVcnt=PerEVcnt;
1175    }
1176    /***********************************************************************/
1177  }
1178  else if (flag==5) { /* finish */
1179    ptr = *(block->work);
1180    if (ptr!=NULL) {
1181      if (ptr->D!=NULL) {
1182        scicos_free(ptr->D);
1183      }
1184      if (ptr->work!=NULL) {
1185        scicos_free(ptr->work);
1186      }
1187      if (ptr->workt!=NULL) {
1188        scicos_free(ptr->workt);
1189      }
1190      scicos_free(ptr);
1191    }
1192  }
1193  /*************************************************************************/
1194 }
1195
1196 int Ishm(int *fd,int *Ytype,int *nPoints,int *my,int *ny,int *YsubType)
1197 {
1198  int *ptr_i;
1199  int j,ierr;
1200
1201  /*work array to store header of hypermat*/
1202  if((ptr_i=(int *) scicos_malloc(37*sizeof(int)))==NULL) {
1203    return 0;
1204  }
1205
1206  C2F(mgetnc) (fd, ptr_i, (j=37,&j), fmti, &ierr);  /* read sci id */
1207  if (ierr!=0) {
1208    return 0;
1209  }
1210
1211  if ((ptr_i[0]!=3)  || \
1212      (ptr_i[1]!=1)  || \
1213      (ptr_i[5]!=10) || \
1214      (ptr_i[6]!=1)  || \
1215      (ptr_i[7]!=3)  || \
1216      (ptr_i[8]!=0)  || \
1217      (ptr_i[9]!=1)  || \
1218      (ptr_i[10]!=ptr_i[9]+2)  || \
1219      (ptr_i[11]!=ptr_i[10]+4) || \
1220      (ptr_i[12]!=ptr_i[11]+7) || \
1221      (ptr_i[13]!=17) || \
1222      (ptr_i[14]!=22) || \
1223      (ptr_i[15]!=13) || \
1224      (ptr_i[16]!=18) || \
1225      (ptr_i[17]!=22) || \
1226      (ptr_i[18]!=28) || \
1227      (ptr_i[19]!=14) || \
1228      (ptr_i[20]!=23) || \
1229      (ptr_i[21]!=29) || \
1230      (ptr_i[22]!=27) || \
1231      (ptr_i[23]!=18) || \
1232      (ptr_i[24]!=14) || \
1233      (ptr_i[25]!=28) || \
1234      (ptr_i[26]!=8)  || \
1235      (ptr_i[27]!=1)  || \
1236      (ptr_i[28]!=3)  || \
1237      (ptr_i[29]!=4))
1238   {
1239    sciprint("Error in hypermat scilab coding.\n");
1240    return 0;
1241   }
1242
1243  *my      = ptr_i[30];    /*37*/
1244  *ny      = ptr_i[31];    /*38*/
1245  *nPoints = ptr_i[32];    /*39*/
1246  *Ytype   = ptr_i[33];    /*40*/
1247
1248  if ((ptr_i[34]!=ptr_i[30]*ptr_i[31]*ptr_i[32]) || \
1249      (ptr_i[35]!=1))
1250   {
1251    sciprint("Error in hypermat scilab coding.\n");
1252    return 0;
1253   }
1254
1255  *YsubType = ptr_i[36];   /*43*/
1256
1257  scicos_free(ptr_i);
1258  return 1;
1259 }
1260
1261 int Mytridiagldltsolve(double *d, double * l, double * b, int n)
1262 {
1263
1264   /* variable declaration */
1265   double temp;
1266   int i;
1267
1268   --b;
1269   --l;
1270   --d;
1271
1272   for (i = 2; i <= n; ++i) {
1273     temp = l[i - 1];
1274     l[i - 1] /= d[i - 1];
1275     d[i] -= temp * l[i - 1];
1276     b[i] -= l[i - 1] * b[i - 1];
1277   }
1278   b[n] /= d[n];
1279
1280   for (i = n - 1; i >= 1; --i) {
1281     b[i] = b[i] / d[i] - l[i] * b[i + 1];
1282   }
1283
1284   return 0;
1285 }
1286
1287 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)
1288 {
1289   /* variable declaration */
1290   double tmxa, p, c2, c3, dx;
1291
1292   /* if (old_i != *i) {*/
1293   /*   compute the following Newton form : */
1294   /*    h(t) = ya + da*(t-xa) + c2*(t-xa)^2 + c3*(t-xa)^2*(t-xb) */
1295     dx = 1. / (*xb - *xa);
1296      p = (*yb - *ya) * dx;
1297     c2 = (p - *da) * dx;
1298     c3 = (*db - p + (*da - p)) * (dx * dx);
1299   /*        }         old_i = *i;*/
1300
1301   /*   eval h(t), h'(t), h"(t) and h"'(t), by a generalised Horner 's scheme */
1302   tmxa = *t - *xa;
1303   *h = c2 + c3 * (*t - *xb);
1304   *dh = *h + c3 * tmxa;
1305   *ddh = (*dh + c3 * tmxa) * 2.;
1306   *dddh = c3 * 6.;
1307   *h = *da + *h * tmxa;
1308   *dh = *h + *dh * tmxa;
1309   *h = *ya + *h * tmxa;
1310   return 0;
1311 } /* evalhermite_ */