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 }