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