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