Coverity #1099169, #1320913, #1321304 fixed
[scilab.git] / scilab / modules / scicos_blocks / src / cpp / fromws_c.cpp
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2015 - Scilab Enterprises - Paul Bignier
4  *  Copyright (C) INRIA - METALAU Project <scicos@inria.fr>
5  *
6  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14  */
15
16 #include <cmath>
17 #include <cstring>
18
19 #include "Helpers.hxx"
20
21 extern "C"
22 {
23 #include "api_scilab.h"
24 #include "dynlib_scicos_blocks.h"
25 #include "expandPathVariable.h"
26 #include "scicos_block4.h"
27 #include "scicos_evalhermite.h"
28 #include "scicos.h"
29 #include "sci_malloc.h"
30 #include "h5_fileManagement.h"
31 #include "h5_readDataFromFile.h"
32 #include "h5_attributeConstants.h"
33
34 #include "localization.h"
35
36     SCICOS_BLOCKS_IMPEXP void fromws_c(scicos_block* block, int flag);
37 }
38
39 /*--------------------------------------------------------------------------*/
40 /* work struct for a block */
41 typedef struct
42 {
43     int nPoints;
44     int Hmat;
45     int Yt;
46     int Yst;
47     int cnt1;
48     int cnt2;
49     int EVindex;
50     int PerEVcnt;
51     int firstevent;
52     double* D;
53     void* work;
54     double* workt;
55 } fromwork_struct;
56 /*--------------------------------------------------------------------------*/
57 static int Mytridiagldltsolve(double* &dA, double* &lA, double* &B, int N)
58 {
59     for (int j = 1; j <= N - 1; ++j)
60     {
61         double Temp = lA[j - 1];
62         lA[j - 1] /= dA[j - 1];
63         B[j] -= lA[j - 1] * B[j - 1];
64         dA[j] -= Temp * lA[j - 1];
65     }
66
67     B[N - 1] /= dA[N - 1];
68     for (int j = N - 2; j >= 0; --j)
69     {
70         B[j] = - lA[j] * B[j + 1] + B[j] / dA[j];
71     }
72
73     return 0;
74 }
75
76 /*--------------------------------------------------------------------------*/
77 using namespace org_scilab_modules_xcos_block;
78
79 /*--------------------------------------------------------------------------*/
80 SCICOS_BLOCKS_IMPEXP void fromws_c(scicos_block* block, int flag)
81 {
82     /* Retrieve dimensions of output port */
83     int my    = GetOutPortRows(block, 1); /* Number of rows of the output */
84     int ny    = GetOutPortCols(block, 1); /* Number of cols of the output */
85     int ytype = GetOutType(block, 1);     /* Output type */
86
87     /* Generic pointers */
88     double *y_d, *y_cd, *ptr_d = nullptr, *ptr_T, *ptr_D;
89     char *y_c, *ptr_c;
90     unsigned char *y_uc, *ptr_uc;
91     short int *y_s, *ptr_s;
92     unsigned short int *y_us, *ptr_us;
93     int *y_l, *ptr_l;
94     unsigned int *y_ul, *ptr_ul;
95
96     /* The struct pointer of the block */
97     fromwork_struct** work = (fromwork_struct**) block->work;
98     fromwork_struct* ptr = nullptr;
99
100     int Fnlength = block->ipar[0];
101     int Method   = block->ipar[1 + Fnlength];
102     int ZC       = block->ipar[2 + Fnlength];
103     int OutEnd   = block->ipar[3 + Fnlength];
104
105     switch (flag)
106     {
107         case 4 :
108         {
109             /* Init */
110
111             /* Convert Scilab code of the variable name to C string */
112
113             /* Path to "TMPDIR/Workspace/" */
114             const char* filePrefix = "TMPDIR/Workspace/";
115             const int prefixSize = static_cast<int>(strlen(filePrefix));
116
117             char FName[100];
118             for (int i = 0; i < Fnlength; ++i)
119             {
120                 FName[i] = static_cast<char>(block->ipar[1 + i]);
121             }
122             FName[Fnlength] = '\0';
123
124             char* env = new char[prefixSize + Fnlength + 1];
125             strcpy(env, filePrefix);
126             strcpy(env + prefixSize, FName);
127             env[prefixSize + Fnlength] = '\0';
128
129             char* filename = expandPathVariable(env);
130             delete[] env;
131             int fd = 0, ierr = 0;
132             if (filename)
133             {
134                 /* Open tmp file */
135                 fd = openHDF5File(filename, 0);
136             }
137             FREE(filename);
138             filename = nullptr;
139
140             if (fd < 0)
141             {
142                 Coserror(_("The '%s' variable does not exist.\n"), FName);
143                 return;
144             }
145
146             // Manage version information
147             int iVersion = getSODFormatAttribute(fd);
148             if (iVersion != SOD_FILE_VERSION)
149             {
150                 if (iVersion > SOD_FILE_VERSION)
151                 {
152                     // Cannot read file with version more recent that me!
153                     Coserror(_("%s: Wrong SOD file format version. Max Expected: %d Found: %d\n"), "fromws_c", SOD_FILE_VERSION, iVersion);
154                     return;
155                 }
156             }
157
158             /* Read the variables contained in the file */
159             char* pstNameList[2];
160             int nbVar = getVariableNames(fd, pstNameList);
161             if (nbVar != 2)
162             {
163                 Coserror(_("Erroneous saved file.\n"));
164                 return;
165             }
166             if (strcmp(pstNameList[0], "t") != 0 || strcmp(pstNameList[1], "x") != 0)
167             {
168                 Coserror(_("Erroneous saved file.\n"));
169                 return;
170             }
171
172             /* Read x */
173             int xSetId = getDataSetIdFromName(fd, pstNameList[1]);
174             int xType = getScilabTypeFromDataSet(xSetId);
175
176             int nPoints, mX, nX = 1;
177             int xDims, xSubType;
178             if ((xType == 1) || (xType == 8))
179             {
180                 int* pxDims = nullptr;
181                 getDatasetInfo(xSetId, &xSubType, &xDims, pxDims);
182                 pxDims = new int[xDims];
183                 getDatasetInfo(xSetId, &xSubType, &xDims, pxDims);
184                 nPoints = pxDims[0]; /* Number of data  */
185                 mX = pxDims[1];      /* First dimension  */
186                 if (xDims > 3)
187                 {
188                     nX = pxDims[2];  /* Second dimension */
189                 }
190                 delete[] pxDims;
191                 if (xType == 8)
192                 {
193                     getDatasetPrecision(xSetId, &xSubType); /* For the integer case, overwrite xSubType */
194                 }
195             }
196             else
197             {
198                 Coserror(_("Invalid variable type.\n"));
199                 /*scicos_print(_("Invalid variable type.\n"));
200                 set_block_error(-3);*/
201                 closeHDF5File(fd);
202                 return;
203             }
204
205             /* Check dimension for output port and variable */
206             if ((mX != my) || (nX != ny))
207             {
208                 Coserror(_("Data dimensions are inconsistent:\n Variable size=[%d,%d] \nBlock output size=[%d,%d].\n"), mX, nX, my, ny);
209                 /*set_block_error(-3);*/
210                 closeHDF5File(fd);
211                 return;
212             }
213
214             /* Check variable data type and output block data type */
215             if (xType == 1)
216             {
217                 /* real/complex cases */
218                 switch (xSubType)
219                 {
220                     case 0:
221                         if (ytype != 10)
222                         {
223                             Coserror(_("Output should be of Real type.\n"));
224                             /*set_block_error(-3);*/
225                             closeHDF5File(fd);
226                             return;
227                         }
228                         break;
229
230                     case 1:
231                         if (ytype != 11)
232                         {
233                             Coserror(_("Output should be of complex type.\n"));
234                             /*set_block_error(-3);*/
235                             closeHDF5File(fd);
236                             return;
237                         }
238                         break;
239                 }
240             }
241             else if (xType == 8)
242             {
243                 /* int cases */
244                 switch (xSubType)
245                 {
246                     case SCI_INT8:
247                         if (ytype != 81)
248                         {
249                             Coserror(_("Output should be of int8 type.\n"));
250                             set_block_error(-3);
251                             closeHDF5File(fd);
252                             return;
253                         }
254                         break;
255
256                     case SCI_INT16:
257                         if (ytype != 82)
258                         {
259                             Coserror(_("Output should be of int16 type.\n"));
260                             /*set_block_error(-3);*/
261                             closeHDF5File(fd);
262                             return;
263                         }
264                         break;
265
266                     case SCI_INT32:
267                         if (ytype != 84)
268                         {
269                             Coserror(_("Output should be of int32 type.\n"));
270                             /*set_block_error(-3);*/
271                             closeHDF5File(fd);
272                             return;
273                         }
274                         break;
275
276                     case SCI_UINT8:
277                         if (ytype != 811)
278                         {
279                             Coserror(_("Output should be of uint8 type.\n"));
280                             /*set_block_error(-3);*/
281                             closeHDF5File(fd);
282                             return;
283                         }
284                         break;
285
286                     case SCI_UINT16:
287                         if (ytype != 812)
288                         {
289                             Coserror(_("Output should be of uint16 type.\n"));
290                             /*set_block_error(-3);*/
291                             closeHDF5File(fd);
292                             return;
293                         }
294                         break;
295
296                     case SCI_UINT32:
297                         if (ytype != 814)
298                         {
299                             Coserror(_("Output should be of uint32 type.\n"));
300                             /*set_block_error(-3);*/
301                             closeHDF5File(fd);
302                             return;
303                         }
304                         break;
305                 }
306             }
307
308             /* Allocation of the work structure of that block */
309
310             *work = new fromwork_struct[sizeof(fromwork_struct)];
311             ptr = *work;
312             ptr->D = nullptr;
313             ptr->workt = nullptr;
314             ptr->work = nullptr;
315
316             if (xType == 1)
317             {
318                 /* real/complex case */
319                 switch (xSubType)
320                 {
321                     case 0 : /* Real */
322                         ptr->work = new double[nPoints * mX * nX];
323                         ptr_d = (double*) ptr->work;
324                         ierr = readDoubleMatrix(xSetId, ptr_d);
325                         break;
326                     case 1 :  /* Complex */
327                         ptr->work = new double[nPoints * mX * nX];
328                         ptr_d = (double*) ptr->work;
329                         ierr = readDoubleComplexMatrix(xSetId, ptr_d, ptr_d + nPoints * mX * nX);
330                         break;
331                 }
332             }
333             else if (xType == 8)
334             {
335                 /* int case */
336                 switch (xSubType)
337                 {
338                     case SCI_INT8 :
339                         ptr->work = new char[nPoints * mX * nX];
340                         ptr_c = (char*) ptr->work;
341                         ierr = readInteger8Matrix(xSetId, ptr_c);
342                         break;
343                     case SCI_INT16 :
344                         ptr->work = new short int[nPoints * mX * nX];
345                         ptr_s = (short int*) ptr->work;
346                         ierr = readInteger16Matrix(xSetId, ptr_s);
347                         break;
348                     case SCI_INT32 :
349                         ptr->work = new int[nPoints * mX * nX];
350                         ptr_l = (int*) ptr->work;
351                         ierr = readInteger32Matrix(xSetId, ptr_l);
352                         break;
353                     case SCI_UINT8 :
354                         ptr->work = new unsigned char[nPoints * mX * nX];
355                         ptr_uc = (unsigned char*) ptr->work;
356                         ierr = readUnsignedInteger8Matrix(xSetId, ptr_uc);
357                         break;
358                     case SCI_UINT16 :
359                         ptr->work = new unsigned short int[nPoints * mX * nX];
360                         ptr_us = (unsigned short int*) ptr->work;
361                         ierr = readUnsignedInteger16Matrix(xSetId, ptr_us);
362                         break;
363                     case SCI_UINT32 :
364                         ptr->work = new unsigned int[nPoints * mX * nX];
365                         ptr_ul = (unsigned int*) ptr->work;
366                         ierr = readUnsignedInteger32Matrix(xSetId, ptr_ul);
367                         break;
368                 }
369             }
370             if (ierr != 0)
371             {
372                 Coserror(_("Cannot read the values field.\n"));
373                 delete[] (char*) ptr->work;
374                 delete[] ptr;
375                 closeHDF5File(fd);
376             }
377
378             /* Check Hmat */
379             if (xDims > 2)
380             {
381                 ptr->Hmat = 1;
382             }
383             else
384             {
385                 ptr->Hmat = 0;
386             }
387
388             /* Read t */
389             int tSetId = getDataSetIdFromName(fd, pstNameList[0]);
390             int tType = getScilabTypeFromDataSet(tSetId);
391             if (tType != 1)
392             {
393                 Coserror(_("The Time vector type is not ""double"".\n"));
394                 set_block_error(-3);
395                 *work = nullptr;
396                 delete[] (char*) ptr->work;
397                 delete[] ptr;
398                 closeHDF5File(fd);
399                 return;
400             }
401
402             int tDims, tComplex, *ptDims = nullptr;
403             getDatasetInfo(tSetId, &tComplex, &tDims, ptDims);
404             ptDims = new int[tDims];
405             getDatasetInfo(tSetId, &tComplex, &tDims, ptDims);
406             if (tComplex != 0)
407             {
408                 Coserror(_("The Time vector type is complex.\n"));
409                 set_block_error(-3);
410                 *work = nullptr;
411                 delete[] (char*) ptr->work;
412                 delete[] ptr;
413                 delete[] ptDims;
414                 closeHDF5File(fd);
415                 return;
416             }
417             if (nPoints != ptDims[0])
418             {
419                 Coserror(_("The Time vector has a wrong size, expecting [%d, %d] and getting [%d, %d].\n"), nPoints, 1, ptDims[0], ptDims[1]);
420                 /*set_block_error(-3);*/
421                 *work = nullptr;
422                 delete[] (char*) ptr->work;
423                 delete[] ptr;
424                 delete[] ptDims;
425                 closeHDF5File(fd);
426                 return;
427             }
428             delete[] ptDims;
429
430             ptr->workt = new double[nPoints];
431             ptr_T = (double*) ptr->workt;
432             ierr = readDoubleMatrix(tSetId, ptr_T); /* Read t data */
433             if (ierr != 0)
434             {
435                 Coserror(_("Cannot read the time field.\n"));
436                 delete[] ptr->workt;
437                 delete[] (char*) ptr->work;
438                 delete[] ptr;
439                 closeHDF5File(fd);
440                 return;
441             }
442
443             /* Close the file */
444             closeHDF5File(fd);
445
446             /*================================*/
447             /* Check for an increasing time data */
448             for (int j = 0; j < nPoints - 1; ++j)
449             {
450                 if (ptr_T[j] > ptr_T[j + 1])
451                 {
452                     Coserror(_("The time vector should be an increasing vector.\n"));
453                     /*set_block_error(-3);*/
454                     *work = nullptr;
455                     delete[] ptr->workt;
456                     delete[] (char*) ptr->work;
457                     delete[] ptr;
458                     return;
459                 }
460             }
461             /*=================================*/
462             if ((Method > 1) && (xType == 1) && (!ptr->Hmat))
463             {
464                 /* double or complex */
465                 if (xSubType == 0) /* real */
466                 {
467                     ptr->D = new double[nPoints * mX];
468                 }
469                 else /* complex */
470                 {
471                     ptr->D = new double[2 * nPoints * mX];
472                 }
473
474                 double* spline = new double[3 * nPoints - 2];
475
476                 double* A_d  = spline;
477                 double* A_sd = A_d  + nPoints;
478                 double* qdy  = A_sd + nPoints - 1;
479
480                 for (int j = 0; j < mX; ++j)
481                 {
482                     /* real part */
483                     for (int i = 0; i <= nPoints - 2; ++i)
484                     {
485                         A_sd[i] = 1 / (ptr_T[i + 1] - ptr_T[i]);
486                         qdy[i]  = (ptr_d[i + 1 + j * nPoints] - ptr_d[i + j * nPoints]) * A_sd[i] * A_sd[i];
487                     }
488
489                     for (int i = 1; i <= nPoints - 2; ++i)
490                     {
491                         A_d[i] = 2 * (A_sd[i - 1] + A_sd[i]);
492                         ptr->D[i + j * nPoints] = 3 * (qdy[i - 1] + qdy[i]);
493                     }
494
495                     if (Method == 2)
496                     {
497                         A_d[0] =  2 * A_sd[0];
498                         ptr->D[0 + j * nPoints] = 3 * qdy[0];
499                         A_d[nPoints - 1] =  2 * A_sd[nPoints - 2];
500                         ptr->D[nPoints - 1 + j * nPoints] =  3 * qdy[nPoints - 2];
501                         double* res = &ptr->D[j * nPoints];
502                         Mytridiagldltsolve(A_d, A_sd, res, nPoints);
503                     }
504
505                     if (Method == 3)
506                     {
507                         /*  s'''(x(2)-) = s'''(x(2)+) */
508                         double r = A_sd[1] / A_sd[0];
509                         A_d[0] = A_sd[0] / (1 + r);
510                         ptr->D[j * nPoints] = ((3 * r + 2) * qdy[0] + r * qdy[1]) / ((1 + r) * (1 + r));
511                         /*  s'''(x(n-1)-) = s'''(x(n-1)+) */
512                         r = A_sd[nPoints - 3] / A_sd[nPoints - 2];
513                         A_d[nPoints - 1] = A_sd[nPoints - 2] / (1 + r);
514                         ptr->D[nPoints - 1 + j * nPoints] = ((3 * r + 2) * qdy[nPoints - 2] + r * qdy[nPoints - 3]) / ((1 + r) * (1 + r));
515                         double* res = &ptr->D[j * nPoints];
516                         Mytridiagldltsolve(A_d, A_sd, res, nPoints);
517                     }
518                 }
519
520                 if (xSubType == 1)
521                 {
522                     /* imag part */
523                     for (int j = 0; j < mX; ++j)
524                     {
525                         for (int i = 0; i <= nPoints - 2; ++i)
526                         {
527                             A_sd[i] = 1 / (ptr_T[i + 1] - ptr_T[i]);
528                             qdy[i]  = (ptr_d[nPoints + i + 1 + j * nPoints] - ptr_d[nPoints + i + j * nPoints]) * A_sd[i] * A_sd[i];
529                         }
530
531                         for (int i = 1; i <= nPoints - 2; ++i)
532                         {
533                             A_d[i] = 2 * (A_sd[i - 1] + A_sd[i]);
534                             ptr->D[i + j * nPoints + nPoints] = 3 * (qdy[i - 1] + qdy[i]);
535                         }
536
537                         if (Method == 2)
538                         {
539                             A_d[0] =  2 * A_sd[0];
540                             ptr->D[nPoints + 0 + j * nPoints] = 3 * qdy[0];
541                             A_d[nPoints - 1] =  2 * A_sd[nPoints - 2];
542                             ptr->D[nPoints + nPoints - 1 + j * nPoints] =  3 * qdy[nPoints - 2];
543                             double* res = &ptr->D[nPoints + j * nPoints];
544                             Mytridiagldltsolve(A_d, A_sd, res, nPoints);
545                         }
546
547                         if (Method == 3)
548                         {
549                             /*  s'''(x(2)-) = s'''(x(2)+) */
550                             double r = A_sd[1] / A_sd[0];
551                             A_d[0] = A_sd[0] / (1 + r);
552                             ptr->D[nPoints + j * nPoints] = ((3 * r + 2) * qdy[0] + r * qdy[1]) / ((1 + r) * (1 + r));
553                             /*  s'''(x(n-1)-) = s'''(x(n-1)+) */
554                             r = A_sd[nPoints - 3] / A_sd[nPoints - 2];
555                             A_d[nPoints - 1] = A_sd[nPoints - 2] / (1 + r);
556                             ptr->D[nPoints + nPoints - 1 + j * nPoints] = ((3 * r + 2) * qdy[nPoints - 2] + r * qdy[nPoints - 3]) / ((1 + r) * (1 + r));
557                             double* res = &ptr->D[nPoints + j * nPoints];
558                             Mytridiagldltsolve(A_d, A_sd, res, nPoints);
559                         }
560                     }
561                 }
562
563                 delete[] spline;
564             }
565             /*===================================*/
566             int cnt1 = nPoints - 1;
567             int cnt2 = nPoints;
568             for (int i = 0; i < nPoints; ++i)
569             {
570                 /* finding the first positive time instant */
571                 if (ptr->workt[i] >= 0)
572                 {
573                     cnt1 = i - 1;
574                     cnt2 = i;
575                     break;
576                 }
577             }
578             ptr->nPoints = nPoints;
579             ptr->Yt = xType;
580             ptr->Yst = xSubType;
581             ptr->cnt1 = cnt1;
582             ptr->cnt2 = cnt2;
583             ptr->EVindex = 0;
584             ptr->PerEVcnt = 0;
585             ptr->firstevent = 1;
586             break;
587             /*******************************************************/
588             /*******************************************************/
589         }
590         case 1 :
591         {
592             /* Output computation */
593
594             /* Retrieve 'ptr' of the structure of the block */
595             ptr = *work;
596             int nPoints = ptr->nPoints;
597             int cnt1 = ptr->cnt1;
598             int cnt2 = ptr->cnt2;
599             int EVindex = ptr->EVindex;
600             int PerEVcnt = ptr->PerEVcnt;
601
602             /* Get current simulation time */
603             double t = get_scicos_time();
604             double t1 = t, t2;
605
606             double TNm1  = ptr->workt[nPoints - 1];
607             double TP    = TNm1 - 0;
608
609             int inow;
610             if (ZC == 1)
611             {
612                 /* Zero-crossing enabled */
613                 if (OutEnd == 2)
614                 {
615                     if (PerEVcnt > 0)
616                     {
617                         // We ran out of value and OutEnd is 2 (Repeat)
618                         // Use fake time within our range.
619                         t -= (PerEVcnt) * TP;
620                     }
621                     inow = nPoints - 1;
622                 }
623                 else
624                 {
625                     inow = nPoints + 1; // Arbitrary value more than nPoints, will be overwritten if needed.
626                 }
627                 for (int i = cnt1 ; i < nPoints ; ++i)
628                 {
629                     if (i == -1)
630                     {
631                         continue;
632                     }
633                     if (t <= ptr->workt[i])
634                     {
635                         if (t < ptr->workt[i])
636                         {
637                             inow = i - 1;
638                         }
639                         else
640                         {
641                             inow = i;
642                         }
643                         if (inow < cnt2)
644                         {
645                             cnt2 = inow;
646                         }
647                         else
648                         {
649                             cnt1 = cnt2;
650                             cnt2 = inow;
651                         }
652                         break;
653                     }
654                 }
655             }
656             else   /* Zero-crossing disabled */
657             {
658                 if (OutEnd == 2)
659                 {
660                     double r = 0;
661                     if (TP != 0)
662                     {
663                         r = floor((t / TP));
664                     }
665                     t -= static_cast<int>(r) * TP;
666                     inow = nPoints - 1;
667                 }
668                 else
669                 {
670                     inow = nPoints + 1; // Arbitrary value more than nPoints, will be overwritten if needed.
671                 }
672                 // Look in time value table a range to have current time in.
673                 // Beware exact values.
674                 for (int i = 0 ; i < nPoints ; ++i)
675                 {
676                     if (t <= ptr->workt[i])
677                     {
678                         if (t < ptr->workt[i])
679                         {
680                             inow = i - 1;
681                         }
682                         else
683                         {
684                             inow = i;
685                         }
686                         break;
687                     }
688                 }
689             }
690
691             ptr->cnt1 = cnt1;
692             ptr->cnt2 = cnt2;
693             ptr->EVindex = EVindex;
694             ptr->PerEVcnt = PerEVcnt;
695
696             /***************************/
697             /* Hypermatrix case */
698             if (ptr->Hmat)
699             {
700                 for (int j = 0; j < my * ny; ++j)
701                 {
702                     if (ptr->Yt == 1)
703                     {
704                         if (ptr->Yst == 0)
705                         {
706                             /* real case */
707                             y_d = GetRealOutPortPtrs(block, 1);
708                             ptr_d = (double*) ptr->work;
709
710                             if (inow > nPoints)
711                             {
712                                 if (OutEnd == 0)
713                                 {
714                                     y_d[j] = 0; /* Outputs set to zero */
715                                 }
716                                 else if (OutEnd == 1)
717                                 {
718                                     y_d[j] = ptr_d[(nPoints - 1) * ny * my + j]; /* Hold outputs at the end */
719                                 }
720                             }
721                             else
722                             {
723                                 if (inow < 0)
724                                 {
725                                     y_d[j] = 0;
726                                 }
727                                 else
728                                 {
729                                     y_d[j] = ptr_d[inow * ny * my + j];
730                                 }
731                             }
732                         }
733                         else
734                         {
735                             /* complexe case */
736                             y_d = GetRealOutPortPtrs(block, 1);
737                             y_cd = GetImagOutPortPtrs(block, 1);
738                             ptr_d = (double*) ptr->work;
739
740                             if (inow > nPoints)
741                             {
742                                 if (OutEnd == 0)
743                                 {
744                                     y_d[j] = 0; /* Outputs set to zero */
745                                     y_cd[j] = 0; /* Outputs set to zero */
746                                 }
747                                 else if (OutEnd == 1)
748                                 {
749                                     y_d[j] = ptr_d[(nPoints - 1) * ny * my + j]; /* Hold outputs at the end */
750                                     y_cd[j] = ptr_d[nPoints * my * ny + (nPoints - 1) * ny * my + j]; /* Hold outputs at the end */
751                                 }
752                             }
753                             else
754                             {
755                                 if (inow < 0)
756                                 {
757                                     y_d[j] = 0; /* Outputs set to zero */
758                                     y_cd[j] = 0; /* Outputs set to zero */
759                                 }
760                                 else
761                                 {
762                                     y_d[j] = ptr_d[inow * ny * my + j];
763                                     y_cd[j] = ptr_d[nPoints * my * ny + inow * ny * my + j];
764                                 }
765                             }
766                         }
767                     }
768                     else if (ptr->Yt == 8)
769                     {
770                         switch (ptr->Yst)
771                         {
772                             case 1:
773                                 /* --------------------- int8 char  ----------------------------*/
774                                 y_c = Getint8OutPortPtrs(block, 1);
775                                 ptr_c = (char*) ptr->work;
776                                 if (inow > nPoints)
777                                 {
778                                     if (OutEnd == 0)
779                                     {
780                                         y_c[j] = 0; /* Outputs set to zero */
781                                     }
782                                     else if (OutEnd == 1)
783                                     {
784                                         y_c[j] = ptr_c[(nPoints - 1) * ny * my + j]; /* Hold outputs at the end */
785                                     }
786                                 }
787                                 else
788                                 {
789                                     if (inow < 0)
790                                     {
791                                         y_c[j] = 0;
792                                     }
793                                     else
794                                     {
795                                         y_c[j] = ptr_c[inow * ny * my + j];
796                                     }
797                                 }
798                                 break;
799
800                             case 2:
801                                 /* --------------------- int16 short int ---------------------*/
802                                 y_s = Getint16OutPortPtrs(block, 1);
803                                 ptr_s = (short int*) ptr->work;
804                                 if (inow > nPoints)
805                                 {
806                                     if (OutEnd == 0)
807                                     {
808                                         y_s[j] = 0; /* Outputs set to zero */
809                                     }
810                                     else if (OutEnd == 1)
811                                     {
812                                         y_s[j] = ptr_s[(nPoints - 1) * ny * my + j]; /* Hold outputs at the end */
813                                     }
814                                 }
815                                 else
816                                 {
817                                     if (inow < 0)
818                                     {
819                                         y_s[j] = 0;
820                                     }
821                                     else
822                                     {
823                                         y_s[j] = ptr_s[inow * ny * my + j];
824                                     }
825                                 }
826                                 break;
827
828                             case 4:
829                                 /* --------------------- int32 long ---------------------*/
830                                 y_l = Getint32OutPortPtrs(block, 1);
831                                 ptr_l = (int*) ptr->work;
832                                 if (inow > nPoints)
833                                 {
834                                     if (OutEnd == 0)
835                                     {
836                                         y_l[j] = 0; /* Outputs set to zero */
837                                     }
838                                     else if (OutEnd == 1)
839                                     {
840                                         y_l[j] = ptr_l[(nPoints - 1) * ny * my + j]; /* Hold outputs at the end */
841                                     }
842                                 }
843                                 else
844                                 {
845                                     if (inow < 0)
846                                     {
847                                         y_l[j] = 0;
848                                     }
849                                     else
850                                     {
851                                         y_l[j] = ptr_l[inow * ny * my + j];
852                                     }
853                                 }
854                                 break;
855
856                             case 11:
857                                 /*--------------------- uint8 uchar ---------------------*/
858                                 y_uc = Getuint8OutPortPtrs(block, 1);
859                                 ptr_uc = (unsigned char*) ptr->work;
860                                 if (inow > nPoints)
861                                 {
862                                     if (OutEnd == 0)
863                                     {
864                                         y_uc[j] = 0; /* Outputs set to zero */
865                                     }
866                                     else if (OutEnd == 1)
867                                     {
868                                         y_uc[j] = ptr_uc[(nPoints - 1) * ny * my + j]; /* Hold outputs at the end */
869                                     }
870                                 }
871                                 else
872                                 {
873                                     if (inow < 0)
874                                     {
875                                         y_uc[j] = 0;
876                                     }
877                                     else
878                                     {
879                                         y_uc[j] = ptr_uc[inow * ny * my + j];
880                                     }
881                                 }
882                                 break;
883
884                             case 12:
885                                 /* --------------------- uint16 ushort ---------------------*/
886                                 y_us = Getuint16OutPortPtrs(block, 1);
887                                 ptr_us = (unsigned short int*) ptr->work;
888                                 if (inow > nPoints)
889                                 {
890                                     if (OutEnd == 0)
891                                     {
892                                         y_us[j] = 0; /* Outputs set to zero */
893                                     }
894                                     else if (OutEnd == 1)
895                                     {
896                                         y_us[j] = ptr_us[(nPoints - 1) * ny * my + j]; /* Hold outputs at the end */
897                                     }
898                                 }
899                                 else
900                                 {
901                                     if (inow < 0)
902                                     {
903                                         y_us[j] = 0;
904                                     }
905                                     else
906                                     {
907                                         y_us[j] = ptr_us[inow * ny * my + j];
908                                     }
909                                 }
910                                 break;
911
912                             case 14:
913                                 /* --------------------- uint32 ulong ---------------------*/
914                                 y_ul = Getuint32OutPortPtrs(block, 1);
915                                 ptr_ul = (unsigned int*) ptr->work;
916                                 if (inow > nPoints)
917                                 {
918                                     if (OutEnd == 0)
919                                     {
920                                         y_ul[j] = 0; /* Outputs set to zero */
921                                     }
922                                     else if (OutEnd == 1)
923                                     {
924                                         y_ul[j] = ptr_ul[(nPoints - 1) * ny * my + j]; /* Hold outputs at the end */
925                                     }
926                                 }
927                                 else
928                                 {
929                                     if (inow < 0)
930                                     {
931                                         y_ul[j] = 0;
932                                     }
933                                     else
934                                     {
935                                         y_ul[j] = ptr_ul[inow * ny * my + j];
936                                     }
937                                 }
938                                 break;
939                         }
940                     }
941                 } /* for j loop */
942             }
943             /****************************/
944             /* Scalar or vectorial case */
945             else
946             {
947                 for (int j = 0; j < my; ++j)
948                 {
949                     double y1, y2, d1, d2, h, dh, ddh, dddh;
950                     if (ptr->Yt == 1)
951                     {
952                         if ((ptr->Yst == 0) || (ptr->Yst == 1))
953                         {
954                             /*  If real or complex*/
955                             y_d = GetRealOutPortPtrs(block, 1);
956                             ptr_d = (double*) ptr->work;
957                             ptr_D = (double*) ptr->D;
958
959                             if (inow > nPoints)
960                             {
961                                 if (OutEnd == 0)
962                                 {
963                                     y_d[j] = 0.0; /* Outputs set to zero */
964                                 }
965                                 else if (OutEnd == 1)
966                                 {
967                                     y_d[j] = ptr_d[nPoints - 1 + (j) * nPoints]; /* Hold outputs at the end */
968                                 }
969                             }
970                             else if (Method == 0)
971                             {
972                                 if (inow < 0)
973                                 {
974                                     y_d[j] = 0.0;
975                                 }
976                                 else
977                                 {
978                                     y_d[j] = ptr_d[inow + (j) * nPoints];
979                                 }
980                             }
981                             else if (Method == 1)
982                             {
983                                 if (inow < 0)
984                                 {
985                                     inow = 0;
986                                 }
987                                 t1 = ptr->workt[inow];
988                                 t2 = ptr->workt[inow + 1];
989                                 y1 = ptr_d[inow + j * nPoints];
990                                 y2 = ptr_d[inow + 1 + j * nPoints];
991                                 y_d[j] = (y2 - y1) * (t - t1) / (t2 - t1) + y1;
992                             }
993                             else if (Method >= 2)
994                             {
995                                 t1 = ptr->workt[inow];
996                                 t2 = ptr->workt[inow + 1];
997                                 y1 = ptr_d[inow + j * nPoints];
998                                 y2 = ptr_d[inow + 1 + j * nPoints];
999                                 d1 = ptr_D[inow + j * nPoints];
1000                                 d2 = ptr_D[inow + 1 + j * nPoints];
1001                                 scicos_evalhermite(&t, &t1, &t2, &y1, &y2, &d1, &d2, &h, &dh, &ddh, &dddh, &inow);
1002                                 y_d[j] = h;
1003                             }
1004                         }
1005                         if (ptr->Yst == 1)
1006                         {
1007                             /*  -------------- complex ----------------------*/
1008                             y_cd = GetImagOutPortPtrs(block, 1);
1009                             if (inow > nPoints)
1010                             {
1011                                 if (OutEnd == 0)
1012                                 {
1013                                     y_cd[j] = 0.0; /* Outputs set to zero*/
1014                                 }
1015                                 else if (OutEnd == 1)
1016                                 {
1017                                     y_cd[j] = ptr_d[nPoints * my + nPoints - 1 + (j) * nPoints]; // Hold outputs at the end
1018                                 }
1019                             }
1020                             else if (Method == 0)
1021                             {
1022                                 if (inow < 0)
1023                                 {
1024                                     y_cd[j] = 0.0; /* Outputs set to zero */
1025                                 }
1026                                 else
1027                                 {
1028                                     y_cd[j] = ptr_d[nPoints * my + inow + (j) * nPoints];
1029                                 }
1030                             }
1031                             else if (Method == 1)
1032                             {
1033                                 if (inow < 0)
1034                                 {
1035                                     inow = 0;
1036                                 } /* Extrapolation for 0<t<X(0) */
1037                                 t1 = ptr->workt[inow];
1038                                 t2 = ptr->workt[inow + 1];
1039                                 y1 = ptr_d[nPoints * my + inow + j * nPoints];
1040                                 y2 = ptr_d[nPoints * my + inow + 1 + j * nPoints];
1041                                 y_cd[j] = (y2 - y1) * (t - t1) / (t2 - t1) + y1;
1042                             }
1043                             else if (Method >= 2)
1044                             {
1045                                 t1 = ptr->workt[inow];
1046                                 t2 = ptr->workt[inow + 1];
1047                                 y1 = ptr_d[inow + j * nPoints + nPoints];
1048                                 y2 = ptr_d[inow + 1 + j * nPoints + nPoints];
1049                                 d1 = ptr_D[inow + j * nPoints + nPoints];
1050                                 d2 = ptr_D[inow + 1 + j * nPoints + nPoints];
1051                                 scicos_evalhermite(&t, &t1, &t2, &y1, &y2, &d1, &d2, &h, &dh, &ddh, &dddh, &inow);
1052                                 y_cd[j] = h;
1053                             }
1054                         }
1055                     }
1056                     else if (ptr->Yt == 8)
1057                     {
1058                         switch (ptr->Yst)
1059                         {
1060                             case 1: /* --------------------- int8 char  ----------------------------*/
1061                                 y_c = Getint8OutPortPtrs(block, 1);
1062                                 ptr_c = (char*) ptr->work;
1063                                 /* y_c[j]=ptr_c[inow+(j)*nPoints]; */
1064                                 if (inow > nPoints)
1065                                 {
1066                                     if (OutEnd == 0)
1067                                     {
1068                                         y_c[j] = 0; /* Outputs set to zero */
1069                                     }
1070                                     else if (OutEnd == 1)
1071                                     {
1072                                         y_c[j] = ptr_c[nPoints - 1 + (j) * nPoints]; /* Hold outputs at the end */
1073                                     }
1074                                 }
1075                                 else if (Method == 0)
1076                                 {
1077                                     if (inow < 0)
1078                                     {
1079                                         y_c[j] = 0;
1080                                     }
1081                                     else
1082                                     {
1083                                         y_c[j] = ptr_c[inow + (j) * nPoints];
1084                                     }
1085                                 }
1086                                 else if (Method >= 1)
1087                                 {
1088                                     if (inow < 0)
1089                                     {
1090                                         inow = 0;
1091                                     }
1092                                     t1 = ptr->workt[inow];
1093                                     t2 = ptr->workt[inow + 1];
1094                                     y1 = (double)ptr_c[inow + j * nPoints];
1095                                     y2 = (double)ptr_c[inow + 1 + j * nPoints];
1096                                     y_c[j] = (char)((y2 - y1) * (t - t1) / (t2 - t1) + y1);
1097                                 }
1098                                 break;
1099                             case 2:
1100                                 /* --------------------- int16 short ---------------------*/
1101                                 y_s = Getint16OutPortPtrs(block, 1);
1102                                 ptr_s = (short int*) ptr->work;
1103                                 /* y_s[j]=ptr_s[inow+(j)*nPoints]; */
1104                                 if (inow > nPoints)
1105                                 {
1106                                     if (OutEnd == 0)
1107                                     {
1108                                         y_s[j] = 0; /* Outputs set to zero */
1109                                     }
1110                                     else if (OutEnd == 1)
1111                                     {
1112                                         y_s[j] = ptr_s[nPoints - 1 + (j) * nPoints]; // Hold outputs at the end
1113                                     }
1114                                 }
1115                                 else if (Method == 0)
1116                                 {
1117                                     if (inow < 0)
1118                                     {
1119                                         y_s[j] = 0;
1120                                     }
1121                                     else
1122                                     {
1123                                         y_s[j] = ptr_s[inow + (j) * nPoints];
1124                                     }
1125                                 }
1126                                 else if (Method >= 1)
1127                                 {
1128                                     if (inow < 0)
1129                                     {
1130                                         inow = 0;
1131                                     }
1132                                     t1 = ptr->workt[inow];
1133                                     t2 = ptr->workt[inow + 1];
1134                                     y1 = (double)ptr_s[inow + j * nPoints];
1135                                     y2 = (double)ptr_s[inow + 1 + j * nPoints];
1136                                     y_s[j] = (short int)((y2 - y1) * (t - t1) / (t2 - t1) + y1);
1137                                 }
1138                                 break;
1139                             case 4:
1140                                 /* --------------------- int32 long ---------------------*/
1141                                 y_l = Getint32OutPortPtrs(block, 1);
1142                                 ptr_l = (int*) ptr->work;
1143                                 /* y_l[j]=ptr_l[inow+(j)*nPoints]; */
1144                                 if (inow > nPoints)
1145                                 {
1146                                     if (OutEnd == 0)
1147                                     {
1148                                         y_l[j] = 0; /* Outputs set to zero */
1149                                     }
1150                                     else if (OutEnd == 1)
1151                                     {
1152                                         y_l[j] = ptr_l[nPoints - 1 + (j) * nPoints]; /* Hold outputs at the end */
1153                                     }
1154                                 }
1155                                 else if (Method == 0)
1156                                 {
1157                                     if (inow < 0)
1158                                     {
1159                                         y_l[j] = 0;
1160                                     }
1161                                     else
1162                                     {
1163                                         y_l[j] = ptr_l[inow + (j) * nPoints];
1164                                     }
1165                                 }
1166                                 else if (Method >= 1)
1167                                 {
1168                                     t1 = ptr->workt[inow];
1169                                     t2 = ptr->workt[inow + 1];
1170                                     y1 = (double)ptr_l[inow + j * nPoints];
1171                                     y2 = (double)ptr_l[inow + 1 + j * nPoints];
1172                                     y_l[j] = (int)((y2 - y1) * (t - t1) / (t2 - t1) + y1);
1173                                 }
1174                                 break;
1175                             case 11: /*--------------------- uint8 uchar ---------------------*/
1176                                 y_uc = Getuint8OutPortPtrs(block, 1);
1177                                 ptr_uc = (unsigned char*) ptr->work;
1178                                 /* y_uc[j]=ptr_uc[inow+(j)*nPoints]; */
1179                                 if (inow > nPoints)
1180                                 {
1181                                     if (OutEnd == 0)
1182                                     {
1183                                         y_uc[j] = 0; /* Outputs set to zero */
1184                                     }
1185                                     else if (OutEnd == 1)
1186                                     {
1187                                         y_uc[j] = ptr_uc[nPoints - 1 + (j) * nPoints]; /* Hold outputs at the end */
1188                                     }
1189                                 }
1190                                 else if (Method == 0)
1191                                 {
1192                                     if (inow < 0)
1193                                     {
1194                                         y_uc[j] = 0;
1195                                     }
1196                                     else
1197                                     {
1198                                         y_uc[j] = ptr_uc[inow + (j) * nPoints];
1199                                     }
1200                                 }
1201                                 else if (Method >= 1)
1202                                 {
1203                                     t1 = ptr->workt[inow];
1204                                     t2 = ptr->workt[inow + 1];
1205                                     y1 = (double)ptr_uc[inow + j * nPoints];
1206                                     y2 = (double)ptr_uc[inow + 1 + j * nPoints];
1207                                     y_uc[j] = (unsigned char)((y2 - y1) * (t - t1) / (t2 - t1) + y1);
1208                                 }
1209                                 break;
1210                             case 12:
1211                                 /* --------------------- uint16 ushort int ---------------------*/
1212                                 y_us = Getuint16OutPortPtrs(block, 1);
1213                                 ptr_us = (unsigned short int*) ptr->work;
1214                                 /* y_us[j]=ptr_us[inow+(j)*nPoints]; */
1215                                 if (inow > nPoints)
1216                                 {
1217                                     if (OutEnd == 0)
1218                                     {
1219                                         y_us[j] = 0; /* Outputs set to zero */
1220                                     }
1221                                     else if (OutEnd == 1)
1222                                     {
1223                                         y_us[j] = ptr_us[nPoints - 1 + (j) * nPoints]; /* Hold outputs at the end */
1224                                     }
1225                                 }
1226                                 else if (Method == 0)
1227                                 {
1228                                     if (inow < 0)
1229                                     {
1230                                         y_us[j] = 0;
1231                                     }
1232                                     else
1233                                     {
1234                                         y_us[j] = ptr_us[inow + (j) * nPoints];
1235                                     }
1236                                 }
1237                                 else if (Method >= 1)
1238                                 {
1239                                     t1 = ptr->workt[inow];
1240                                     t2 = ptr->workt[inow + 1];
1241                                     y1 = (double)ptr_us[inow + j * nPoints];
1242                                     y2 = (double)ptr_us[inow + 1 + j * nPoints];
1243                                     y_us[j] = (unsigned short int)((y2 - y1) * (t - t1) / (t2 - t1) + y1);
1244                                 }
1245                                 break;
1246                             case 14:
1247                                 /* --------------------- uint32 ulong ---------------------*/
1248                                 y_ul = Getuint32OutPortPtrs(block, 1);
1249                                 ptr_ul = (unsigned int*) ptr->work;
1250                                 /* y_ul[j]=ptr_ul[inow+(j)*nPoints]; */
1251                                 if (inow > nPoints)
1252                                 {
1253                                     if (OutEnd == 0)
1254                                     {
1255                                         y_ul[j] = 0; /* Outputs set to zero */
1256                                     }
1257                                     else if (OutEnd == 1)
1258                                     {
1259                                         y_ul[j] = ptr_ul[nPoints - 1 + (j) * nPoints]; /* Hold outputs at the end */
1260                                     }
1261                                 }
1262                                 else if (Method == 0)
1263                                 {
1264                                     if (inow < 0)
1265                                     {
1266                                         y_ul[j] = 0;
1267                                     }
1268                                     else
1269                                     {
1270                                         y_ul[j] = ptr_ul[inow + (j) * nPoints];
1271                                     }
1272                                 }
1273                                 else if (Method >= 1)
1274                                 {
1275                                     t1 = ptr->workt[inow];
1276                                     t2 = ptr->workt[inow + 1];
1277                                     y1 = (double)ptr_ul[inow + j * nPoints];
1278                                     y2 = (double)ptr_ul[inow + 1 + j * nPoints];
1279                                     y_ul[j] = (unsigned int)((y2 - y1) * (t - t1) / (t2 - t1) + y1);
1280                                 }
1281                                 break;
1282                         }
1283                     }
1284                 } /* for j loop */
1285             }
1286             /********************************************************************/
1287             break;
1288         }
1289         case 3 :
1290         {
1291             /* Event date computation */
1292             /* Retrieve 'ptr' of the structure of the block */
1293             ptr = *work;
1294             int nPoints = ptr->nPoints;
1295             int cnt1 = ptr->cnt1;
1296             int cnt2 = ptr->cnt2;
1297             int EVindex = ptr->EVindex;
1298             int PerEVcnt = ptr->PerEVcnt;
1299
1300             /* Get current simulation time */
1301             //double t = get_scicos_time();
1302
1303             double TNm1  = ptr->workt[nPoints - 1];
1304             double TP    = TNm1 - 0;
1305
1306             int jfirst;
1307             if (ZC == 1)
1308             {
1309                 /* Generate Events only if ZC is active */
1310                 if ((Method == 1) || (Method == 0))
1311                 {
1312                     /*-------------------------*/
1313                     if (ptr->firstevent == 1)
1314                     {
1315                         jfirst = nPoints - 1; /* Finding first positive time instant */
1316                         for (int j = 0; j < nPoints; ++j)
1317                         {
1318                             if (ptr->workt[j] > 0)
1319                             {
1320                                 jfirst = j;
1321                                 break;
1322                             }
1323                         }
1324                         block->evout[0] = ptr->workt[jfirst];
1325                         EVindex = jfirst;
1326                         ptr->EVindex = EVindex;
1327                         ptr->firstevent = 0;
1328                         return;
1329                     }
1330                     /*------------------------*/
1331                     int i = EVindex;
1332                     /*------------------------*/
1333                     if (i < nPoints - 1)
1334                     {
1335                         block->evout[0] = ptr->workt[i + 1] - ptr->workt[i];
1336                         EVindex = i + 1;
1337                     }
1338                     /*------------------------*/
1339                     if (i == nPoints - 1)
1340                     {
1341                         if (OutEnd == 2)
1342                         {
1343                             /*  Periodic */
1344                             cnt1 = -1;
1345                             cnt2 = 0;
1346                             PerEVcnt++; /* When OutEnd==2 (perodic output) */
1347                             jfirst = nPoints - 1; /* Finding first positive time instant */
1348                             for (int j = 0; j < nPoints; ++j)
1349                             {
1350                                 if (ptr->workt[j] >= 0)
1351                                 {
1352                                     jfirst = j;
1353                                     break;
1354                                 }
1355                             }
1356                             block->evout[0] = ptr->workt[jfirst];
1357                             EVindex = jfirst;
1358                         }
1359                     }
1360                     /*-------------------------- */
1361                 }
1362                 else if (Method <= 3)
1363                 {
1364                     if (ptr->firstevent == 1)
1365                     {
1366                         block->evout[0] = TP;
1367                         ptr->firstevent = 0;
1368                     }
1369                     else
1370                     {
1371                         if (OutEnd == 2)
1372                         {
1373                             block->evout[0] = TP;
1374                         }
1375                         PerEVcnt++;
1376                     }
1377                     cnt1 = -1;
1378                     cnt2 = 0;
1379                 }
1380                 ptr->cnt1 = cnt1;
1381                 ptr->cnt2 = cnt2;
1382                 ptr->EVindex = EVindex;
1383                 ptr->PerEVcnt = PerEVcnt;
1384             }
1385             /***********************************************************************/
1386             break;
1387         }
1388         case 5 :
1389         {
1390             /* Finish */
1391             ptr = *work;
1392             if (ptr != nullptr)
1393             {
1394                 if (ptr->D != nullptr)
1395                 {
1396                     delete[] ptr->D;
1397                 }
1398                 if (ptr->work != nullptr)
1399                 {
1400                     delete[] (char*) ptr->work;
1401                 }
1402                 if (ptr->workt != nullptr)
1403                 {
1404                     delete[] ptr->workt;
1405                 }
1406                 delete[] ptr;
1407             }
1408             break;
1409             /***********************************************************************/
1410         }
1411         case 6 :
1412         {
1413             /* Re-init */
1414
1415             // Finish then init
1416             fromws_c(block, 5);
1417             fromws_c(block, 4);
1418             break;
1419         }
1420         /*************************************************************************/
1421         default :
1422             return;
1423     }
1424 }