Xcos solvers : Sundials explicit Runge-Kutta 4(5) implem and interface
[scilab.git] / scilab / modules / scicos / src / c / scicos.c
1 /*  Scicos
2 *
3 *  Copyright (C) INRIA - METALAU Project <scicos@inria.fr>
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 *
19 * See the file ./license.txt
20 */
21 /* 11-03-2005,  Masoud
22 * adding A-Jacobian
23 * istate =-1 case;
24 *
25 * xx/03/06, Alan
26 * enable matrix typed
27 * input/output regular ports
28 *
29 * xx/02/07, Alan
30 * added object paramaters/states
31 *
32 * 20/07/07,  Masoud
33 * CVODE (sundials) replaced LSODAR
34 * IDA  (Sundials)  replaced DASKR
35 */
36 /*--------------------------------------------------------------------------*/
37 #include <stdlib.h>
38 #include <string.h>
39 #include <math.h>
40
41 /* Sundials includes */
42 #include <cvode/cvode.h>           /* prototypes for CVODES fcts. and consts. */
43 #include <cvode/cvode_dense.h>     /* prototype for CVDense */
44 #include <ida/ida.h>
45 #include <ida/ida_dense.h>
46 #include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., and macros */
47 #include <sundials/sundials_dense.h> /* definitions DenseMat and DENSE_ELEM */
48 #include <sundials/sundials_types.h> /* definition of type realtype */
49 #include <sundials/sundials_math.h>
50 #include <kinsol/kinsol.h>
51 #include <kinsol/kinsol_dense.h>
52 #include <sundials/sundials_extension.h> /* uses extension for scicos */
53 #include "ida_impl.h"
54
55 #include "machine.h" /* C2F */
56 #include "dynamic_link.h"
57 #include "scicos-def.h"
58 #include "stack-def.h"
59 #include "sciprint.h"
60 #include "scicos.h"
61 #include "import.h"
62 #include "blocks.h"
63 #include "core_math.h"
64 #include "storeCommand.h"
65 #include "syncexec.h"
66 #include "realtime.h"
67 #include "math_graphics.h"
68 #include "MALLOC.h"
69 #include "cvstr.h"
70 #include "ezxml.h"
71 #include "xscion.h"
72
73 #include "sciblk2.h"
74 #include "sciblk4.h"
75 #include "dynlib_scicos.h"
76
77 #include "lsodar.h"           /* prototypes for lsodar fcts. and consts. */
78
79 #if defined(linux) && defined(__i386__)
80 #include "setPrecisionFPU.h"
81 #endif
82
83 #include "localization.h"
84 #include "charEncoding.h"
85 /*--------------------------------------------------------------------------*/
86 typedef struct
87 {
88     void *ida_mem;
89     N_Vector ewt;
90     double *rwork;
91     int *iwork;
92     double *gwork; /* just added for a very special use: a
93                                    space passing to grblkdakr for zero crossing surfaces
94                                    when updating mode variables during initialization */
95 } *UserData;
96
97 SCICOS_IMPEXP SCSPTR_struct C2F(scsptr);
98 /*--------------------------------------------------------------------------*/
99
100 #define freeall                                 \
101         if (*neq>0) if (C2F(cmsolver).solver) (CVodeFree(&cvode_mem)); else (LSodarFree(&cvode_mem));           \
102         if (*neq>0) N_VDestroy_Serial(y);               \
103         if ( ng>0 ) FREE(jroot);                        \
104         if ( ng>0 ) FREE(zcros);
105
106
107 /* TJacque allocates by sundials */
108 #define freeallx                                \
109         if (*neq>0) free(TJacque);      \
110         if (*neq>0) FREE(data->rwork);          \
111         if (( ng>0 )&& (*neq>0)) FREE(data->gwork);     \
112         if (*neq>0) N_VDestroy_Serial(data->ewt);       \
113         if (*neq>0) FREE(data);                 \
114         if (*neq>0) IDAFree(&ida_mem);          \
115         if (*neq>0) N_VDestroy_Serial(IDx);             \
116         if (*neq>0) N_VDestroy_Serial(yp);              \
117         if (*neq>0) N_VDestroy_Serial(yy);              \
118         if ( ng>0 ) FREE(jroot);                        \
119         if ( ng>0 ) FREE(zcros);                        \
120         if (nmod>0) FREE(Mode_save);
121
122 #define freeouttbptr                            \
123         FREE(outtbd);                                   \
124         FREE(outtbc);                                   \
125         FREE(outtbs);                                   \
126         FREE(outtbl);                                   \
127         FREE(outtbuc);                          \
128         FREE(outtbus);                          \
129         FREE(outtbul);
130
131 #define freekinsol                              \
132         FREE(Mode_save);                                \
133         N_VDestroy_Serial(y);                           \
134         N_VDestroy_Serial(fscale);                      \
135         N_VDestroy_Serial(yscale);                      \
136         KINFree(&kin_mem);
137
138
139 #define ONE   RCONST(1.0)
140 #define ZERO  RCONST(0.0)
141 #define T0    RCONST(0.0)
142 /*--------------------------------------------------------------------------*/
143 /* Table of constant values */
144 static int c__90 = 90;
145 static int c__91 = 91;
146 static int c__0 = 0;
147 static double c_b14 = 0.;
148 static int c__1 = 1;
149
150 int TCritWarning = 0;
151
152 static double *t0 = NULL, *tf = NULL;
153 static double *x = NULL, *tevts = NULL, *xd = NULL, *g = NULL;
154 static  double Atol = 0., rtol = 0., ttol = 0., deltat = 0., hmax = 0.;
155 static int *ierr = NULL;
156 static int *pointi = NULL;
157 static int *xptr = NULL, *modptr = NULL, *evtspt = NULL;
158 static int *funtyp = NULL, *inpptr = NULL, *outptr = NULL, *inplnk = NULL, *outlnk = NULL;
159 static int *clkptr = NULL, *ordptr = NULL, *ordclk = NULL, *cord = NULL, *iord = NULL, *oord = NULL,  *zord = NULL,  *critev = NULL,  *zcptr = NULL;
160 static int nblk = 0, nordptr = 0, nlnk = 0, nx = 0, ng = 0, ncord = 0, noord = 0, nzord = 0;
161 static int nordclk = 0, niord = 0, nmod = 0;
162 static int nelem = 0;
163 static int *mod = NULL;
164 static int *neq = NULL;
165 static int *xprop = NULL; /* xproperty */
166 static int debug_block = 0;
167 static int *iwa = NULL;
168 static int hot;
169
170 /* declaration of ptr for typed port */
171 static void **outtbptr = NULL; /*pointer array of object of outtb*/
172 static int *outtbsz = NULL;    /*size of object of outtb*/
173 static int *outtbtyp = NULL;   /*type of object of outtb*/
174
175 SCSREAL_COP *outtbdptr = NULL;     /*to store double of outtb*/
176 SCSINT8_COP *outtbcptr = NULL;     /*to store int8 of outtb*/
177 SCSINT16_COP *outtbsptr = NULL;    /*to store int16 of outtb*/
178 SCSINT32_COP *outtblptr = NULL;    /*to store int32 of outtb*/
179 SCSUINT8_COP *outtbucptr = NULL;   /*to store unsigned int8 of outtb */
180 SCSUINT16_COP *outtbusptr = NULL;  /*to store unsigned int16 of outtb */
181 SCSUINT32_COP *outtbulptr = NULL;  /*to store unsigned int32 of outtb */
182
183 static outtb_el *outtb_elem = NULL;
184
185 static scicos_block *Blocks = NULL;
186
187 /* pass to external variable for code generation */
188 /* reserved variable name */
189 int *block_error = NULL;
190 double scicos_time = 0.;
191 int phase = 0;
192 int Jacobian_Flag = 0;
193 // double CI = 0., CJ = 0.;  // doubles returned by Get_Jacobian_ci and Get_Jacobian_cj respectively
194 double  CJJ = 0.;            // returned by Get_Jacobian_parameter
195 double SQuround = 0.;
196 /* Jacobian*/
197 static int AJacobian_block = 0;
198
199
200 /* Variable declaration moved to scicos.c because it was in the scicos-def.h therefore
201 * multiple declaration of the variable and linkers were complaining about duplicate
202 * symbols
203 */
204 SCICOS_IMPEXP COSDEBUGCOUNTER_struct C2F(cosdebugcounter);
205 SCICOS_IMPEXP RTFACTOR_struct C2F(rtfactor);
206 SCICOS_IMPEXP SOLVER_struct C2F(cmsolver);
207 SCICOS_IMPEXP CURBLK_struct C2F(curblk);
208 SCICOS_IMPEXP COSDEBUG_struct C2F(cosdebug);
209 SCICOS_IMPEXP COSHLT_struct C2F(coshlt);
210 SCICOS_IMPEXP DBCOS_struct C2F(dbcos);
211 SCICOS_IMPEXP COSTOL_struct C2F(costol);
212 SCICOS_IMPEXP COSERR_struct coserr;
213 /*--------------------------------------------------------------------------*/
214 static void FREE_blocks();
215 static void cosini(double *told);
216 static void cosend(double *told);
217 static void cossim(double *told);
218 static void cossimdaskr(double *told);
219 static void doit(double *told);
220 static void idoit(double *told);
221 static void zdoit(double *told, double *xt, double *xtd, double *g);
222 static void cdoit(double *told);
223 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa);
224 static void odoit(double *told, double *xt, double *xtd, double *residual);
225 static void ddoit(double *told);
226 static void edoit(double *told, int *kiwa);
227 static void reinitdoit(double *told);
228 static int CallKinsol(double *told);
229 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data);
230 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);
231 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data);
232 static int simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res);
233 static int grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngroot, realtype * res);
234 static void addevs(double t, int *evtnb, int *ierr1);
235 static int synchro_g_nev(ScicosImport *scs_imp, double *g, int kf, int *ierr);
236 static void Multp(double *A, double *B, double *R, int ra, int rb, int ca, int cb);
237 static int read_id(ezxml_t *elements, char *id, double *value);
238 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata);
239 static int Jacobians(long int Neq, realtype tt, N_Vector yy, N_Vector yp,
240                      N_Vector resvec, realtype cj, void *jdata, DenseMat Jacque,
241                      N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
242 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk);
243 static int synchro_nev(ScicosImport *scs_imp, int kf, int *ierr);
244 /*--------------------------------------------------------------------------*/
245 extern int C2F(dset)(int *n, double *dx, double *dy, int *incy);
246 extern int C2F(dcopy)(int *, double *, int *, double *, int *);
247 extern int C2F(msgs)();
248 extern int C2F(getscsmax)();
249 extern int C2F(makescicosimport)();
250 extern int C2F(clearscicosimport)();
251 /*--------------------------------------------------------------------------*/
252 void putevs(double *t, int *evtnb, int *ierr1);
253 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job);
254 int simblkKinsol(N_Vector yy, N_Vector resval, void *rdata);
255 /*--------------------------------------------------------------------------*/
256 int C2F(scicos)(double *x_in, int *xptr_in, double *z__,
257                 void **work, int *zptr, int *modptr_in,
258                 void **oz, int *ozsz, int *oztyp, int *ozptr,
259                 int *iz, int *izptr, double *t0_in,
260                 double *tf_in, double *tevts_in, int *evtspt_in,
261                 int *nevts, int *pointi_in, void **outtbptr_in,
262                 int *outtbsz_in, int *outtbtyp_in,
263                 outtb_el *outtb_elem_in, int *nelem1, int *nlnk1,
264                 int *funptr, int *funtyp_in, int *inpptr_in,
265                 int *outptr_in, int *inplnk_in, int *outlnk_in,
266                 double *rpar, int *rpptr, int *ipar, int *ipptr,
267                 void **opar, int *oparsz, int *opartyp, int *opptr,
268                 int *clkptr_in, int *ordptr_in, int *nordptr1,
269                 int *ordclk_in, int *cord_in, int *ncord1,
270                 int *iord_in, int *niord1, int *oord_in,
271                 int *noord1, int *zord_in, int *nzord1,
272                 int *critev_in, int *nblk1, int *ztyp,
273                 int *zcptr_in, int *subscr, int *nsubs,
274                 double *simpar, int *flag__, int *ierr_out)
275 {
276     int i1, kf, lprt, in, out, job = 1;
277
278
279     static int mxtb = 0, ierr0 = 0, kfun0 = 0, i = 0, j = 0, k = 0, jj = 0;
280     static int ni = 0, no = 0;
281     static int nz = 0, noz = 0, nopar = 0;
282     double *W = NULL;
283
284     // Set FPU Flag to Extended for scicos simulation
285     // in order to override Java setting it to Double.
286 #if defined(linux) && defined(__i386__)
287     setFPUToExtended();
288 #endif
289     /*     Copyright INRIA */
290     /* iz,izptr are used to pass block labels */
291     TCritWarning = 0;
292
293     t0 = t0_in;
294     tf = tf_in;
295     ierr = ierr_out;
296
297     /* Parameter adjustments */
298     pointi = pointi_in;
299     x = x_in;
300     xptr = xptr_in - 1;
301     modptr = modptr_in - 1;
302     --zptr;
303     --izptr;
304     --ozptr;
305     evtspt = evtspt_in - 1;
306     tevts = tevts_in - 1;
307     outtbptr = outtbptr_in;
308     outtbsz = outtbsz_in;
309     outtbtyp = outtbtyp_in;
310     outtb_elem = outtb_elem_in;
311     funtyp = funtyp_in - 1;
312     inpptr = inpptr_in - 1;
313     outptr = outptr_in - 1;
314     inplnk = inplnk_in - 1;
315     outlnk = outlnk_in - 1;
316     --rpptr;
317     --ipptr;
318     --opptr;
319     clkptr = clkptr_in - 1;
320     ordptr = ordptr_in - 1;
321     ordclk = ordclk_in - 1;
322     cord = cord_in - 1;
323     iord = iord_in - 1;
324     oord = oord_in - 1;
325     zord = zord_in - 1;
326
327     critev = critev_in - 1;
328     --ztyp;
329     zcptr = zcptr_in - 1;
330     --simpar;
331
332     /* Function Body */
333     Atol = simpar[1];
334     rtol = simpar[2];
335     ttol = simpar[3];
336     deltat = simpar[4];
337     C2F(rtfactor).scale = simpar[5];
338     C2F(cmsolver).solver = (int) simpar[6];
339     hmax = simpar[7];
340
341     nordptr = *nordptr1;
342     nblk  = *nblk1;
343     ncord = *ncord1;
344     noord = *noord1;
345     nzord = *nzord1;
346     niord = *niord1;
347     nlnk  = *nlnk1;
348     nelem = *nelem1;
349     *ierr = 0;
350
351     nordclk = ordptr[nordptr] - 1;  /* number of rows in ordclk is ordptr(nclkp1)-1 */
352     ng      = zcptr[nblk + 1] - 1;  /* computes number of zero crossing surfaces */
353     nmod    = modptr[nblk + 1] - 1; /* computes number of modes */
354     nz      = zptr[nblk + 1] - 1;   /* number of discrete real states */
355     noz     = ozptr[nblk + 1] - 1;  /* number of discrete object states */
356     nopar   = opptr[nblk + 1] - 1;  /* number of object parameters */
357     nx      = xptr[nblk + 1] - 1;   /* number of object parameters */
358     neq     = &nx;
359
360     xd = &x[xptr[nblk + 1] - 1];
361
362     /* check for hard coded maxsize */
363     for (i = 1; i <= nblk; ++i)
364     {
365         if (funtyp[i] < 10000)
366         {
367             funtyp[i] %= 1000;
368         }
369         else
370         {
371             funtyp[i] = funtyp[i] % 1000 + 10000;
372         }
373         ni = inpptr[i + 1] - inpptr[i];
374         no = outptr[i + 1] - outptr[i];
375         if (funtyp[i] == 1)
376         {
377             if (ni + no > 11)
378             {
379                 /*     hard coded maxsize in callf.c */
380                 C2F(msgs)(&c__90, &c__0);
381                 C2F(curblk).kfun = i;
382                 *ierr = i + 1005;
383                 return 0;
384             }
385         }
386         else if (funtyp[i] == 2 || funtyp[i] == 3)
387         {
388             /*     hard coded maxsize in scicos.h */
389             if (ni + no > SZ_SIZE)
390             {
391                 C2F(msgs)(&c__90, &c__0);
392                 C2F(curblk).kfun = i;
393                 *ierr = i + 1005;
394                 return 0;
395             }
396         }
397         mxtb = 0;
398         if (funtyp[i] == 0)
399         {
400             if (ni > 1)
401             {
402                 for (j = 1; j <= ni; ++j)
403                 {
404                     k = inplnk[inpptr[i] - 1 + j];
405                     mxtb = mxtb + (outtbsz[k - 1] * outtbsz[(k - 1) + nlnk]);
406                 }
407             }
408             if (no > 1)
409             {
410                 for (j = 1; j <= no; ++j)
411                 {
412                     k = outlnk[outptr[i] - 1 + j];
413                     mxtb = mxtb + (outtbsz[k - 1] * outtbsz[(k - 1) + nlnk]);
414                 }
415             }
416             if (mxtb > TB_SIZE)
417             {
418                 C2F(msgs)(&c__91, &c__0);
419                 C2F(curblk).kfun = i;
420                 *ierr = i + 1005;
421                 return 0;
422             }
423         }
424     }
425
426     if (nx > 0) /* xprop */
427     {
428         if ((xprop = MALLOC(sizeof(int) * nx)) == NULL )
429         {
430             *ierr = 5;
431             return 0;
432         }
433     }
434     for (i = 0; i < nx; i++) /* initialize */
435     {
436         xprop[i] = 1;
437     }
438     if (nmod > 0) /* mod */
439     {
440         if ((mod = MALLOC(sizeof(int) * nmod)) == NULL )
441         {
442             *ierr = 5;
443             if (nx > 0) FREE(xprop);
444             return 0;
445         }
446     }
447     if (ng > 0) /* g becomes global */
448     {
449         if ((g = MALLOC(sizeof(double) * ng)) == NULL )
450         {
451             *ierr = 5;
452             if (nmod > 0) FREE(mod);
453             if (nx > 0) FREE(xprop);
454             return 0;
455         }
456     }
457
458     debug_block = -1; /* no debug block for start */
459     C2F(cosdebugcounter).counter = 0;
460
461     /** Create Block's array **/
462     if ((Blocks = MALLOC(sizeof(scicos_block) * nblk)) == NULL )
463     {
464         *ierr = 5;
465         if (nx > 0) FREE(xprop);
466         if (nmod > 0) FREE(mod);
467         if (ng > 0) FREE(g);
468         return 0;
469     }
470
471     /** Setting blocks properties for each entry in Block's array **/
472
473     /* 1 : type and pointer on simulation function */
474     for (kf = 0; kf < nblk; ++kf)   /*for each block */
475     {
476         C2F(curblk).kfun = kf + 1;
477         i = funptr[kf];
478         Blocks[kf].type = funtyp[kf + 1];
479         if (i < 0)
480         {
481             switch (funtyp[kf + 1])
482             {
483                 case 0:
484                     Blocks[kf].funpt = F2C(sciblk);
485                     break;
486                 case 1:
487                     sciprint(_("type 1 function not allowed for scilab blocks\n"));
488                     *ierr = 1000 + kf + 1;
489                     FREE_blocks();
490                     return 0;
491                 case 2:
492                     sciprint(_("type 2 function not allowed for scilab blocks\n"));
493                     *ierr = 1000 + kf + 1;
494                     FREE_blocks();
495                     return 0;
496                 case 3:
497                     Blocks[kf].funpt = sciblk2;
498                     Blocks[kf].type = 2;
499                     break;
500                 case 5:
501                     Blocks[kf].funpt = sciblk4;
502                     Blocks[kf].type = 4;
503                     break;
504                 case 99: /* debugging block */
505                     Blocks[kf].funpt = sciblk4;
506                     /*Blocks[kf].type=4;*/
507                     debug_block = kf;
508                     break;
509
510                 case 10005:
511                     Blocks[kf].funpt = sciblk4;
512                     Blocks[kf].type = 10004;
513                     break;
514                 default :
515                     sciprint(_("Undefined Function type\n"));
516                     *ierr = 1000 + kf + 1;
517                     FREE_blocks();
518                     return 0;
519             }
520             Blocks[kf].scsptr = -i; /* set scilab function adress for sciblk */
521         }
522         else if (i <= ntabsim)
523         {
524             Blocks[kf].funpt = *(tabsim[i - 1].fonc);
525             Blocks[kf].scsptr = 0;     /* this is done for being able to test if a block
526                                                                          is a scilab block in the debugging phase when
527                                                                          sciblk4 is called */
528         }
529         else
530         {
531             i -= (ntabsim + 1);
532             GetDynFunc(i, &Blocks[kf].funpt);
533             if ( Blocks[kf].funpt == (voidf) 0)
534             {
535                 sciprint(_("Function not found\n"));
536                 *ierr = 1000 + kf + 1;
537                 FREE_blocks();
538                 return 0;
539             }
540             Blocks[kf].scsptr = 0;   /* this is done for being able to test if a block
541                                                                    is a scilab block in the debugging phase when
542                                                                    sciblk4 is called */
543         }
544
545         /* 2 : Dimension properties */
546         Blocks[kf].ztyp = ztyp[kf + 1];
547         Blocks[kf].nx = xptr[kf + 2] - xptr[kf + 1]; /* continuuous state dimension*/
548         Blocks[kf].ng = zcptr[kf + 2] - zcptr[kf + 1]; /* number of zero crossing surface*/
549         Blocks[kf].nz = zptr[kf + 2] - zptr[kf + 1]; /* number of double discrete state*/
550         Blocks[kf].noz = ozptr[kf + 2] - ozptr[kf + 1]; /* number of other discrete state*/
551         Blocks[kf].nrpar = rpptr[kf + 2] - rpptr[kf + 1]; /* size of double precision parameter vector*/
552         Blocks[kf].nipar = ipptr[kf + 2] - ipptr[kf + 1]; /* size of integer precision parameter vector*/
553         Blocks[kf].nopar = opptr[kf + 2] - opptr[kf + 1]; /* number of other parameters (matrix, data structure,..)*/
554         Blocks[kf].nin = inpptr[kf + 2] - inpptr[kf + 1]; /* number of input ports */
555         Blocks[kf].nout = outptr[kf + 2] - outptr[kf + 1]; /* number of output ports */
556
557         /* 3 : input port properties */
558         /* in insz, we store :
559         *  - insz[0..nin-1] : first dimension of input ports
560         *  - insz[nin..2*nin-1] : second dimension of input ports
561         *  - insz[2*nin..3*nin-1] : type of data of input ports
562         */
563         /* allocate size and pointer arrays (number of input ports)*/
564         Blocks[kf].insz = NULL;
565         Blocks[kf].inptr = NULL;
566         if (Blocks[kf].nin != 0)
567         {
568             if ((Blocks[kf].insz = MALLOC(Blocks[kf].nin * 3 * sizeof(int))) == NULL )
569             {
570                 FREE_blocks();
571                 *ierr = 5;
572                 return 0;
573             }
574             if ((Blocks[kf].inptr = MALLOC(Blocks[kf].nin * sizeof(double*))) == NULL )
575             {
576                 FREE_blocks();
577                 *ierr = 5;
578                 return 0;
579             }
580         }
581         for (in = 0; in < Blocks[kf].nin; in++)
582         {
583             lprt = inplnk[inpptr[kf + 1] + in];
584             Blocks[kf].inptr[in] = outtbptr[lprt - 1]; /* pointer on the data*/
585             Blocks[kf].insz[in] = outtbsz[lprt - 1]; /* row dimension of the input port*/
586             Blocks[kf].insz[Blocks[kf].nin + in] = outtbsz[(lprt - 1) + nlnk]; /* column dimension of the input port*/
587             Blocks[kf].insz[2 * Blocks[kf].nin + in] = outtbtyp[lprt - 1]; /*type of data of the input port*/
588         }
589         /* 4 : output port properties */
590         /* in outsz, we store :
591         *  - outsz[0..nout-1] : first dimension of output ports
592         *  - outsz[nout..2*nout-1] : second dimension of output ports
593         *  - outsz[2*nout..3*nout-1] : type of data of output ports
594         */
595         /* allocate size and pointer arrays (number of output ports)*/
596         Blocks[kf].outsz = NULL;
597         Blocks[kf].outptr = NULL;
598         if (Blocks[kf].nout != 0)
599         {
600             if ((Blocks[kf].outsz = MALLOC(Blocks[kf].nout * 3 * sizeof(int))) == NULL )
601             {
602                 FREE_blocks();
603                 *ierr = 5;
604                 return 0;
605             }
606             if ((Blocks[kf].outptr = MALLOC(Blocks[kf].nout * sizeof(double*))) == NULL )
607             {
608                 FREE_blocks();
609                 *ierr = 5;
610                 return 0;
611             }
612         }
613         /* set the values */
614         for (out = 0; out < Blocks[kf].nout; out++) /*for each output port */
615         {
616             lprt = outlnk[outptr[kf + 1] + out];
617             Blocks[kf].outptr[out] = outtbptr[lprt - 1]; /*pointer on data */
618             Blocks[kf].outsz[out] = outtbsz[lprt - 1]; /*row dimension of output port*/
619             Blocks[kf].outsz[Blocks[kf].nout + out] = outtbsz[(lprt - 1) + nlnk]; /*column dimension of output ports*/
620             Blocks[kf].outsz[2 * Blocks[kf].nout + out] = outtbtyp[lprt - 1]; /*type of data of output port */
621         }
622
623         /* 5 : event output port properties */
624         Blocks[kf].evout = NULL;
625         Blocks[kf].nevout = clkptr[kf + 2] - clkptr[kf + 1];
626         if (Blocks[kf].nevout != 0)
627         {
628             if ((Blocks[kf].evout = CALLOC(Blocks[kf].nevout, sizeof(double))) == NULL )
629             {
630                 FREE_blocks();
631                 *ierr = 5;
632                 return 0;
633             }
634         }
635
636         /* 6 : pointer on the begining of the double discrete state array ( z) */
637         Blocks[kf].z = &(z__[zptr[kf + 1] - 1]);
638
639         /* 7 : type, size and pointer on the other discrete states  data structures (oz) */
640         Blocks[kf].ozsz = NULL;
641         if (Blocks[kf].noz == 0)
642         {
643             Blocks[kf].ozptr = NULL;
644             Blocks[kf].oztyp = NULL;
645         }
646         else
647         {
648             Blocks[kf].ozptr = &(oz[ozptr[kf + 1] - 1]);
649             if ((Blocks[kf].ozsz = MALLOC(Blocks[kf].noz * 2 * sizeof(int))) == NULL )
650             {
651                 FREE_blocks();
652                 *ierr = 5;
653                 return 0;
654             }
655             for (i = 0; i < Blocks[kf].noz; i++)
656             {
657                 Blocks[kf].ozsz[i] = ozsz[(ozptr[kf + 1] - 1) + i];
658                 Blocks[kf].ozsz[i + Blocks[kf].noz] = ozsz[(ozptr[kf + 1] - 1 + noz) + i];
659             }
660             Blocks[kf].oztyp = &(oztyp[ozptr[kf + 1] - 1]);
661         }
662
663         /* 8 : pointer on the begining of the double parameter array ( rpar ) */
664         Blocks[kf].rpar = &(rpar[rpptr[kf + 1] - 1]);
665
666         /* 9 : pointer on the begining of the integer parameter array ( ipar ) */
667         Blocks[kf].ipar = &(ipar[ipptr[kf + 1] - 1]);
668
669         /* 10 : type, size and pointer on the other parameters  data structures (opar) */
670         Blocks[kf].oparsz = NULL;
671         if (Blocks[kf].nopar == 0)
672         {
673             Blocks[kf].oparptr = NULL;
674             Blocks[kf].opartyp = NULL;
675         }
676         else
677         {
678             Blocks[kf].oparptr = &(opar[opptr[kf + 1] - 1]);
679             if ((Blocks[kf].oparsz = MALLOC(Blocks[kf].nopar * 2 * sizeof(int))) == NULL)
680             {
681                 FREE_blocks();
682                 *ierr = 5;
683                 return 0;
684             }
685             for (i = 0; i < Blocks[kf].nopar; i++)
686             {
687                 Blocks[kf].oparsz[i] = oparsz[(opptr[kf + 1] - 1) + i];
688                 Blocks[kf].oparsz[i + Blocks[kf].nopar] = oparsz[(opptr[kf + 1] - 1 + nopar) + i];
689             }
690             Blocks[kf].opartyp = &(opartyp[opptr[kf + 1] - 1]);
691         }
692
693         /* 10 : pointer on the beginning of the residual array (res) */
694         Blocks[kf].res = NULL;
695         if (Blocks[kf].nx != 0)
696         {
697             if ((Blocks[kf].res = MALLOC(Blocks[kf].nx * sizeof(double))) == NULL)
698             {
699                 FREE_blocks();
700                 *ierr = 5;
701                 return 0;
702             }
703         }
704
705         /* 11 : block label (label) */
706         i1 = izptr[kf + 2] - izptr[kf + 1];
707         if ((Blocks[kf].label = MALLOC(sizeof(char) * (i1 + 1))) == NULL)
708         {
709             FREE_blocks();
710             *ierr = 5;
711             return 0;
712         }
713         Blocks[kf].label[i1] = '\0';
714         C2F(cvstr)(&i1, &(iz[izptr[kf + 1] - 1]), Blocks[kf].label, &job, i1);
715
716         /* 12 : block array of crossed surfaces (jroot) */
717         Blocks[kf].jroot = NULL;
718         if (Blocks[kf].ng != 0)
719         {
720             if ((Blocks[kf].jroot = CALLOC(Blocks[kf].ng, sizeof(int))) == NULL)
721             {
722                 FREE_blocks();
723                 *ierr = 5;
724                 return 0;
725             }
726         }
727
728         /* 13 : block work array  (work) */
729         Blocks[kf].work = (void **)(((double *)work) + kf);
730
731         /* 14 : block modes  array  (mode) */
732         Blocks[kf].nmode = modptr[kf + 2] - modptr[kf + 1];
733         if (Blocks[kf].nmode != 0)
734         {
735             Blocks[kf].mode = &(mod[modptr[kf + 1] - 1]);
736         }
737
738         /* 15 : block xprop  array  (xprop) */
739         Blocks[kf].xprop = NULL;
740         if (Blocks[kf].nx != 0)
741         {
742             Blocks[kf].xprop = &(xprop[xptr[kf + 1] - 1]);
743         }
744
745         /* 16 : pointer on the zero crossing surface computation function of the block (g) */
746         Blocks[kf].g = NULL;
747         if (Blocks[kf].ng != 0)
748         {
749             Blocks[kf].g = &(g[zcptr[kf + 1] - 1]);
750         }
751     }
752     /** all block properties are stored in the Blocks array **/
753
754     /* iwa */
755     iwa = NULL;
756     if ((*nevts) != 0)
757     {
758         if ((iwa = MALLOC(sizeof(int) * (*nevts))) == NULL)
759         {
760             FREE_blocks();
761             *ierr = 5;
762             return 0;
763         }
764     }
765
766     /* save ptr of scicos in import structure */
767     C2F(makescicosimport)(x, &nx, &xptr[1], &zcptr[1], z__, &nz, &zptr[1],
768                           &noz, oz, ozsz, oztyp, &ozptr[1],
769                           g, &ng, mod, &nmod, &modptr[1], iz, &izptr[1],
770                           &inpptr[1], &inplnk[1], &outptr[1], &outlnk[1],
771                           outtbptr, outtbsz, outtbtyp,
772                           outtb_elem, &nelem,
773                           &nlnk, rpar, &rpptr[1], ipar, &ipptr[1],
774                           opar, oparsz, opartyp, &opptr[1],
775                           &nblk, subscr, nsubs,
776                           &tevts[1], &evtspt[1], nevts, pointi,
777                           &iord[1], &niord, &oord[1], &noord, &zord[1], &nzord,
778                           funptr, &funtyp[1], &ztyp[1],
779                           &cord[1], &ncord, &ordclk[1], &nordclk, &clkptr[1],
780                           &ordptr[1], &nordptr, &critev[1], iwa, Blocks,
781                           t0, tf, &Atol, &rtol, &ttol, &deltat, &hmax,
782                           xprop, xd);
783
784     if (*flag__ == 1)   /*start*/
785     {
786         /*     initialisation des blocks */
787         for (kf = 0; kf < nblk; ++kf)
788         {
789             *(Blocks[kf].work) = NULL;
790         }
791         cosini(t0);
792         if (*ierr != 0)
793         {
794             ierr0 = *ierr;
795             kfun0 = C2F(curblk).kfun;
796             cosend(t0);
797             *ierr = ierr0;
798             C2F(curblk).kfun = kfun0;
799         }
800
801     }
802     else if (*flag__ == 2)     /*run*/
803     {
804
805         /*     integration */
806         if (C2F(cmsolver).solver == 0)        /*  LSODAR: Method: DYNAMIC, Nonlinear solver= DYNAMIC */
807         {
808             cossim(t0);
809         }
810         else if (C2F(cmsolver).solver == 1)   /*  CVODE: Method: BDF,   Nonlinear solver= NEWTON     */
811         {
812             cossim(t0);
813         }
814         else if (C2F(cmsolver).solver == 2)   /*  CVODE: Method: BDF,   Nonlinear solver= FUNCTIONAL */
815         {
816             cossim(t0);
817         }
818         else if (C2F(cmsolver).solver == 3)   /*  CVODE: Method: ADAMS, Nonlinear solver= NEWTON     */
819         {
820             cossim(t0);
821         }
822         else if (C2F(cmsolver).solver == 4)   /*  CVODE: Method: ADAMS, Nonlinear solver= FUNCTIONAL */
823         {
824             cossim(t0);
825         }
826         else if (C2F(cmsolver).solver == 5)   /*  DOPRI: Method: Dormand-Price, Nonlinear solver=  */
827         {
828             cossim(t0);
829         }
830         else if (C2F(cmsolver).solver == 6)   /*  RK45: Method: Runge-Kutta, Nonlinear solver=  */
831         {
832             cossim(t0);
833         }
834         else if (C2F(cmsolver).solver == 100)  /* IDA  : Method:       , Nonlinear solver=  */
835         {
836             cossimdaskr(t0);
837         }
838         else
839         {
840             /*     add a warning message please */
841         }
842         if (*ierr != 0)
843         {
844             ierr0 = *ierr;
845             kfun0 = C2F(curblk).kfun;
846             cosend(t0);
847             *ierr = ierr0;
848             C2F(curblk).kfun = kfun0;
849         }
850
851     }
852     else if (*flag__ == 3)     /*finish*/
853     {
854         /*     fermeture des blocks */
855         cosend(t0);
856     }
857     else if (*flag__ == 4)     /*linear*/
858     {
859         phase = 1;
860         idoit(t0);
861         if (*ierr == 0)
862         {
863             if ((W = MALLOC(sizeof(double) * (Max(nx, ng)))) == NULL )
864             {
865                 FREE(iwa);
866                 FREE_blocks();
867                 *ierr = 5;
868                 return 0;
869             }
870
871             /*---------à la place de old simblk--------*/
872             /*  C2F(simblk)(&nx, t0, x, W);  */
873
874             if (ng > 0 && nmod > 0)
875             {
876                 zdoit(t0, x, x + nx, W); /* updating modes as a function of state values; this was necessary in iGUI*/
877             }
878             for (jj = 0; jj < nx; jj++) W[jj] = 0.0;
879             C2F(ierode).iero = 0;
880             *ierr = 0;
881             if (C2F(cmsolver).solver < 100)
882             {
883                 odoit(t0, x, W, W);
884             }
885             else
886             {
887                 odoit(t0, x, x + nx, W);
888             }
889             C2F(ierode).iero = *ierr;
890             /*-----------------------------------------*/
891             for (i = 0; i < nx; ++i)
892             {
893                 x[i] = W[i];
894             }
895             FREE(W);
896         }
897     }
898     else if (*flag__ == 5)     /* initial_KINSOL= "Kinsol" */
899     {
900         C2F(ierode).iero = 0;
901         *ierr = 0;
902         idoit(t0);
903         CallKinsol(t0);
904         *ierr = C2F(ierode).iero;
905     }
906
907
908     FREE(iwa);
909     FREE_blocks();
910
911     C2F(clearscicosimport)();
912     return 0;
913 } /* scicos_ */
914 /*--------------------------------------------------------------------------*/
915 /* check_flag */
916 static int check_flag(void *flagvalue, char *funcname, int opt)
917 {
918     int *errflag = NULL;
919
920     /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
921     if (opt == 0 && flagvalue == NULL)
922     {
923         sciprint(_("\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n"), funcname);
924         return(1);
925     }
926     /* Check if flag < 0 */
927     else if (opt == 1)
928     {
929         errflag = (int *) flagvalue;
930         if (*errflag < 0)
931         {
932             sciprint(_("\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n"),
933                      funcname, *errflag);
934             return(1);
935         }
936     }
937     /* Check if function returned NULL pointer - no memory allocated */
938     else if (opt == 2 && flagvalue == NULL)
939     {
940         sciprint(_("\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n"), funcname);
941         return(1);
942     }
943
944     return(0);
945 } /* check_flag */
946
947 /*--------------------------------------------------------------------------*/
948 static void cosini(double *told)
949 {
950     static scicos_flag flag__ = 0;
951     static int i = 0;
952
953     static int kfune = 0;
954     static int jj = 0;
955
956     SCSREAL_COP *outtbd = NULL;    /*to save double of outtb*/
957     SCSINT8_COP *outtbc = NULL;    /*to save int8 of outtb*/
958     SCSINT16_COP *outtbs = NULL;   /*to save int16 of outtb*/
959     SCSINT32_COP *outtbl = NULL;   /*to save int32 of outtb*/
960     SCSUINT8_COP *outtbuc = NULL;  /*to save unsigned int8 of outtb*/
961     SCSUINT16_COP *outtbus = NULL; /*to save unsigned int16 of outtb*/
962     SCSUINT32_COP *outtbul = NULL; /*to save unsigned int32 of outtb*/
963     int szouttbd = 0;  /*size of arrays*/
964     int szouttbc = 0,  szouttbs = 0,  szouttbl = 0;
965     int szouttbuc = 0, szouttbus = 0, szouttbul = 0;
966     int curouttbd = 0; /*current position in arrays*/
967     int curouttbc = 0,  curouttbs = 0,  curouttbl = 0;
968     int curouttbuc = 0, curouttbus = 0, curouttbul = 0;
969
970     int ii = 0, kk = 0; /*local counters*/
971     int sszz = 0;  /*local size of element of outtb*/
972     /*Allocation of arrays for outtb*/
973     for (ii = 0; ii < nlnk; ii++)
974     {
975         switch (outtbtyp[ii])
976         {
977             case SCSREAL_N    :
978                 szouttbd += outtbsz[ii] * outtbsz[ii + nlnk]; /*double real matrix*/
979                 outtbd = (SCSREAL_COP *) REALLOC (outtbd, szouttbd * sizeof(SCSREAL_COP));
980                 break;
981
982             case SCSCOMPLEX_N :
983                 szouttbd += 2 * outtbsz[ii] * outtbsz[ii + nlnk]; /*double complex matrix*/
984                 outtbd = (SCSCOMPLEX_COP *) REALLOC (outtbd, szouttbd * sizeof(SCSCOMPLEX_COP));
985                 break;
986
987             case SCSINT8_N    :
988                 szouttbc += outtbsz[ii] * outtbsz[ii + nlnk]; /*int8*/
989                 outtbc = (SCSINT8_COP *) REALLOC (outtbc, szouttbc * sizeof(SCSINT8_COP));
990                 break;
991
992             case SCSINT16_N   :
993                 szouttbs += outtbsz[ii] * outtbsz[ii + nlnk]; /*int16*/
994                 outtbs = (SCSINT16_COP *) REALLOC (outtbs, szouttbs * sizeof(SCSINT16_COP));
995                 break;
996
997             case SCSINT32_N   :
998                 szouttbl += outtbsz[ii] * outtbsz[ii + nlnk]; /*int32*/
999                 outtbl = (SCSINT32_COP *) REALLOC (outtbl, szouttbl * sizeof(SCSINT32_COP));
1000                 break;
1001
1002             case SCSUINT8_N   :
1003                 szouttbuc += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint8*/
1004                 outtbuc = (SCSUINT8_COP *) REALLOC (outtbuc, szouttbuc * sizeof(SCSUINT8_COP));
1005                 break;
1006
1007             case SCSUINT16_N  :
1008                 szouttbus += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint16*/
1009                 outtbus = (SCSUINT16_COP *) REALLOC (outtbus, szouttbus * sizeof(SCSUINT16_COP));
1010                 break;
1011
1012             case SCSUINT32_N  :
1013                 szouttbul += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint32*/
1014                 outtbul = (SCSUINT32_COP *) REALLOC (outtbul, szouttbul * sizeof(SCSUINT32_COP));
1015                 break;
1016
1017             default  : /* Add a message here */
1018                 break;
1019         }
1020     }
1021
1022     /* Jacobian*/
1023     AJacobian_block = 0;
1024
1025     /* Function Body */
1026     *ierr = 0;
1027
1028     /*     initialization (flag 4) */
1029     /*     loop on blocks */
1030     C2F(dset)(&ng, &c_b14, g, &c__1);
1031
1032     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1033     {
1034         flag__ = 4;
1035         if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1036         {
1037             Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
1038             Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
1039         }
1040         Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1041         if (funtyp[C2F(curblk).kfun] >= 0)   /* debug_block is not called here */
1042         {
1043             /*callf(told, xd, x, x,g,&flag__);*/
1044             Jacobian_Flag = 0;
1045             callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1046             if (flag__ < 0 && *ierr == 0)
1047             {
1048                 *ierr = 5 - flag__;
1049                 kfune = C2F(curblk).kfun;
1050             }
1051             if ((Jacobian_Flag == 1) && (AJacobian_block == 0)) AJacobian_block = C2F(curblk).kfun;
1052         }
1053     }
1054     if (*ierr != 0)
1055     {
1056         C2F(curblk).kfun = kfune;
1057         freeouttbptr;
1058         return;
1059     }
1060
1061     /*     initialization (flag 6) */
1062     flag__ = 6;
1063     for (jj = 1; jj <= ncord; ++jj)
1064     {
1065         C2F(curblk).kfun = cord[jj];
1066         Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1067         if (funtyp[C2F(curblk).kfun] >= 0)
1068         {
1069             /*callf(told, xd, x, x,g,&flag__);*/
1070             callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1071             if (flag__ < 0)
1072             {
1073                 *ierr = 5 - flag__;
1074                 freeouttbptr;
1075                 return;
1076             }
1077         }
1078     }
1079
1080     /*     point-fix iterations */
1081     flag__ = 6;
1082     for (i = 1; i <= nblk + 1; ++i)   /*for each block*/
1083     {
1084         /*     loop on blocks */
1085         for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1086         {
1087             Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1088             if (funtyp[C2F(curblk).kfun] >= 0)
1089             {
1090                 /*callf(told, xd, x, x,g,&flag__);*/
1091                 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1092                 if (flag__ < 0)
1093                 {
1094                     *ierr = 5 - flag__;
1095                     freeouttbptr;
1096                     return;
1097                 }
1098             }
1099         }
1100
1101         flag__ = 6;
1102         for (jj = 1; jj <= ncord; ++jj)   /*for each continous block*/
1103         {
1104             C2F(curblk).kfun = cord[jj];
1105             if (funtyp[C2F(curblk).kfun] >= 0)
1106             {
1107                 /*callf(told, xd, x, x,g,&flag__);*/
1108                 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1109                 if (flag__ < 0)
1110                 {
1111                     *ierr = 5 - flag__;
1112                     freeouttbptr;
1113                     return;
1114                 }
1115             }
1116         }
1117
1118         /*comparison between outtb and arrays*/
1119         curouttbd = 0;
1120         curouttbc = 0;
1121         curouttbs = 0;
1122         curouttbl = 0;
1123         curouttbuc = 0;
1124         curouttbus = 0;
1125         curouttbul = 0;
1126         for (jj = 0; jj < nlnk; jj++)
1127         {
1128             switch (outtbtyp[jj]) /*for each type of ports*/
1129             {
1130                 case SCSREAL_N    :
1131                     outtbdptr = (SCSREAL_COP *)outtbptr[jj]; /*double real matrix*/
1132                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1133                     for (kk = 0; kk < sszz; kk++)
1134                     {
1135                         if (outtbdptr[kk] != (SCSREAL_COP)outtbd[curouttbd + kk]) goto L30;
1136                     }
1137                     curouttbd += sszz;
1138                     break;
1139
1140                 case SCSCOMPLEX_N :
1141                     outtbdptr = (SCSCOMPLEX_COP *)outtbptr[jj]; /*double complex matrix*/
1142                     sszz = 2 * outtbsz[jj] * outtbsz[jj + nlnk];
1143                     for (kk = 0; kk < sszz; kk++)
1144                     {
1145                         if (outtbdptr[kk] != (SCSCOMPLEX_COP)outtbd[curouttbd + kk]) goto L30;
1146                     }
1147                     curouttbd += sszz;
1148                     break;
1149
1150                 case SCSINT8_N    :
1151                     outtbcptr = (SCSINT8_COP *)outtbptr[jj]; /*int8*/
1152                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1153                     for (kk = 0; kk < sszz; kk++)
1154                     {
1155                         if (outtbcptr[kk] != (SCSINT8_COP)outtbc[curouttbc + kk]) goto L30;
1156                     }
1157                     curouttbc += sszz;
1158                     break;
1159
1160                 case SCSINT16_N   :
1161                     outtbsptr = (SCSINT16_COP *)outtbptr[jj]; /*int16*/
1162                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1163                     for (kk = 0; kk < sszz; kk++)
1164                     {
1165                         if (outtbsptr[kk] != (SCSINT16_COP)outtbs[curouttbs + kk]) goto L30;
1166                     }
1167                     curouttbs += sszz;
1168                     break;
1169
1170                 case SCSINT32_N   :
1171                     outtblptr = (SCSINT32_COP *)outtbptr[jj]; /*int32*/
1172                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1173                     for (kk = 0; kk < sszz; kk++)
1174                     {
1175                         if (outtblptr[kk] != (SCSINT32_COP)outtbl[curouttbl + kk]) goto L30;
1176                     }
1177                     curouttbl += sszz;
1178                     break;
1179
1180                 case SCSUINT8_N   :
1181                     outtbucptr = (SCSUINT8_COP *)outtbptr[jj]; /*uint8*/
1182                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1183                     for (kk = 0; kk < sszz; kk++)
1184                     {
1185                         if (outtbucptr[kk] != (SCSUINT8_COP)outtbuc[curouttbuc + kk]) goto L30;
1186                     }
1187                     curouttbuc += sszz;
1188                     break;
1189
1190                 case SCSUINT16_N  :
1191                     outtbusptr = (SCSUINT16_COP *)outtbptr[jj]; /*uint16*/
1192                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1193                     for (kk = 0; kk < sszz; kk++)
1194                     {
1195                         if (outtbusptr[kk] != (SCSUINT16_COP)outtbus[curouttbus + kk]) goto L30;
1196                     }
1197                     curouttbus += sszz;
1198                     break;
1199
1200                 case SCSUINT32_N  :
1201                     outtbulptr = (SCSUINT32_COP *)outtbptr[jj]; /*uint32*/
1202                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1203                     for (kk = 0; kk < sszz; kk++)
1204                     {
1205                         if (outtbulptr[kk] != (SCSUINT32_COP)outtbul[curouttbul + kk]) goto L30;
1206                     }
1207                     curouttbul += sszz;
1208                     break;
1209
1210                 default  : /* Add a message here */
1211                     break;
1212             }
1213         }
1214         freeouttbptr;
1215         return;
1216
1217 L30:
1218         /*Save data of outtb in arrays*/
1219         curouttbd = 0;
1220         curouttbc = 0;
1221         curouttbs = 0;
1222         curouttbl = 0;
1223         curouttbuc = 0;
1224         curouttbus = 0;
1225         curouttbul = 0;
1226         for (ii = 0; ii < nlnk; ii++) /*for each link*/
1227         {
1228             switch (outtbtyp[ii])  /*switch to type of outtb object*/
1229             {
1230                 case SCSREAL_N    :
1231                     outtbdptr = (SCSREAL_COP *)outtbptr[ii]; /*double real matrix*/
1232                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1233                     C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1234                     curouttbd += sszz;
1235                     break;
1236
1237                 case SCSCOMPLEX_N :
1238                     outtbdptr = (SCSCOMPLEX_COP *)outtbptr[ii]; /*double complex matrix*/
1239                     sszz = 2 * outtbsz[ii] * outtbsz[ii + nlnk];
1240                     C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1241                     curouttbd += sszz;
1242                     break;
1243
1244                 case SCSINT8_N    :
1245                     outtbcptr = (SCSINT8_COP *)outtbptr[ii];  /*int8*/
1246                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1247                     for (kk = 0; kk < sszz; kk++) outtbc[curouttbc + kk] = (SCSINT8_COP)outtbcptr[kk];
1248                     curouttbc += sszz;
1249                     break;
1250
1251                 case SCSINT16_N   :
1252                     outtbsptr = (SCSINT16_COP *)outtbptr[ii]; /*int16*/
1253                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1254                     for (kk = 0; kk < sszz; kk++) outtbs[curouttbs + kk] = (SCSINT16_COP)outtbsptr[kk];
1255                     curouttbs += sszz;
1256                     break;
1257
1258                 case SCSINT32_N   :
1259                     outtblptr = (SCSINT32_COP *)outtbptr[ii];  /*int32*/
1260                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1261                     for (kk = 0; kk < sszz; kk++) outtbl[curouttbl + kk] = (SCSINT32_COP)outtblptr[kk];
1262                     curouttbl += sszz;
1263                     break;
1264
1265                 case SCSUINT8_N   :
1266                     outtbucptr = (SCSUINT8_COP *)outtbptr[ii]; /*uint8*/
1267                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1268                     for (kk = 0; kk < sszz; kk++) outtbuc[curouttbuc + kk] = (SCSUINT8_COP)outtbucptr[kk];
1269                     curouttbuc += sszz;
1270                     break;
1271
1272                 case SCSUINT16_N  :
1273                     outtbusptr = (SCSUINT16_COP *)outtbptr[ii]; /*uint16*/
1274                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1275                     for (kk = 0; kk < sszz; kk++) outtbus[curouttbus + kk] = (SCSUINT16_COP)outtbusptr[kk];
1276                     curouttbus += sszz;
1277                     break;
1278
1279                 case SCSUINT32_N  :
1280                     outtbulptr = (SCSUINT32_COP *)outtbptr[ii]; /*uint32*/
1281                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1282                     for (kk = 0; kk < sszz; kk++) outtbul[curouttbul + kk] = (SCSUINT32_COP)outtbulptr[kk];
1283                     curouttbul += sszz;
1284                     break;
1285
1286                 default  : /* Add a message here */
1287                     break;
1288             }
1289         }
1290     }
1291     *ierr = 20;
1292     freeouttbptr;
1293 } /* cosini_ */
1294
1295 /*--------------------------------------------------------------------------*/
1296 static void cossim(double *told)
1297 {
1298     /* System generated locals */
1299     int i3 = 0;
1300
1301     //** used for the [stop] button
1302     static char CommandToUnstack[1024];
1303     static int CommandLength = 0;
1304     static int SeqSync = 0;
1305     static int one = 1;
1306
1307     /* Local variables */
1308     static scicos_flag flag__ = 0;
1309     static int ierr1 = 0;
1310     static int j = 0, k = 0;
1311     static double t = 0.;
1312     static int jj = 0;
1313     static double rhotmp = 0., tstop = 0.;
1314     static int inxsci = 0;
1315     static int kpo = 0, kev = 0;
1316     int Discrete_Jump = 0;
1317     int *jroot = NULL, *zcros = NULL;
1318     realtype reltol = 0., abstol = 0.;
1319     N_Vector y = NULL;
1320     void *cvode_mem = NULL;
1321     int flag = 0, flagr = 0;
1322     int cnt = 0;
1323     jroot = NULL;
1324     if (ng != 0)
1325     {
1326         if ((jroot = MALLOC(sizeof(int) * ng)) == NULL )
1327         {
1328             *ierr = 10000;
1329             return;
1330         }
1331     }
1332
1333     for ( jj = 0 ; jj < ng ; jj++ )
1334         jroot[jj] = 0 ;
1335
1336     zcros = NULL;
1337     if (ng != 0)
1338     {
1339         if ((zcros = MALLOC(sizeof(int) * ng)) == NULL )
1340         {
1341             *ierr = 10000;
1342             if (ng > 0) FREE(jroot);
1343             return;
1344         }
1345     }
1346
1347     reltol = (realtype) rtol;
1348     abstol = (realtype) Atol;  /* Ith(abstol,1) = realtype) Atol;*/
1349
1350     if (*neq > 0) /* Unfortunately CVODE does not work with NEQ==0 */
1351     {
1352         y = N_VNewEmpty_Serial(*neq);
1353         if (check_flag((void *)y, "N_VNewEmpty_Serial", 0))
1354         {
1355             *ierr = 10000;
1356             if (ng > 0) FREE(jroot);
1357             if (ng > 0) FREE(zcros);
1358             return;
1359         }
1360
1361         NV_DATA_S(y) = x;
1362
1363         cvode_mem = NULL;
1364
1365         /* Set extension of Sundials for scicos */
1366         set_sundials_with_extension(TRUE);
1367
1368         switch (C2F(cmsolver).solver)
1369         {
1370             case 0:
1371                 cvode_mem = LSodarCreate(neq, ng); /* Create the lsodar problem */
1372                 break;
1373             case 1:
1374                 cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
1375                 break;
1376             case 2:
1377                 cvode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);
1378                 break;
1379             case 3:
1380                 cvode_mem = CVodeCreate(CV_ADAMS, CV_NEWTON);
1381                 break;
1382             case 4:
1383                 cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
1384                 break;
1385             case 5:
1386                 cvode_mem = CVodeCreate(CV_DOPRI, CV_FUNCTIONAL);
1387                 break;
1388             case 6:
1389                 cvode_mem = CVodeCreate(CV_ExpRK, CV_FUNCTIONAL);
1390                 break;
1391         }
1392
1393         /*    cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);*/
1394
1395         if (check_flag((void *)cvode_mem, "CVodeCreate", 0))
1396         {
1397             *ierr = 10000;
1398             N_VDestroy_Serial(y);
1399             FREE(jroot);
1400             FREE(zcros);
1401             return;
1402         }
1403
1404         if (!C2F(cmsolver).solver)
1405         {
1406             flag = LSodarMalloc(cvode_mem, simblklsodar, T0, y, CV_SS, reltol, &abstol);
1407         }
1408         else
1409         flag = CVodeMalloc(cvode_mem, simblk, T0, y, CV_SS, reltol, &abstol);
1410         if (check_flag(&flag, "CVodeMalloc", 1))
1411         {
1412             *ierr = 300 + (-flag);
1413             freeall
1414             return;
1415         }
1416
1417         if (!C2F(cmsolver).solver)
1418         {
1419             flag = LSodarRootInit(cvode_mem, ng, grblklsodar, NULL);
1420         }
1421         else
1422         flag = CVodeRootInit(cvode_mem, ng, grblk, NULL);
1423         if (check_flag(&flag, "CVodeRootInit", 1))
1424         {
1425             *ierr = 300 + (-flag);
1426             freeall
1427             return;
1428         }
1429
1430         if (C2F(cmsolver).solver)
1431         /* Call CVDense to specify the CVDENSE dense linear solver */
1432         flag = CVDense(cvode_mem, *neq);
1433         if (check_flag(&flag, "CVDense", 1))
1434         {
1435             *ierr = 300 + (-flag);
1436             freeall
1437             return;
1438         }
1439
1440         if (hmax > 0)
1441         {
1442             if (!C2F(cmsolver).solver)
1443             {
1444                 flag = LSodarSetMaxStep(cvode_mem, (realtype) hmax);
1445             }
1446             else
1447             flag = CVodeSetMaxStep(cvode_mem, (realtype) hmax);
1448             if (check_flag(&flag, "CVodeSetMaxStep", 1))
1449             {
1450                 *ierr = 300 + (-flag);
1451                 freeall;
1452                 return;
1453             }
1454         }
1455         /* Set the Jacobian routine to Jac (user-supplied)
1456         flag = CVDenseSetJacFn(cvode_mem, Jac, NULL);
1457         if (check_flag(&flag, "CVDenseSetJacFn", 1)) return(1);  */
1458
1459     }/* testing if neq>0 */
1460
1461     /* Function Body */
1462     C2F(coshlt).halt = 0;
1463     *ierr = 0;
1464
1465     C2F(xscion)(&inxsci);
1466     /*     initialization */
1467     C2F(realtimeinit)(told, &C2F(rtfactor).scale);
1468
1469     phase = 1;
1470     hot = 0;
1471
1472     jj = 0;
1473     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1474     {
1475         if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
1476         {
1477             zcros[jj] = C2F(curblk).kfun;
1478             ++jj;
1479         }
1480     }
1481     /*     . Il faut:  ng >= jj */
1482     if (jj != ng)
1483     {
1484         zcros[jj] = -1;
1485     }
1486     /*     initialisation (propagation of constant blocks outputs) */
1487     idoit(told);
1488     if (*ierr != 0)
1489     {
1490         freeall;
1491         return;
1492     }
1493     /*--discrete zero crossings----dzero--------------------*/
1494     if (ng > 0) /* storing ZC signs just after a solver call*/
1495     {
1496         /*zdoit(told, g, x, x);*/
1497         zdoit(told, x, x, g);
1498         if (*ierr != 0)
1499         {
1500             freeall;
1501             return;
1502         }
1503         for (jj = 0; jj < ng; ++jj)
1504             if (g[jj] >= 0) jroot[jj] = 5;
1505             else jroot[jj] = -5;
1506     }
1507     /*--discrete zero crossings----dzero--------------------*/
1508
1509     /*     main loop on time */
1510     while (*told < *tf)
1511     {
1512         while (ismenu()) //** if the user has done something, do the actions
1513         {
1514             int ierr2 = 0;
1515             SeqSync = GetCommand(CommandToUnstack); //** get at the action
1516             CommandLength = (int)strlen(CommandToUnstack);
1517             syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
1518         }
1519         if (C2F(coshlt).halt != 0)
1520         {
1521             if (C2F(coshlt).halt == 2) *told = *tf; /* end simulation */
1522             C2F(coshlt).halt = 0;
1523             freeall;
1524             return;
1525         }
1526         if (*pointi == 0)
1527         {
1528             t = *tf;
1529         }
1530         else
1531         {
1532             t = tevts[*pointi];
1533         }
1534         if (fabs(t - *told) < ttol)
1535         {
1536             t = *told;
1537             /*     update output part */
1538         }
1539         if (*told > t)
1540         {
1541             /*     !  scheduling problem */
1542             *ierr = 1;
1543             freeall;
1544             return;
1545         }
1546         if (*told != t)
1547         {
1548             if (xptr[nblk + 1] == 1)
1549             {
1550                 /*     .     no continuous state */
1551                 if (*told + deltat + ttol > t)
1552                 {
1553                     *told = t;
1554                 }
1555                 else
1556                 {
1557                     *told += deltat;
1558                 }
1559                 /*     .     update outputs of 'c' type blocks with no continuous state */
1560                 if (*told >= *tf)
1561                 {
1562                     /*     .     we are at the end, update continuous part before leaving */
1563                     if (ncord > 0)
1564                     {
1565                         cdoit(told);
1566                         freeall;
1567                         return;
1568                     }
1569                 }
1570             }
1571             else
1572             {
1573                 /*     integrate */
1574                 rhotmp = *tf + ttol;
1575                 if (*pointi != 0)
1576                 {
1577                     kpo = *pointi;
1578 L20:
1579                     if (critev[kpo] == 1)
1580                     {
1581                         rhotmp = tevts[kpo];
1582                         goto L30;
1583                     }
1584                     kpo = evtspt[kpo];
1585                     if (kpo != 0)
1586                     {
1587                         goto L20;
1588                     }
1589 L30:
1590                     if (rhotmp < tstop)
1591                     {
1592                         hot = 0;
1593                     }
1594                 }
1595                 tstop = rhotmp;
1596                 t = Min(*told + deltat, Min(t, *tf + ttol));
1597
1598                 if (ng > 0 &&  hot == 0 && nmod > 0)
1599                 {
1600                     zdoit(told, x, x, g);
1601                     if (*ierr != 0)
1602                     {
1603                         freeall;
1604                         return;
1605                     }
1606                 }
1607
1608                 if (hot == 0) /* hot==0 : cold restart*/
1609                 {
1610                     if (!C2F(cmsolver).solver)
1611                     {
1612                         flag = LSodarSetStopTime(cvode_mem, (realtype)tstop);
1613                     }
1614                     else
1615                     flag = CVodeSetStopTime(cvode_mem, (realtype)tstop);  /* Setting the stop time*/
1616                     if (check_flag(&flag, "CVodeSetStopTime", 1))
1617                     {
1618                         *ierr = 300 + (-flag);
1619                         freeall;
1620                         return;
1621                     }
1622
1623                     if (!C2F(cmsolver).solver)
1624                     {
1625                         flag = LSodarReInit(cvode_mem, simblklsodar, (realtype)(*told), y, CV_SS, reltol, &abstol);
1626                     }
1627                     else
1628                     flag = CVodeReInit(cvode_mem, simblk, (realtype)(*told), y, CV_SS, reltol, &abstol);
1629                     if (check_flag(&flag, "CVodeReInit", 1))
1630                     {
1631                         *ierr = 300 + (-flag);
1632                         freeall;
1633                         return;
1634                     }
1635                 }
1636
1637                 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1638                 {
1639                     sciprint(_("****SUNDIALS.Cvode from: %f to %f hot= %d  \n"), *told, t, hot);
1640                 }
1641
1642                 /*--discrete zero crossings----dzero--------------------*/
1643                 /*--check for Dzeros after Mode settings or ddoit()----*/
1644                 Discrete_Jump = 0;
1645
1646                 if (ng > 0 && hot == 0)
1647                 {
1648                     zdoit(told, x, x, g);
1649                     if (*ierr != 0)
1650                     {
1651                         freeall;
1652                         return;
1653                     }
1654                     for (jj = 0; jj < ng; ++jj)
1655                     {
1656                         if ((g[jj] >= 0.0) && (jroot[jj] == -5))
1657                         {
1658                             Discrete_Jump = 1;
1659                             jroot[jj] = 1;
1660                         }
1661                         else if ((g[jj] < 0.0) && (jroot[jj] == 5))
1662                         {
1663                             Discrete_Jump = 1;
1664                             jroot[jj] = -1;
1665                         }
1666                         else jroot[jj] = 0;
1667                     }
1668                 }
1669                 /*--discrete zero crossings----dzero--------------------*/
1670
1671                 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
1672                 {
1673                     phase = 2;
1674                     if (!C2F(cmsolver).solver)
1675                     {
1676                         flag = LSodar(cvode_mem, t, y, told, LS_NORMAL);
1677                     }
1678                     else
1679                     flag = CVode(cvode_mem, t, y, told, CV_NORMAL_TSTOP);
1680                     if (*ierr != 0)
1681                     {
1682                         freeall;
1683                         return;
1684                     }
1685                     phase = 1;
1686                 }
1687                 else
1688                 {
1689                     flag = CV_ROOT_RETURN; /* in order to handle discrete jumps */
1690                 }
1691
1692                 /*     .     update outputs of 'c' type  blocks if we are at the end*/
1693                 if (*told >= *tf)
1694                 {
1695                     if (ncord > 0)
1696                     {
1697                         cdoit(told);
1698                         freeall;
1699                         return;
1700                     }
1701                 }
1702
1703                 if (flag >= 0)
1704                 {
1705                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1706                         sciprint(_("****SUNDIALS.Cvode reached: %f\n"), *told);
1707                     hot = 1;
1708                     cnt = 0;
1709                 }
1710                 else if ( flag == CV_TOO_MUCH_WORK ||  flag == CV_CONV_FAILURE || flag == CV_ERR_FAILURE)
1711                 {
1712                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1713                         sciprint(_("****SUNDIALS.Cvode: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
1714                     hot = 0;
1715                     cnt++;
1716                     if (cnt > 5)
1717                     {
1718                         *ierr = 300 + (-flag);
1719                         freeall;
1720                         return;
1721                     }
1722                 }
1723                 else
1724                 {
1725                     if (flag < 0) *ierr = 300 + (-flag); /* raising errors due to internal errors, other wise erros due to flagr*/
1726                     freeall;
1727                     return;
1728                 }
1729
1730                 if (flag == CV_ZERO_DETACH_RETURN)
1731                 {
1732                     hot = 0;
1733                 };  /* new feature of sundials, detects zero-detaching */
1734
1735                 if (flag == CV_ROOT_RETURN)
1736                 {
1737                     /*     .        at a least one root has been found */
1738                     hot = 0;
1739                     if (Discrete_Jump == 0)
1740                     {
1741                         if (!C2F(cmsolver).solver)
1742                         {
1743                             flagr = LSodarGetRootInfo(cvode_mem, jroot);
1744                         }
1745                         else
1746                         flagr = CVodeGetRootInfo(cvode_mem, jroot);
1747                         if (check_flag(&flagr, "CVodeGetRootInfo", 1))
1748                         {
1749                             *ierr = 300 + (-flagr);
1750                             freeall;
1751                             return;
1752                         }
1753                     }
1754                     /*     .        at a least one root has been found */
1755                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1756                     {
1757                         sciprint(_("root found at t=: %f\n"), *told);
1758                     }
1759                     /*     .        update outputs affecting ztyp blocks ONLY FOR OLD BLOCKS */
1760                     zdoit(told, x, xd, g);
1761                     if (*ierr != 0)
1762                     {
1763                         freeall;
1764                         return;
1765                     }
1766                     for (jj = 0; jj < ng; ++jj)
1767                     {
1768                         C2F(curblk).kfun = zcros[ jj];
1769                         if (C2F(curblk).kfun == -1)
1770                         {
1771                             break;
1772                         }
1773                         kev = 0;
1774
1775                         for (j = zcptr[C2F(curblk).kfun] - 1 ;
1776                                 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
1777                         {
1778                             if (jroot[j] != 0)
1779                             {
1780                                 kev = 1;
1781                                 break;
1782                             }
1783                         }
1784                         /*   */
1785                         if (kev != 0)
1786                         {
1787                             Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
1788                             if (funtyp[C2F(curblk).kfun] > 0)
1789                             {
1790
1791                                 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
1792                                 {
1793                                     flag__ = 3;
1794                                     if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1795                                     {
1796                                         Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
1797                                         Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
1798                                     }
1799                                     /* call corresponding block to determine output event (kev) */
1800                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
1801                                     /*callf(told, xd, x, x,g,&flag__);*/
1802                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1803                                     if (flag__ < 0)
1804                                     {
1805                                         *ierr = 5 - flag__;
1806                                         freeall;
1807                                         return;
1808                                     }
1809                                     /*     .              update event agenda */
1810                                     for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
1811                                     {
1812                                         if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0.)
1813                                         {
1814                                             i3 = k + clkptr[C2F(curblk).kfun] ;
1815                                             addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
1816                                             if (ierr1 != 0)
1817                                             {
1818                                                 /*     .                       nevts too small */
1819                                                 *ierr = 3;
1820                                                 freeall;
1821                                                 return;
1822                                             }
1823                                         }
1824                                     }
1825                                 }
1826                                 /*     .              update state */
1827                                 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1828                                 {
1829                                     /*     .              call corresponding block to update state */
1830                                     flag__ = 2;
1831                                     Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
1832                                     Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
1833                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
1834                                     /*callf(told, xd, x, x,g,&flag__);*/
1835                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1836
1837                                     if (flag__ < 0)
1838                                     {
1839                                         *ierr = 5 - flag__;
1840                                         freeall;
1841                                         return;
1842                                     }
1843                                 }
1844                             }
1845                         }
1846                     }
1847                 }
1848             }
1849             /*--discrete zero crossings----dzero--------------------*/
1850             if (ng > 0) /* storing ZC signs just after a sundials call*/
1851             {
1852                 zdoit(told, x, x, g);
1853                 if (*ierr != 0)
1854                 {
1855                     freeall;
1856                     return;
1857                 }
1858                 for (jj = 0; jj < ng; ++jj)
1859                 {
1860                     if (g[jj] >= 0)
1861                     {
1862                         jroot[jj] = 5;
1863                     }
1864                     else
1865                     {
1866                         jroot[jj] = -5;
1867                     }
1868                 }
1869             }
1870             /*--discrete zero crossings----dzero--------------------*/
1871
1872             C2F(realtime)(told);
1873         }
1874         else
1875         {
1876             /*     .  t==told */
1877             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1878             {
1879                 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
1880                 for (kev = 0; kev < nblk; kev++)
1881                 {
1882                     if (Blocks[kev].nmode > 0)
1883                     {
1884                         sciprint(_("mode of block %d=%d, "), kev, Blocks[kev].mode[0]);
1885                     }
1886                 }
1887                 sciprint(_("**mod**\n"));
1888             }
1889
1890             ddoit(told);
1891             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1892             {
1893                 sciprint(_("End of activation\n"));
1894             }
1895             if (*ierr != 0)
1896             {
1897                 freeall;
1898                 return;
1899             }
1900
1901         }
1902         /*     end of main loop on time */
1903     }
1904     freeall;
1905 } /* cossim_ */
1906
1907 /*--------------------------------------------------------------------------*/
1908 static void cossimdaskr(double *told)
1909 {
1910     /* System generated locals */
1911     int i3;
1912     //** used for the [stop] button
1913     static char CommandToUnstack[1024];
1914     static int CommandLength = 0;
1915     static int SeqSync = 0;
1916     static int one = 1;
1917
1918     /* Local variables */
1919     static scicos_flag flag__ = 0;
1920     static int ierr1 = 0;
1921     static int j = 0, k = 0;
1922     static double t = 0.;
1923     static int jj = 0, jt = 0;
1924     static double rhotmp = 0., tstop = 0.;
1925     static int inxsci = 0;
1926     static int kpo = 0, kev = 0;
1927
1928     int *jroot = NULL, *zcros = NULL;
1929     int maxord = 0;
1930     int *Mode_save = NULL;
1931     int Mode_change = 0;
1932
1933     int flag = 0, flagr = 0;
1934     N_Vector   yy = NULL, yp = NULL;
1935     realtype reltol = 0., abstol = 0.;
1936     int Discrete_Jump = 0;
1937     N_Vector IDx = NULL;
1938     realtype *scicos_xproperty = NULL;
1939     DenseMat TJacque = NULL;
1940
1941     void *ida_mem = NULL;
1942     UserData data = NULL;
1943     IDAMem copy_IDA_mem = NULL;
1944     int maxnj = 0, maxnit = 0;
1945     /*-------------------- Analytical Jacobian memory allocation ----------*/
1946     int  Jn = 0, Jnx = 0, Jno = 0, Jni = 0, Jactaille = 0;
1947     double uround = 0.;
1948     int cnt = 0, N_iters = 0;
1949     maxord = 5;
1950
1951
1952     /* Set extension of Sundials for scicos */
1953     set_sundials_with_extension(TRUE);
1954
1955     // CI=1.0;   /* for function Get_Jacobian_ci */
1956     jroot = NULL;
1957     if (ng != 0)
1958     {
1959         if ((jroot = MALLOC(sizeof(int) * ng)) == NULL )
1960         {
1961             *ierr = 10000;
1962             return;
1963         }
1964     }
1965     for ( jj = 0 ; jj < ng ; jj++ )  jroot[jj] = 0 ;
1966
1967     zcros = NULL;
1968     if (ng != 0)
1969     {
1970         if ((zcros = MALLOC(sizeof(int) * ng)) == NULL )
1971         {
1972             *ierr = 10000;
1973             if (ng != 0) FREE(jroot);
1974             return;
1975         }
1976     }
1977
1978     Mode_save = NULL;
1979     if (nmod != 0)
1980     {
1981         if ((Mode_save = MALLOC(sizeof(int) * nmod)) == NULL )
1982         {
1983             *ierr = 10000;
1984             if (ng != 0) FREE(jroot);
1985             if (ng != 0) FREE(zcros);
1986             return;
1987         }
1988     }
1989
1990     reltol = (realtype) rtol;
1991     abstol = (realtype) Atol;  /*  Ith(abstol,1) = realtype) Atol;*/
1992
1993     if (*neq > 0)
1994     {
1995         yy = NULL;
1996         yy = N_VNewEmpty_Serial(*neq);
1997         if (check_flag((void *)yy, "N_VNew_Serial", 0))
1998         {
1999             if (ng != 0) FREE(jroot);
2000             if (ng != 0) FREE(zcros);
2001             if (nmod != 0) FREE(Mode_save);
2002         }
2003         NV_DATA_S(yy) = x;
2004
2005         yp = NULL;
2006         yp = N_VNewEmpty_Serial(*neq);
2007         if (check_flag((void *)yp, "N_VNew_Serial", 0))
2008         {
2009             if (*neq > 0) N_VDestroy_Serial(yy);
2010             if (ng != 0)  FREE(jroot);
2011             if (ng != 0) FREE(zcros);
2012             if (nmod != 0) FREE(Mode_save);
2013             return;
2014         }
2015         NV_DATA_S(yp) = xd;
2016
2017         IDx = NULL;
2018         IDx = N_VNew_Serial(*neq);
2019         if (check_flag((void *)IDx, "N_VNew_Serial", 0))
2020         {
2021             *ierr = 10000;
2022             if (*neq > 0) N_VDestroy_Serial(yp);
2023             if (*neq > 0) N_VDestroy_Serial(yy);
2024             if (ng != 0)  FREE(jroot);
2025             if (ng != 0) FREE(zcros);
2026             if (nmod != 0) FREE(Mode_save);
2027             return;
2028         }
2029
2030         /* Call IDACreate and IDAMalloc to initialize IDA memory */
2031         ida_mem = NULL;
2032         ida_mem = IDACreate();
2033         if (check_flag((void *)ida_mem, "IDACreate", 0))
2034         {
2035             if (*neq > 0) N_VDestroy_Serial(IDx);
2036             if (*neq > 0) N_VDestroy_Serial(yp);
2037             if (*neq > 0) N_VDestroy_Serial(yy);
2038             if (ng != 0) FREE(jroot);
2039             if (ng != 0) FREE(zcros);
2040             if (nmod != 0)  FREE(Mode_save);
2041             return;
2042         }
2043         copy_IDA_mem = (IDAMem) ida_mem;
2044
2045         flag = IDAMalloc(ida_mem, simblkdaskr, T0, yy, yp, IDA_SS, reltol, &abstol);
2046         if (check_flag(&flag, "IDAMalloc", 1))
2047         {
2048             *ierr = 200 + (-flag);
2049             if (*neq > 0)IDAFree(&ida_mem);
2050             if (*neq > 0)N_VDestroy_Serial(IDx);
2051             if (*neq > 0) N_VDestroy_Serial(yp);
2052             if (*neq > 0)N_VDestroy_Serial(yy);
2053             if (ng != 0) FREE(jroot);
2054             if (ng != 0) FREE(zcros);
2055             if (nmod != 0) FREE(Mode_save);
2056             return;
2057         }
2058
2059
2060         flag = IDARootInit(ida_mem, ng, grblkdaskr, NULL);
2061         if (check_flag(&flag, "IDARootInit", 1))
2062         {
2063             *ierr = 200 + (-flag);
2064             if (*neq > 0)IDAFree(&ida_mem);
2065             if (*neq > 0)N_VDestroy_Serial(IDx);
2066             if (*neq > 0)N_VDestroy_Serial(yp);
2067             if (*neq > 0)N_VDestroy_Serial(yy);
2068             if (ng != 0) FREE(jroot);
2069             if (ng != 0) FREE(zcros);
2070             if (nmod != 0) FREE(Mode_save);
2071             return;
2072         }
2073
2074
2075         flag = IDADense(ida_mem, *neq);
2076         if (check_flag(&flag, "IDADense", 1))
2077         {
2078             *ierr = 200 + (-flag);
2079             if (*neq > 0)IDAFree(&ida_mem);
2080             if (*neq > 0)N_VDestroy_Serial(IDx);
2081             if (*neq > 0)N_VDestroy_Serial(yp);
2082             if (*neq > 0)N_VDestroy_Serial(yy);
2083             if (ng != 0) FREE(jroot);
2084             if (ng != 0) FREE(zcros);
2085             if (nmod != 0) FREE(Mode_save);
2086             return;
2087         }
2088
2089         data = NULL;
2090         if ((data = (UserData) MALLOC(sizeof(*data))) == NULL)
2091         {
2092             *ierr = 10000;
2093             if (*neq > 0)IDAFree(&ida_mem);
2094             if (*neq > 0)N_VDestroy_Serial(IDx);
2095             if (*neq > 0)N_VDestroy_Serial(yp);
2096             if (*neq > 0)N_VDestroy_Serial(yy);
2097             if (ng != 0) FREE(jroot);
2098             if (ng != 0) FREE(zcros);
2099             if (nmod != 0) FREE(Mode_save);
2100             return;
2101         }
2102         data->ida_mem = ida_mem;
2103         data->ewt   = NULL;
2104         data->iwork = NULL;
2105         data->rwork = NULL;
2106         data->gwork = NULL;
2107
2108         data->ewt   = N_VNew_Serial(*neq);
2109         if (check_flag((void *)data->ewt, "N_VNew_Serial", 0))
2110         {
2111             *ierr = 200 + (-flag);
2112             if (*neq > 0)FREE(data);
2113             if (*neq > 0)IDAFree(&ida_mem);
2114             if (*neq > 0) N_VDestroy_Serial(IDx);
2115             if (*neq > 0) N_VDestroy_Serial(yp);
2116             if (*neq > 0) N_VDestroy_Serial(yy);
2117             if (ng != 0) FREE(jroot);
2118             if (ng != 0) FREE(zcros);
2119             if (nmod != 0) FREE(Mode_save);
2120             return;
2121         }
2122         if ( ng > 0 )
2123         {
2124             if ((data->gwork = (double *) MALLOC(ng * sizeof(double))) == NULL)
2125             {
2126                 if (*neq > 0) N_VDestroy_Serial(data->ewt);
2127                 if (*neq > 0)FREE(data);
2128                 if (*neq > 0)IDAFree(&ida_mem);
2129                 if (*neq > 0)N_VDestroy_Serial(IDx);
2130                 if (*neq > 0)N_VDestroy_Serial(yp);
2131                 if (*neq > 0)N_VDestroy_Serial(yy);
2132                 if (ng != 0) FREE(jroot);
2133                 if (ng != 0) FREE(zcros);
2134                 if (nmod != 0) FREE(Mode_save);
2135                 return;
2136             }
2137         }
2138         /*Jacobian_Flag=0; */
2139         if (AJacobian_block > 0) /* set by the block with A-Jac in flag-4 using Set_Jacobian_flag(1); */
2140         {
2141             Jn = *neq;
2142             Jnx = Blocks[AJacobian_block - 1].nx;
2143             Jno = Blocks[AJacobian_block - 1].nout;
2144             Jni = Blocks[AJacobian_block - 1].nin;
2145         }
2146         else
2147         {
2148             Jn = *neq;
2149             Jnx = 0;
2150             Jno = 0;
2151             Jni = 0;
2152         }
2153         Jactaille = 3 * Jn + (Jn + Jni) * (Jn + Jno) + Jnx * (Jni + 2 * Jn + Jno) + (Jn - Jnx) * (2 * (Jn - Jnx) + Jno + Jni) + 2 * Jni * Jno;
2154
2155         if ((data->rwork = (double *) MALLOC(Jactaille * sizeof(double))) == NULL)
2156         {
2157             if ( ng > 0 ) FREE(data->gwork);
2158             if (*neq > 0)N_VDestroy_Serial(data->ewt);
2159             if (*neq > 0)FREE(data);
2160             if (*neq > 0)IDAFree(&ida_mem);
2161             if (*neq > 0)N_VDestroy_Serial(IDx);
2162             if (*neq > 0)N_VDestroy_Serial(yp);
2163             if (*neq > 0)N_VDestroy_Serial(yy);
2164             if (ng != 0) FREE(jroot);
2165             if (ng != 0) FREE(zcros);
2166             if (nmod != 0) FREE(Mode_save);
2167             *ierr = 10000;
2168             return;
2169         }
2170
2171         flag = IDADenseSetJacFn(ida_mem, Jacobians, data);
2172         if (check_flag(&flag, "IDADenseSetJacFn", 1))
2173         {
2174             *ierr = 200 + (-flag);
2175             freeallx
2176             return;
2177         }
2178
2179         TJacque = (DenseMat) DenseAllocMat(*neq, *neq);
2180
2181         flag = IDASetRdata(ida_mem, data);
2182         if (check_flag(&flag, "IDASetRdata", 1))
2183         {
2184             *ierr = 200 + (-flag);
2185             freeallx
2186             return;
2187         }
2188
2189         if (hmax > 0)
2190         {
2191             flag = IDASetMaxStep(ida_mem, (realtype) hmax);
2192             if (check_flag(&flag, "IDASetMaxStep", 1))
2193             {
2194                 *ierr = 200 + (-flag);
2195                 freeallx
2196                 return;
2197             }
2198         }
2199
2200         maxnj = 100; /* setting the maximum number of Jacobian evaluation during a Newton step */
2201         flag = IDASetMaxNumJacsIC(ida_mem, maxnj);
2202         if (check_flag(&flag, "IDASetMaxNumJacsIC", 1))
2203         {
2204             *ierr = 200 + (-flag);
2205             freeallx
2206             return;
2207         }
2208
2209         maxnit = 10; /* setting the maximum number of Newton iterations in any one attemp to solve CIC */
2210         flag = IDASetMaxNumItersIC(ida_mem, maxnit);
2211         if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
2212         {
2213             *ierr = 200 + (-flag);
2214             freeallx
2215             return;
2216         }
2217
2218         /* setting the maximum number of steps in an integration interval */
2219         flag = IDASetMaxNumSteps(ida_mem, 2000);
2220         if (check_flag(&flag, "IDASetMaxNumSteps", 1))
2221         {
2222             *ierr = 200 + (-flag);
2223             freeallx
2224             return;
2225         }
2226
2227     } /* testing if neq>0 */
2228
2229     uround = 1.0;
2230     do
2231     {
2232         uround = uround * 0.5;
2233     }
2234     while ( 1.0 + uround != 1.0);
2235     uround = uround * 2.0;
2236     SQuround = sqrt(uround);
2237     /* Function Body */
2238
2239     C2F(coshlt).halt = 0;
2240     *ierr = 0;
2241     /*     hot = .false. */
2242     phase = 1;
2243     hot = 0;
2244
2245     jt = 2;
2246
2247     /*      stuck=.false. */
2248     C2F(xscion)(&inxsci);
2249     /*     initialization */
2250     C2F(realtimeinit)(told, &C2F(rtfactor).scale);
2251     /*     ATOL and RTOL are scalars */
2252
2253     jj = 0;
2254     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
2255     {
2256         if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
2257         {
2258             zcros[jj] = C2F(curblk).kfun;
2259             ++jj;
2260         }
2261     }
2262     /*     . Il faut:  ng >= jj */
2263     if (jj != ng)
2264     {
2265         zcros[jj] = -1;
2266     }
2267     /*     initialisation (propagation of constant blocks outputs) */
2268     idoit(told);
2269     if (*ierr != 0)
2270     {
2271         freeallx;
2272         return;
2273     }
2274
2275     /*--discrete zero crossings----dzero--------------------*/
2276     if (ng > 0) /* storing ZC signs just after a solver call*/
2277     {
2278         zdoit(told, x, x, g);
2279         if (*ierr != 0)
2280         {
2281             freeallx;
2282             return;
2283         }
2284         for (jj = 0; jj < ng; ++jj)
2285             if (g[jj] >= 0) jroot[jj] = 5;
2286             else jroot[jj] = -5;
2287     }
2288     /*     main loop on time */
2289     while (*told < *tf)
2290     {
2291         while (ismenu()) //** if the user has done something, do the actions
2292         {
2293             int ierr2 = 0;
2294             SeqSync = GetCommand(CommandToUnstack); //** get at the action
2295             CommandLength = (int)strlen(CommandToUnstack);
2296             syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2297         }
2298         if (C2F(coshlt).halt != 0)
2299         {
2300             if (C2F(coshlt).halt == 2) *told = *tf; /* end simulation */
2301             C2F(coshlt).halt = 0;
2302             freeallx;
2303             return;
2304         }
2305         if (*pointi == 0)
2306         {
2307             t = *tf;
2308         }
2309         else
2310         {
2311             t = tevts[*pointi];
2312         }
2313         if (fabs(t - *told) < ttol)
2314         {
2315             t = *told;
2316             /*     update output part */
2317         }
2318         if (*told > t)
2319         {
2320             /*     !  scheduling problem */
2321             *ierr = 1;
2322             freeallx;
2323             return;
2324         }
2325         if (*told != t)
2326         {
2327             if (xptr[nblk + 1] == 1)
2328             {
2329                 /*     .     no continuous state */
2330                 if (*told + deltat + ttol > t)
2331                 {
2332                     *told = t;
2333                 }
2334                 else
2335                 {
2336                     *told += deltat;
2337                 }
2338                 /*     .     update outputs of 'c' type blocks with no continuous state */
2339                 if (*told >= *tf)
2340                 {
2341                     /*     .     we are at the end, update continuous part before leaving */
2342                     cdoit(told);
2343                     freeallx;
2344                     return;
2345                 }
2346             }
2347             else
2348             {
2349                 rhotmp = *tf + ttol;
2350                 if (*pointi != 0)
2351                 {
2352                     kpo = *pointi;
2353 L20:
2354                     if (critev[kpo] == 1)
2355                     {
2356                         rhotmp = tevts[kpo];
2357                         goto L30;
2358                     }
2359                     kpo = evtspt[kpo];
2360                     if (kpo != 0)
2361                     {
2362                         goto L20;
2363                     }
2364 L30:
2365                     if (rhotmp < tstop)
2366                     {
2367                         hot = 0;/* Do cold-restart the solver:if the new TSTOP isn't beyong the previous one*/
2368                     }
2369                 }
2370                 tstop = rhotmp;
2371                 t = Min(*told + deltat, Min(t, *tf + ttol));
2372
2373                 if (hot == 0)  /* CIC calculation when hot==0 */
2374                 {
2375
2376                     /* Setting the stop time*/
2377                     flag = IDASetStopTime(ida_mem, (realtype)tstop);
2378                     if (check_flag(&flag, "IDASetStopTime", 1))
2379                     {
2380                         *ierr = 200 + (-flag);
2381                         freeallx;
2382                         return;
2383                     }
2384
2385                     if (ng > 0 && nmod > 0)
2386                     {
2387                         phase = 1;
2388                         zdoit(told, x, xd, g);
2389                         if (*ierr != 0)
2390                         {
2391                             freeallx;
2392                             return;
2393                         }
2394                     }
2395
2396                     /*----------ID setting/checking------------*/
2397                     N_VConst(ONE, IDx); /* Initialize id to 1's. */
2398                     scicos_xproperty = NV_DATA_S(IDx);
2399                     reinitdoit(told);
2400                     if (*ierr > 0)
2401                     {
2402                         freeallx;
2403                         return;
2404                     }
2405                     for (jj = 0; jj < *neq; jj++)
2406                     {
2407                         if (xprop[jj] ==  1) scicos_xproperty[jj] = ONE;
2408                         if (xprop[jj] == -1) scicos_xproperty[jj] = ZERO;
2409                     }
2410                     /* CI=0.0;CJ=100.0; // for functions Get_Jacobian_ci and Get_Jacobian_cj
2411                     Jacobians(*neq, (realtype) (*told), yy, yp, bidon, (realtype) CJ, data, TJacque, tempv1, tempv2, tempv3);
2412                     for (jj=0;jj<*neq;jj++){
2413                     Jacque_col=DENSE_COL(TJacque,jj);
2414                     CI=ZERO;
2415                     for (kk=0;kk<*neq;kk++){
2416                     if ((Jacque_col[kk]-Jacque_col[kk]!=0)) {
2417                     CI=-ONE;
2418                     break;
2419                     }else{
2420                     if (Jacque_col[kk]!=0){
2421                     CI=ONE;
2422                     break;
2423                     }
2424                     }
2425                     }
2426                     if (CI>=ZERO){  scicos_xproperty[jj]=CI;}else{fprintf(stderr,"\nWarinng! Xproperties are not match for i=%d!",jj);}
2427                     } */
2428                     /* printf("\n"); for(jj=0;jj<*neq;jj++) { printf("x%d=%g ",jj,scicos_xproperty[jj]); }*/
2429
2430                     flag = IDASetId(ida_mem, IDx);
2431                     if (check_flag(&flag, "IDASetId", 1))
2432                     {
2433                         *ierr = 200 + (-flag);
2434                         freeallx
2435                         return;
2436                     }
2437                     // CI=1.0;  // for function Get_Jacobian_ci
2438                     /*--------------------------------------------*/
2439                     // maxnj=100; /* setting the maximum number of Jacobian evaluation during a Newton step */
2440                     // flag=IDASetMaxNumJacsIC(ida_mem, maxnj);
2441                     // if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
2442                     //   *ierr=200+(-flag);
2443                     //   freeallx;
2444                     //   return;
2445                     // };
2446                     // flag=IDASetLineSearchOffIC(ida_mem,FALSE);  /* (def=false)  */
2447                     // if (check_flag(&flag, "IDASetLineSearchOffIC", 1)) {
2448                     //   *ierr=200+(-flag);
2449                     //   freeallx;
2450                     //   return;
2451                     // };
2452                     // flag=IDASetMaxNumItersIC(ida_mem, 10);/* (def=10) setting the maximum number of Newton iterations in any one attemp to solve CIC */
2453                     // if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
2454                     //   *ierr=200+(-flag);
2455                     //   freeallx;
2456                     //   return;
2457                     // };
2458
2459                     N_iters = 4 + nmod * 4;
2460                     for (j = 0; j <= N_iters; j++)
2461                     {
2462                         /* counter to reevaluate the
2463                                                                         modes in  mode->CIC->mode->CIC-> loop
2464                                                                         do it once in the absence of mode (nmod=0) */
2465                         /* updating the modes through Flag==9, Phase==1 */
2466
2467                         /* Serge Steer 29/06/2009 */
2468                         while (ismenu()) //** if the user has done something, do the actions
2469                         {
2470                             int ierr2 = 0;
2471                             SeqSync = GetCommand(CommandToUnstack); //** get at the action
2472                             CommandLength = (int)strlen(CommandToUnstack);
2473                             syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2474                         }
2475                         if (C2F(coshlt).halt != 0)
2476                         {
2477                             if (C2F(coshlt).halt == 2) *told = *tf; /* end simulation */
2478                             C2F(coshlt).halt = 0;
2479                             freeallx;
2480                             return;
2481                         }
2482
2483                         /* yy->PH */
2484                         flag = IDAReInit(ida_mem, simblkdaskr, (realtype)(*told), yy, yp, IDA_SS, reltol, &abstol);
2485                         if (check_flag(&flag, "CVodeReInit", 1))
2486                         {
2487                             *ierr = 200 + (-flag);
2488                             freeallx;
2489                             return;
2490                         }
2491
2492                         phase = 2; /* IDACalcIC: PHI-> yy0: if (ok) yy0_cic-> PHI*/
2493                         copy_IDA_mem->ida_kk = 1;
2494
2495                         // the initial conditons y0 and yp0 do not satisfy the DAE
2496                         flagr = IDACalcIC(ida_mem, IDA_YA_YDP_INIT, (realtype)(t));
2497
2498                         phase = 1;
2499                         flag = IDAGetConsistentIC(ida_mem, yy, yp); /* PHI->YY */
2500
2501                         if (*ierr > 5)    /* *ierr>5 => singularity in block */
2502                         {
2503                             freeallx;
2504                             return;
2505                         }
2506
2507                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2508                         {
2509                             if (flagr >= 0)
2510                             {
2511                                 sciprint(_("**** SUNDIALS.IDA successfully initialized *****\n") );
2512                             }
2513                             else
2514                             {
2515                                 sciprint(_("**** SUNDIALS.IDA failed to initialize ->try again *****\n") );
2516                             }
2517                         }
2518                         /*-------------------------------------*/
2519                         /* saving the previous modes*/
2520                         for (jj = 0; jj < nmod; ++jj)
2521                         {
2522                             Mode_save[jj] = mod[jj];
2523                         }
2524                         if (ng > 0 && nmod > 0)
2525                         {
2526                             phase = 1;
2527                             zdoit(told, x, xd, g);
2528                             if (*ierr != 0)
2529                             {
2530                                 freeallx;
2531                                 return;
2532                             }
2533                         }
2534                         /*------------------------------------*/
2535                         Mode_change = 0;
2536                         for (jj = 0; jj < nmod; ++jj)
2537                         {
2538                             if (Mode_save[jj] != mod[jj])
2539                             {
2540                                 Mode_change = 1;
2541                                 break;
2542                             }
2543                         }
2544                         if (Mode_change == 0)
2545                         {
2546                             if (flagr >= 0 )
2547                             {
2548                                 break;  /*   if (flagr>=0) break;  else{ *ierr=200+(-flagr); freeallx;  return; }*/
2549                             }
2550                             else if (j >= (int)( N_iters / 2))
2551                             {
2552                                 /* IDASetMaxNumStepsIC(mem,10); */     /* maxnh (def=5) */
2553                                 IDASetMaxNumJacsIC(ida_mem, 10);      /* maxnj 100 (def=4)*/
2554                                 /* IDASetMaxNumItersIC(mem,100000); */ /* maxnit in IDANewtonIC (def=10) */
2555                                 IDASetLineSearchOffIC(ida_mem, TRUE); /* (def=false)  */
2556                                 /* IDASetNonlinConvCoefIC(mem,1.01);*/ /* (def=0.01-0.33*/
2557                                 flag = IDASetMaxNumItersIC(ida_mem, 1000);
2558                                 if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
2559                                 {
2560                                     *ierr = 200 + (-flag);
2561                                     freeallx;
2562                                     return;
2563                                 };
2564                             }
2565                         }
2566                     }/* mode-CIC  counter*/
2567                     if (Mode_change == 1)
2568                     {
2569                         /* In tghis case, we try again by relaxing all modes and calling IDA_calc again
2570                         /Masoud */
2571                         phase = 1;
2572                         copy_IDA_mem->ida_kk = 1;
2573                         flagr = IDACalcIC(ida_mem, IDA_YA_YDP_INIT, (realtype)(t));
2574                         phase = 1;
2575                         flag = IDAGetConsistentIC(ida_mem, yy, yp); /* PHI->YY */
2576                         if ((flagr < 0) || (*ierr > 5)) /* *ierr>5 => singularity in block */
2577                         {
2578                             *ierr = 23;
2579                             freeallx;
2580                             return;
2581                         }
2582                     }
2583                     /*-----If flagr<0 the initialization solver has not converged-----*/
2584                     if (flagr < 0)
2585                     {
2586                         *ierr = 237;
2587                         freeallx;
2588                         return;
2589                     }
2590
2591                 } /* CIC calculation when hot==0 */
2592
2593                 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2594                 {
2595                     sciprint(_("****daskr from: %f to %f hot= %d  \n"), *told, t, hot);
2596                 }
2597
2598                 /*--discrete zero crossings----dzero--------------------*/
2599                 /*--check for Dzeros after Mode settings or ddoit()----*/
2600                 Discrete_Jump = 0;
2601                 if (ng > 0 && hot == 0)
2602                 {
2603                     zdoit(told, x, xd, g);
2604                     if (*ierr != 0)
2605                     {
2606                         freeallx;
2607                         return;
2608                     }
2609                     for (jj = 0; jj < ng; ++jj)
2610                     {
2611                         if ((g[jj] >= 0.0) && ( jroot[jj] == -5))
2612                         {
2613                             Discrete_Jump = 1;
2614                             jroot[jj] = 1;
2615                         }
2616                         else if ((g[jj] < 0.0) && ( jroot[jj] == 5))
2617                         {
2618                             Discrete_Jump = 1;
2619                             jroot[jj] = -1;
2620                         }
2621                         else jroot[jj] = 0;
2622                     }
2623                 }
2624
2625                 /*--discrete zero crossings----dzero--------------------*/
2626                 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
2627                 {
2628                     phase = 2;
2629                     flagr = IDASolve(ida_mem, t, told, yy, yp, IDA_NORMAL_TSTOP);
2630                     phase = 1;
2631                     if (*ierr != 0)
2632                     {
2633                         freeallx;
2634                         return;
2635                     }
2636                 }
2637                 else
2638                 {
2639                     flagr = IDA_ROOT_RETURN; /* in order to handle discrete jumps */
2640                 }
2641                 if (flagr >= 0)
2642                 {
2643                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2644                         sciprint(_("****SUNDIALS.Ida reached: %f\n"), *told);
2645                     hot = 1;
2646                     cnt = 0;
2647                 }
2648                 else if ( flagr == IDA_TOO_MUCH_WORK ||  flagr == IDA_CONV_FAIL || flagr == IDA_ERR_FAIL)
2649                 {
2650                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2651                         sciprint(_("**** SUNDIALS.Ida: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
2652                     hot = 0;
2653                     cnt++;
2654                     if (cnt > 5)
2655                     {
2656                         *ierr = 200 + (-flagr);
2657                         freeallx;
2658                         return;
2659                     }
2660                 }
2661                 else
2662                 {
2663                     if (flagr < 0) *ierr = 200 + (-flagr); /* raising errors due to internal errors, other wise erros due to flagr*/
2664                     freeallx;
2665                     return;
2666                 }
2667
2668                 /*     update outputs of 'c' type  blocks if we are at the end*/
2669                 if (*told >= *tf)
2670                 {
2671                     cdoit(told);
2672                     freeallx;
2673                     return;
2674                 }
2675
2676                 if (flagr == IDA_ZERO_DETACH_RETURN)
2677                 {
2678                     hot = 0;
2679                 }; /* new feature of sundials, detects unmasking */
2680                 if (flagr == IDA_ROOT_RETURN)
2681                 {
2682                     /*     .        at a least one root has been found */
2683                     hot = 0;
2684                     if (Discrete_Jump == 0)
2685                     {
2686                         flagr = IDAGetRootInfo(ida_mem, jroot);
2687                         if (check_flag(&flagr, "IDAGetRootInfo", 1))
2688                         {
2689                             *ierr = 200 + (-flagr);
2690                             freeallx;
2691                             return;
2692                         }
2693                     }
2694
2695                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2696                     {
2697                         sciprint(_("root found at t=: %f\n"), *told);
2698                     }
2699                     /*     .        update outputs affecting ztyp blocks  ONLY FOR OLD BLOCKS*/
2700                     zdoit(told, x, xd, g);
2701                     if (*ierr != 0)
2702                     {
2703                         freeallx;
2704                         return;
2705                     }
2706                     for (jj = 0; jj < ng; ++jj)
2707                     {
2708                         C2F(curblk).kfun = zcros[jj];
2709                         if (C2F(curblk).kfun == -1)
2710                         {
2711                             break;
2712                         }
2713                         kev = 0;
2714                         for (j = zcptr[C2F(curblk).kfun] - 1 ;
2715                                 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
2716                         {
2717                             if (jroot[j] != 0)
2718                             {
2719                                 kev = 1;
2720                                 break;
2721                             }
2722                         }
2723                         if (kev != 0)
2724                         {
2725                             Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
2726                             if (funtyp[C2F(curblk).kfun] > 0)
2727                             {
2728                                 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
2729                                 {
2730                                     flag__ = 3;
2731                                     if (Blocks[C2F(curblk).kfun - 1].nx > 0)
2732                                     {
2733                                         Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
2734                                         Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
2735                                     }
2736                                     /*     call corresponding block to determine output event (kev) */
2737                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
2738                                     /*callf(told, xd, x, x,g,&flag__);*/
2739                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
2740                                     if (flag__ < 0)
2741                                     {
2742                                         *ierr = 5 - flag__;
2743                                         freeallx;
2744                                         return;
2745                                     }
2746                                     /*     update event agenda */
2747                                     for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
2748                                     {
2749                                         if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0)
2750                                         {
2751                                             i3 = k + clkptr[C2F(curblk).kfun] ;
2752                                             addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
2753                                             if (ierr1 != 0)
2754                                             {
2755                                                 /*     .                       nevts too small */
2756                                                 *ierr = 3;
2757                                                 freeallx;
2758                                                 return;
2759                                             }
2760                                         }
2761                                     }
2762                                 }
2763                                 /* update state */
2764                                 if ((Blocks[C2F(curblk).kfun - 1].nx > 0) || (*Blocks[C2F(curblk).kfun - 1].work != NULL) )
2765                                 {
2766                                     /* call corresponding block to update state */
2767                                     flag__ = 2;
2768                                     if (Blocks[C2F(curblk).kfun - 1].nx > 0)
2769                                     {
2770                                         Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
2771                                         Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
2772                                     }
2773                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
2774
2775                                     Blocks[C2F(curblk).kfun - 1].xprop = &xprop[-1 + xptr[C2F(curblk).kfun]];
2776                                     /*callf(told, xd, x, x,g,&flag__);*/
2777                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
2778
2779                                     if (flag__ < 0)
2780                                     {
2781                                         *ierr = 5 - flag__;
2782                                         freeallx;
2783                                         return;
2784                                     }
2785                                     for (j = 0; j < *neq; j++) /* Adjust xprop for IDx */
2786                                     {
2787                                         if (xprop[j] ==  1) scicos_xproperty[j] = ONE;
2788                                         if (xprop[j] == -1) scicos_xproperty[j] = ZERO;
2789                                     }
2790                                 }
2791                             }
2792                         }
2793                     }
2794                 }
2795                 /* Serge Steer 29/06/2009 */
2796                 while (ismenu()) //** if the user has done something, do the actions
2797                 {
2798                     int ierr2 = 0;
2799                     SeqSync = GetCommand(CommandToUnstack); //** get at the action
2800                     CommandLength = (int)strlen(CommandToUnstack);
2801                     syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2802                 }
2803
2804                 if (C2F(coshlt).halt != 0)
2805                 {
2806                     if (C2F(coshlt).halt == 2) *told = *tf; /* end simulation */
2807                     C2F(coshlt).halt = 0;
2808                     freeallx;
2809                     return;
2810                 }
2811                 /* if(*pointi!=0){
2812                 t=tevts[*pointi];
2813                 if(*told<t-ttol){
2814                 cdoit(told);
2815                 goto L15;
2816                 }
2817                 }else{
2818                 if(*told<*tf){
2819                 cdoit(told);
2820                 goto L15;
2821                 }
2822                 }*/
2823
2824                 /*--discrete zero crossings----dzero--------------------*/
2825                 if (ng > 0) /* storing ZC signs just after a ddaskr call*/
2826                 {
2827                     zdoit(told, x, xd, g);
2828                     if (*ierr != 0)
2829                     {
2830                         freeallx;
2831                         return;
2832                     }
2833                     for (jj = 0; jj < ng; ++jj)
2834                     {
2835                         if (g[jj] >= 0)
2836                         {
2837                             jroot[jj] = 5;
2838                         }
2839                         else
2840                         {
2841                             jroot[jj] = -5;
2842                         }
2843                     }
2844                 }
2845                 /*--discrete zero crossings----dzero--------------------*/
2846             }
2847             C2F(realtime)(told);
2848         }
2849         else
2850         {
2851             /*     .  t==told */
2852             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2853             {
2854                 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
2855             }
2856
2857             ddoit(told);
2858             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2859             {
2860                 sciprint(_("End of activation"));
2861             }
2862             if (*ierr != 0)
2863             {
2864                 freeallx;
2865                 return;
2866             }
2867         }
2868         /*     end of main loop on time */
2869     }
2870     freeallx;
2871 } /* cossimdaskr_ */
2872 /*--------------------------------------------------------------------------*/
2873 /* Subroutine cosend */
2874 static void cosend(double *told)
2875 {
2876     /* Local variables */
2877     static scicos_flag flag__ = 0;
2878
2879     static int kfune = 0;
2880
2881     /* Function Body */
2882     *ierr = 0;
2883     /*     loop on blocks */
2884     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
2885     {
2886         flag__ = 5;
2887         Blocks[C2F(curblk).kfun - 1].nevprt = 0;
2888         if (funtyp[C2F(curblk).kfun] >= 0)
2889         {
2890             if (Blocks[C2F(curblk).kfun - 1].nx > 0)
2891             {
2892                 Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
2893                 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
2894             }
2895             /*callf(told, xd, x, x,x,&flag__);*/
2896             callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
2897             if (flag__ < 0 && *ierr == 0)
2898             {
2899                 *ierr = 5 - flag__;
2900                 kfune = C2F(curblk).kfun;
2901             }
2902         }
2903     }
2904     if (*ierr != 0)
2905     {
2906         C2F(curblk).kfun = kfune;
2907         return;
2908     }
2909 } /* cosend_ */
2910 /*--------------------------------------------------------------------------*/
2911 /* callf */
2912 void callf(double *t, scicos_block *block, scicos_flag *flag)
2913 {
2914     double* args[SZ_SIZE];
2915     int sz[SZ_SIZE];
2916     double intabl[TB_SIZE];
2917     double outabl[TB_SIZE];
2918
2919     int ii = 0, in = 0, out = 0, ki = 0, ko = 0, no = 0, ni = 0, k = 0, j = 0;
2920     int szi = 0, flagi = 0;
2921     double *ptr_d = NULL;
2922
2923     /* function pointers type def */
2924     voidf loc ;
2925     ScicosF0 loc0;
2926     ScicosF loc1;
2927     /*  ScicosFm1 loc3;*/
2928     ScicosF2 loc2;
2929     ScicosF2z loc2z;
2930     ScicosFi loci1;
2931     ScicosFi2 loci2;
2932     ScicosFi2z loci2z;
2933     ScicosF4 loc4;
2934
2935     int solver = C2F(cmsolver).solver;
2936     int cosd   = C2F(cosdebug).cosd;
2937     /*int kf     = C2F(curblk).kfun;*/
2938     scicos_time     = *t;
2939     block_error     = (int*) flag;
2940
2941     /* debug block is never called */
2942     /*if (kf==(debug_block+1)) return;*/
2943     if (block->type == 99) return;
2944
2945     /* flag 7 implicit initialization */
2946     flagi = (int) * flag;
2947     /* change flag to zero if flagi==7 for explicit block */
2948     if (flagi == 7 && block->type < 10000)
2949     {
2950         *flag = 0;
2951     }
2952
2953     /* display information for debugging mode */
2954     if (cosd > 1)
2955     {
2956         if (cosd != 3)
2957         {
2958             sciprint(_("block %d is called "), C2F(curblk).kfun);
2959             sciprint(_("with flag %d "), *flag);
2960             sciprint(_("at time %f \n"), *t);
2961         }
2962         if (debug_block > -1)
2963         {
2964             if (cosd != 3) sciprint(_("Entering the block \n"));
2965             call_debug_scicos(block, flag, flagi, debug_block);
2966             if (*flag < 0) return; /* error in debug block */
2967         }
2968     }
2969
2970     C2F(scsptr).ptr = block->scsptr;
2971
2972     /* get pointer of the function */
2973     loc = block->funpt;
2974
2975     /* continuous state */
2976     if (solver == 100 && block->type < 10000 && *flag == 0)
2977     {
2978         ptr_d = block->xd;
2979         block->xd  = block->res;
2980     }
2981
2982     /* switch loop */
2983     //sciprint("callf type=%d flag=%d\n",block->type,flagi);
2984     switch (block->type)
2985     {
2986             /*******************/
2987             /* function type 0 */
2988             /*******************/
2989         case 0 :
2990         {
2991             /* This is for compatibility */
2992             /* jroot is returned in g for old type */
2993             if (block->nevprt < 0)
2994             {
2995                 for (j = 0; j < block->ng; ++j)
2996                 {
2997                     block->g[j] = (double)block->jroot[j];
2998                 }
2999             }
3000
3001             /* concatenated entries and concatened outputs */
3002             /* catenate inputs if necessary */
3003             ni = 0;
3004             if (block->nin > 1)
3005             {
3006                 ki = 0;
3007                 for (in = 0; in < block->nin; in++)
3008                 {
3009                     szi = block->insz[in] * block->insz[in + block->nin];
3010                     for (ii = 0; ii < szi; ii++)
3011                     {
3012                         intabl[ki++] = *((double *)(block->inptr[in]) + ii);
3013                     }
3014                     ni = ni + szi;
3015                 }
3016                 args[0] = &(intabl[0]);
3017             }
3018             else
3019             {
3020                 if (block->nin == 0)
3021                 {
3022                     args[0] = NULL;
3023                 }
3024                 else
3025                 {
3026                     args[0] = (double *)(block->inptr[0]);
3027                     ni = block->insz[0] * block->insz[1];
3028                 }
3029             }
3030
3031             /* catenate outputs if necessary */
3032             no = 0;
3033             if (block->nout > 1)
3034             {
3035                 ko = 0;
3036                 for (out = 0; out < block->nout; out++)
3037                 {
3038                     szi = block->outsz[out] * block->outsz[out + block->nout];
3039                     for (ii = 0; ii < szi; ii++)
3040                     {
3041                         outabl[ko++] = *((double *)(block->outptr[out]) + ii);
3042                     }
3043                     no = no + szi;
3044                 }
3045                 args[1] = &(outabl[0]);
3046             }
3047             else
3048             {
3049                 if (block->nout == 0)
3050                 {
3051                     args[1] = NULL;
3052                 }
3053                 else
3054                 {
3055                     args[1] = (double *)(block->outptr[0]);
3056                     no = block->outsz[0] * block->outsz[1];
3057                 }
3058             }
3059
3060             loc0 = (ScicosF0) loc;
3061
3062             (*loc0)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3063                     block->z, &block->nz,
3064                     block->evout, &block->nevout, block->rpar, &block->nrpar,
3065                     block->ipar, &block->nipar, (double *)args[0], &ni,
3066                     (double *)args[1], &no);
3067
3068             /* split output vector on each port if necessary */
3069             if (block->nout > 1)
3070             {
3071                 ko = 0;
3072                 for (out = 0; out < block->nout; out++)
3073                 {
3074                     szi = block->outsz[out] * block->outsz[out + block->nout];
3075                     for (ii = 0; ii < szi; ii++)
3076                     {
3077                         *((double *)(block->outptr[out]) + ii) = outabl[ko++];
3078                     }
3079                 }
3080             }
3081
3082             /* adjust values of output register */
3083             for (in = 0; in < block->nevout; ++in)
3084             {
3085                 block->evout[in] = block->evout[in] - *t;
3086             }
3087
3088             break;
3089         }
3090
3091         /*******************/
3092         /* function type 1 */
3093         /*******************/
3094         case 1 :
3095         {
3096             /* This is for compatibility */
3097             /* jroot is returned in g for old type */
3098             if (block->nevprt < 0)
3099             {
3100                 for (j = 0; j < block->ng; ++j)
3101                 {
3102                     block->g[j] = (double)block->jroot[j];
3103                 }
3104             }
3105
3106             /* one entry for each input or output */
3107             for (in = 0 ; in < block->nin ; in++)
3108             {
3109                 args[in] = block->inptr[in];
3110                 sz[in] = block->insz[in];
3111             }
3112             for (out = 0; out < block->nout; out++)
3113             {
3114                 args[in + out] = block->outptr[out];
3115                 sz[in + out] = block->outsz[out];
3116             }
3117             /* with zero crossing */
3118             if (block->ztyp > 0)
3119             {
3120                 args[block->nin + block->nout] = block->g;
3121                 sz[block->nin + block->nout] = block->ng;
3122             }
3123
3124             loc1 = (ScicosF) loc;
3125
3126             (*loc1)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3127                     block->z, &block->nz,
3128                     block->evout, &block->nevout, block->rpar, &block->nrpar,
3129                     block->ipar, &block->nipar,
3130                     (double *)args[0], &sz[0],
3131                     (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3132                     (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3133                     (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3134                     (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3135                     (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3136                     (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3137                     (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3138                     (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3139                     (double *)args[17], &sz[17]);
3140
3141             /* adjust values of output register */
3142             for (in = 0; in < block->nevout; ++in)
3143             {
3144                 block->evout[in] = block->evout[in] - *t;
3145             }
3146
3147             break;
3148         }
3149
3150         /*******************/
3151         /* function type 2 */
3152         /*******************/
3153         case 2 :
3154         {
3155             /* This is for compatibility */
3156             /* jroot is returned in g for old type */
3157             if (block->nevprt < 0)
3158             {
3159                 for (j = 0; j < block->ng; ++j)
3160                 {
3161                     block->g[j] = (double)block->jroot[j];
3162                 }
3163             }
3164
3165             /* no zero crossing */
3166             if (block->ztyp == 0)
3167             {
3168                 loc2 = (ScicosF2) loc;
3169                 (*loc2)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3170                         block->z, &block->nz,
3171                         block->evout, &block->nevout, block->rpar, &block->nrpar,
3172                         block->ipar, &block->nipar, (double **)block->inptr,
3173                         block->insz, &block->nin,
3174                         (double **)block->outptr, block->outsz, &block->nout);
3175             }
3176             /* with zero crossing */
3177             else
3178             {
3179                 loc2z = (ScicosF2z) loc;
3180                 (*loc2z)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3181                          block->z, &block->nz,
3182                          block->evout, &block->nevout, block->rpar, &block->nrpar,
3183                          block->ipar, &block->nipar, (double **)block->inptr,
3184                          block->insz, &block->nin,
3185                          (double **)block->outptr, block->outsz, &block->nout,
3186                          block->g, &block->ng);
3187             }
3188
3189             /* adjust values of output register */
3190             for (in = 0; in < block->nevout; ++in)
3191             {
3192                 block->evout[in] = block->evout[in] - *t;
3193             }
3194
3195             break;
3196         }
3197
3198         /*******************/
3199         /* function type 4 */
3200         /*******************/
3201         case 4 :
3202         {
3203             /* get pointer of the function type 4*/
3204             loc4 = (ScicosF4) loc;
3205
3206             (*loc4)(block, *flag);
3207
3208             break;
3209         }
3210
3211         /***********************/
3212         /* function type 10001 */
3213         /***********************/
3214         case 10001 :
3215         {
3216             /* This is for compatibility */
3217             /* jroot is returned in g for old type */
3218             if (block->nevprt < 0)
3219             {
3220                 for (j = 0; j < block->ng; ++j)
3221                 {
3222                     block->g[j] = (double)block->jroot[j];
3223                 }
3224             }
3225
3226             /* implicit block one entry for each input or output */
3227             for (in = 0 ; in < block->nin ; in++)
3228             {
3229                 args[in] = block->inptr[in];
3230                 sz[in] = block->insz[in];
3231             }
3232             for (out = 0; out < block->nout; out++)
3233             {
3234                 args[in + out] = block->outptr[out];
3235                 sz[in + out] = block->outsz[out];
3236             }
3237             /* with zero crossing */
3238             if (block->ztyp > 0)
3239             {
3240                 args[block->nin + block->nout] = block->g;
3241                 sz[block->nin + block->nout] = block->ng;
3242             }
3243
3244             loci1 = (ScicosFi) loc;
3245             (*loci1)(flag, &block->nevprt, t, block->res, block->xd, block->x,
3246                      &block->nx, block->z, &block->nz,
3247                      block->evout, &block->nevout, block->rpar, &block->nrpar,
3248                      block->ipar, &block->nipar,
3249                      (double *)args[0], &sz[0],
3250                      (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3251                      (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3252                      (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3253                      (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3254                      (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3255                      (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3256                      (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3257                      (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3258                      (double *)args[17], &sz[17]);
3259
3260             /* adjust values of output register */
3261             for (in = 0; in < block->nevout; ++in)
3262             {
3263                 block->evout[in] = block->evout[in] - *t;
3264             }
3265
3266             break;
3267         }
3268
3269         /***********************/
3270         /* function type 10002 */
3271         /***********************/
3272         case 10002 :
3273         {
3274             /* This is for compatibility */
3275             /* jroot is returned in g for old type */
3276             if (block->nevprt < 0)
3277             {
3278                 for (j = 0; j < block->ng; ++j)
3279                 {
3280                     block->g[j] = (double)block->jroot[j];
3281                 }
3282             }
3283
3284             /* implicit block, inputs and outputs given by a table of pointers */
3285             /* no zero crossing */
3286             if (block->ztyp == 0)
3287             {
3288                 loci2 = (ScicosFi2) loc;
3289                 (*loci2)(flag, &block->nevprt, t, block->res,
3290                          block->xd, block->x, &block->nx,
3291                          block->z, &block->nz,
3292                          block->evout, &block->nevout, block->rpar, &block->nrpar,
3293                          block->ipar, &block->nipar, (double **)block->inptr,
3294                          block->insz, &block->nin,
3295                          (double **)block->outptr, block->outsz, &block->nout);
3296             }
3297             /* with zero crossing */
3298             else
3299             {
3300                 loci2z = (ScicosFi2z) loc;
3301                 (*loci2z)(flag, &block->nevprt, t, block->res,
3302                           block->xd, block->x, &block->nx,
3303                           block->z, &block->nz,
3304                           block->evout, &block->nevout, block->rpar, &block->nrpar,
3305                           block->ipar, &block->nipar,
3306                           (double **)block->inptr, block->insz, &block->nin,
3307                           (double **)block->outptr, block->outsz, &block->nout,
3308                           block->g, &block->ng);
3309             }
3310
3311             /* adjust values of output register */
3312             for (in = 0; in < block->nevout; ++in)
3313             {
3314                 block->evout[in] = block->evout[in] - *t;
3315             }
3316
3317             break;
3318         }
3319
3320         /***********************/
3321         /* function type 10004 */
3322         /***********************/
3323         case 10004 :
3324         {
3325             /* get pointer of the function type 4*/
3326             loc4 = (ScicosF4) loc;
3327
3328             (*loc4)(block, *flag);
3329
3330             break;
3331         }
3332
3333         /***********/
3334         /* default */
3335         /***********/
3336         default :
3337         {
3338             sciprint(_("Undefined Function type\n"));
3339             *flag = -1000;
3340             return; /* exit */
3341         }
3342     }
3343     // sciprint("callf end  flag=%d\n",*flag);
3344     /* Implicit Solver & explicit block & flag==0 */
3345     /* adjust continuous state vector after call */
3346     if (solver == 100 && block->type < 10000 && *flag == 0)
3347     {
3348         block->xd  = ptr_d;
3349         if (flagi != 7)
3350         {
3351             for (k = 0; k < block->nx; k++)
3352             {
3353                 block->res[k] = block->res[k] - block->xd[k];
3354             }
3355         }
3356         else
3357         {
3358             for (k = 0; k < block->nx; k++)
3359             {
3360                 block->xd[k] = block->res[k];
3361             }
3362         }
3363     }
3364
3365     /* debug block */
3366     if (cosd > 1)
3367     {
3368         if (debug_block > -1)
3369         {
3370             if (*flag < 0) return; /* error in block */
3371             if (cosd != 3) sciprint(_("Leaving block %d \n"), C2F(curblk).kfun);
3372             call_debug_scicos(block, flag, flagi, debug_block);
3373             /*call_debug_scicos(flag,kf,flagi,debug_block);*/
3374         }
3375     }
3376 } /* callf */
3377 /*--------------------------------------------------------------------------*/
3378 /* call_debug_scicos */
3379 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk)
3380 {
3381     voidf loc ;
3382     int solver = C2F(cmsolver).solver, k = 0;
3383     ScicosF4 loc4;
3384     double *ptr_d = NULL;
3385
3386     C2F(cosdebugcounter).counter = C2F(cosdebugcounter).counter + 1;
3387     C2F(scsptr).ptr = Blocks[deb_blk].scsptr;
3388
3389     loc  = Blocks[deb_blk].funpt; /* GLOBAL */
3390     loc4 = (ScicosF4) loc;
3391
3392     /* continuous state */
3393     if (solver == 100 && block->type < 10000 && *flag == 0)
3394     {
3395         ptr_d = block->xd;
3396         block->xd  = block->res;
3397     }
3398
3399     (*loc4)(block, *flag);
3400
3401     /* Implicit Solver & explicit block & flag==0 */
3402     /* adjust continuous state vector after call */
3403     if (solver == 100 && block->type < 10000 && *flag == 0)
3404     {
3405         block->xd  = ptr_d;
3406         if (flagi != 7)
3407         {
3408             for (k = 0; k < block->nx; k++)
3409             {
3410                 block->res[k] = block->res[k] - block->xd[k];
3411             }
3412         }
3413         else
3414         {
3415             for (k = 0; k < block->nx; k++)
3416             {
3417                 block->xd[k] = block->res[k];
3418             }
3419         }
3420     }
3421
3422     if (*flag < 0) sciprint(_("Error in the Debug block \n"));
3423 } /* call_debug_scicos */
3424 /*--------------------------------------------------------------------------*/
3425 /* simblk */
3426 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data)
3427 {
3428     double tx = 0., *x = NULL, *xd = NULL;
3429     int i = 0, nantest = 0;
3430
3431     tx = (double) t;
3432     x =  (double *) NV_DATA_S(yy);
3433     xd = (double *) NV_DATA_S(yp);
3434
3435     for (i = 0; i < *neq; i++)   xd[i] = 0; /* à la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
3436     C2F(ierode).iero = 0;
3437     *ierr = 0;
3438     odoit(&tx, x, xd, xd);
3439     C2F(ierode).iero = *ierr;
3440
3441     if (*ierr == 0)
3442     {
3443         nantest = 0;
3444         for (i = 0; i < *neq; i++) /* NaN checking */
3445         {
3446             if ((xd[i] - xd[i] != 0))
3447             {
3448                 sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
3449                 nantest = 1;
3450                 break;
3451             }
3452         }
3453         if (nantest == 1) return 349; /* recoverable error; */
3454     }
3455
3456     return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
3457
3458 } /* simblk */
3459 /*--------------------------------------------------------------------------*/
3460 /* grblk */
3461 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data)
3462 {
3463     double tx = 0., *x = NULL;
3464     int jj = 0, nantest = 0;
3465
3466     tx = (double) t;
3467     x =  (double *) NV_DATA_S(yy);
3468
3469     C2F(ierode).iero = 0;
3470     *ierr = 0;
3471
3472     zdoit(&tx, x, x, (double*) gout);
3473
3474     if (*ierr == 0)
3475     {
3476         nantest = 0;
3477         for (jj = 0; jj < ng; jj++)
3478             if (gout[jj] - gout[jj] != 0)
3479             {
3480                 sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
3481                 nantest = 1;
3482                 break;
3483             } /* NaN checking */
3484         if (nantest == 1) return 350; /* recoverable error; */
3485     }
3486     C2F(ierode).iero = *ierr;
3487
3488     return 0;
3489 } /* grblk */
3490 /*--------------------------------------------------------------------------*/
3491 /* simblklsodar */
3492 static int simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res)
3493 {
3494     double tx = 0.;
3495     int i = 0, nantest = 0;
3496
3497     tx = (double) *tOld;
3498
3499     for (i = 0; i < *nequations; ++i) res[i] = 0; /* à la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
3500     C2F(ierode).iero = 0;
3501     *ierr = 0;
3502     odoit(&tx, actual, res, res);
3503     C2F(ierode).iero = *ierr;
3504
3505     if (*ierr == 0)
3506     {
3507         nantest = 0;
3508         for (i = 0; i < *nequations; i++) /* NaN checking */
3509         {
3510             if ((res[i] - res[i] != 0))
3511             {
3512                 sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
3513                 nantest = 1;
3514                 break;
3515             }
3516         }
3517         if (nantest == 1) return 349; /* recoverable error; */
3518     }
3519
3520     return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
3521
3522 } /* simblklsodar */
3523 /*--------------------------------------------------------------------------*/
3524 /* grblklsodar */
3525 static int grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngc, realtype * res)
3526 {
3527     double tx = 0.;
3528     int jj = 0, nantest = 0;
3529
3530    tx = (double) *tOld;
3531
3532     C2F(ierode).iero = 0;
3533     *ierr = 0;
3534
3535     zdoit(&tx, actual, actual, res);
3536
3537     if (*ierr == 0)
3538     {
3539         nantest = 0;
3540         for (jj = 0; jj < *ngc; jj++)
3541             if (res[jj] - res[jj] != 0)
3542             {
3543                 sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
3544                 nantest = 1;
3545                 break;
3546             } /* NaN checking */
3547         if (nantest == 1) return 350; /* recoverable error; */
3548     }
3549     C2F(ierode).iero = *ierr;
3550
3551     return 0;
3552 } /* grblklsodar */
3553 /*--------------------------------------------------------------------------*/
3554 /* simblkdaskr */
3555 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata)
3556 {
3557     double tx = 0.;
3558     double *xc = NULL, *xcdot = NULL, *residual = NULL;
3559     realtype alpha = 0.;
3560
3561     UserData data;
3562
3563     realtype hh = 0.;
3564     int qlast = 0;
3565     int jj = 0, flag = 0, nantest = 0;
3566
3567     data = (UserData) rdata;
3568
3569     if (get_phase_simulation() == 1)
3570     {
3571         /* Just to update mode in a very special case, i.e., when initialization using modes fails.
3572         in this case, we relax all modes and try again one more time.
3573         */
3574         zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), (double *)data->gwork);
3575     }
3576
3577     hh = ZERO;
3578     flag = IDAGetCurrentStep(data->ida_mem, &hh);
3579     if (flag < 0)
3580     {
3581         *ierr = 200 + (-flag);
3582         return (*ierr);
3583     };
3584
3585     qlast = 0;
3586     flag = IDAGetCurrentOrder(data->ida_mem, &qlast);
3587     if (flag < 0)
3588     {
3589         *ierr = 200 + (-flag);
3590         return (*ierr);
3591     };
3592
3593     alpha = ZERO;
3594     for (jj = 0; jj < qlast; jj++)
3595         alpha = alpha - ONE / (jj + 1);
3596     if (hh != 0)
3597         // CJ=-alpha/hh;  // For function Get_Jacobian_cj
3598         CJJ = -alpha / hh;
3599     else
3600     {
3601         *ierr = 217;
3602         return (*ierr);
3603     }
3604     xc = (double *)  NV_DATA_S(yy);
3605     xcdot = (double *) NV_DATA_S(yp);
3606     residual = (double *) NV_DATA_S(resval);
3607     tx = (double) tres;
3608
3609     C2F(dcopy)(neq, xcdot, &c__1, residual, &c__1);
3610     *ierr = 0;
3611     C2F(ierode).iero = 0;
3612     odoit(&tx, xc, xcdot, residual);
3613
3614     C2F(ierode).iero = *ierr;
3615
3616     if (*ierr == 0)
3617     {
3618         nantest = 0;
3619         for (jj = 0; jj < *neq; jj++)
3620             if (residual[jj] - residual[jj] != 0) /* NaN checking */
3621             {
3622                 //sciprint(_("\nWarning: The residual function #%d returns a NaN"),jj);
3623                 nantest = 1;
3624                 break;
3625             }
3626         if (nantest == 1) return 257; /* recoverable error; */
3627     }
3628
3629     return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
3630 }/* simblkdaskr */
3631 /*--------------------------------------------------------------------------*/
3632 /* grblkdaskr */
3633 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data)
3634 {
3635     double tx = 0.;
3636     int jj = 0, nantest = 0;
3637
3638     tx = (double) t;
3639
3640     *ierr = 0;
3641     C2F(ierode).iero = 0;
3642     zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), (double *)gout);
3643     if (*ierr == 0)
3644     {
3645         nantest = 0; /* NaN checking */
3646         for (jj = 0; jj < ng; jj++)
3647         {
3648             if (gout[jj] - gout[jj] != 0)
3649             {
3650                 sciprint(_("\nWarning: The zero-crossing function #%d returns a NaN"), jj);
3651                 nantest = 1;
3652                 break;
3653             }
3654         }
3655         if (nantest == 1)
3656         {
3657             return 258; /* recoverable error; */
3658         }
3659     }
3660     C2F(ierode).iero = *ierr;
3661     return (*ierr);
3662 }/* grblkdaskr */
3663 /*--------------------------------------------------------------------------*/
3664 /* Subroutine addevs */
3665 static void addevs(double t, int *evtnb, int *ierr1)
3666 {
3667     static int i = 0, j = 0;
3668
3669     /* Function Body */
3670     *ierr1 = 0;
3671     if (evtspt[*evtnb] != -1)
3672     {
3673         if ((evtspt[*evtnb] == 0) && (*pointi == *evtnb))
3674         {
3675             tevts[*evtnb] = t;
3676             return;
3677         }
3678         else
3679         {
3680             if (*pointi == *evtnb)
3681             {
3682                 *pointi = evtspt[*evtnb]; /* remove from chain */
3683             }
3684             else
3685             {
3686                 i = *pointi;
3687                 while (*evtnb != evtspt[i])
3688                 {
3689                     i = evtspt[i];
3690                 }
3691                 evtspt[i] = evtspt[*evtnb]; /* remove old evtnb from chain */
3692                 if (TCritWarning == 0)
3693                 {
3694                     sciprint(_("\n Warning: an event is reprogrammed at t=%g by removing another"), t );
3695                     sciprint(_("\n         (already programmed) event. There may be an error in"));
3696                     sciprint(_("\n         your model. Please check your model\n"));
3697                     TCritWarning = 1;
3698                 }
3699                 do_cold_restart(); /* the erased event could be a critical
3700                                                                    event, so do_cold_restart is added to
3701                                                                    refresh the critical event table */
3702             }
3703             evtspt[*evtnb] = 0;
3704             tevts[*evtnb] = t;
3705         }
3706     }
3707     else
3708     {
3709         evtspt[*evtnb] = 0;
3710         tevts[*evtnb] = t;
3711     }
3712     if (*pointi == 0)
3713     {
3714         *pointi = *evtnb;
3715         return;
3716     }
3717     if (t < tevts[*pointi])
3718     {
3719         evtspt[*evtnb] = *pointi;
3720         *pointi = *evtnb;
3721         return;
3722     }
3723     i = *pointi;
3724
3725 L100:
3726     if (evtspt[i] == 0)
3727     {
3728         evtspt[i] = *evtnb;
3729         return;
3730     }
3731     if (t >= tevts[evtspt[i]])
3732     {
3733         j = evtspt[i];
3734         if (evtspt[j] == 0)
3735         {
3736             evtspt[j] = *evtnb;
3737             return;
3738         }
3739         i = j;
3740         goto L100;
3741     }
3742     else
3743     {
3744         evtspt[*evtnb] = evtspt[i];
3745         evtspt[i] = *evtnb;
3746     }
3747 } /* addevs */
3748 /*--------------------------------------------------------------------------*/
3749 /* Subroutine putevs */
3750 void putevs(double *t, int *evtnb, int *ierr1)
3751 {
3752     /* Function Body */
3753     *ierr1 = 0;
3754     if (evtspt[*evtnb] != -1)
3755     {
3756         *ierr1 = 1;
3757         return;
3758     }
3759     else
3760     {
3761         evtspt[*evtnb] = 0;
3762         tevts[*evtnb] = *t;
3763     }
3764     if (*pointi == 0)
3765     {
3766         *pointi = *evtnb;
3767         return;
3768     }
3769     evtspt[*evtnb] = *pointi;
3770     *pointi = *evtnb;
3771 } /* putevs */
3772 /*--------------------------------------------------------------------------*/
3773 /* Subroutine idoit */
3774 static void idoit(double *told)
3775 {
3776     /* initialisation (propagation of constant blocks outputs) */
3777     /*     Copyright INRIA */
3778
3779     int i2 = 0;
3780     scicos_flag flag = 0;
3781     int i = 0, j = 0;
3782     int ierr1 = 0;
3783
3784     ScicosImport *scs_imp = NULL;
3785     int *kf = NULL;
3786
3787     scs_imp = getscicosimportptr();
3788
3789     flag = 1;
3790     for (j = 0; j < * (scs_imp->niord); j++)
3791     {
3792         kf = &scs_imp->iord[j];
3793         C2F(curblk).kfun = *kf; /* */
3794         if (scs_imp->funtyp[*kf - 1] > -1)
3795         {
3796             /* continuous state */
3797             if (scs_imp->blocks[*kf - 1].nx > 0)
3798             {
3799                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
3800                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
3801             }
3802             scs_imp->blocks[*kf - 1].nevprt = scs_imp->iord[j + * (scs_imp->niord)];
3803             /*callf(told, xd, x, x,x,&flag);*/
3804             callf(told, &scs_imp->blocks[*kf - 1], &flag);
3805             if (flag < 0)
3806             {
3807                 *ierr = 5 - flag;
3808                 return;
3809             }
3810         }
3811         if (scs_imp->blocks[*kf - 1].nevout > 0)
3812         {
3813             if (scs_imp->funtyp[*kf - 1] < 0)
3814             {
3815                 i = synchro_nev(scs_imp, *kf, ierr);
3816                 if (*ierr != 0)
3817                 {
3818                     return;
3819                 }
3820                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
3821                 putevs(told, &i2, &ierr1);
3822                 if (ierr1 != 0)
3823                 {
3824                     /* event conflict */
3825                     *ierr = 3;
3826                     return;
3827                 }
3828                 doit(told);
3829                 if (*ierr != 0)
3830                 {
3831                     return;
3832                 }
3833             }
3834         }
3835     }
3836 } /* idoit_ */
3837 /*--------------------------------------------------------------------------*/
3838 /* Subroutine doit */
3839 static void doit(double *told)
3840 {
3841     /* propagation of blocks outputs on discrete activations */
3842     /*     Copyright INRIA */
3843
3844     int i = 0, i2 = 0;
3845     scicos_flag flag = 0;
3846     int nord = 0;
3847     int ierr1 = 0;
3848     int ii = 0, kever = 0;
3849
3850     ScicosImport *scs_imp = NULL;
3851     int *kf = NULL;
3852
3853     scs_imp = getscicosimportptr();
3854
3855     kever = *pointi;
3856     *pointi = evtspt[kever];
3857     evtspt[kever] = -1;
3858
3859     nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
3860     if (nord == 0)
3861     {
3862         return;
3863     }
3864
3865     for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1 ; ii++)
3866     {
3867         kf = &scs_imp->ordclk[ii - 1];
3868         C2F(curblk).kfun = *kf;
3869         if (scs_imp->funtyp[*kf - 1] > -1)
3870         {
3871             /* continuous state */
3872             if (scs_imp->blocks[*kf - 1].nx > 0)
3873             {
3874                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
3875                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
3876             }
3877             scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
3878             flag = 1;
3879             /*callf(told, xd, x, x,x,&flag);*/
3880             callf(told, &scs_imp->blocks[*kf - 1], &flag);
3881             if (flag < 0)
3882             {
3883                 *ierr = 5 - flag;
3884                 return;
3885             }
3886         }
3887
3888         /* Initialize tvec */
3889         if (scs_imp->blocks[*kf - 1].nevout > 0)
3890         {
3891             if (scs_imp->funtyp[*kf - 1] < 0)
3892             {
3893                 i = synchro_nev(scs_imp, *kf, ierr);
3894                 if (*ierr != 0)
3895                 {
3896                     return;
3897                 }
3898                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
3899                 putevs(told, &i2, &ierr1);
3900                 if (ierr1 != 0)
3901                 {
3902                     /* event conflict */
3903                     *ierr = 3;
3904                     return;
3905                 }
3906                 doit(told);
3907                 if (*ierr != 0)
3908                 {
3909                     return;
3910                 }
3911             }
3912         }
3913     }
3914 } /* doit_ */
3915 /*--------------------------------------------------------------------------*/
3916 /* Subroutine cdoit */
3917 static void cdoit(double *told)
3918 {
3919     /* propagation of continuous blocks outputs */
3920     /*     Copyright INRIA */
3921
3922     int i2 = 0;
3923     scicos_flag flag = 0;
3924     int ierr1 = 0;
3925     int i = 0, j = 0;
3926
3927     ScicosImport *scs_imp = NULL;
3928     int *kf = NULL;
3929
3930     scs_imp = getscicosimportptr();
3931
3932     /* Function Body */
3933     for (j = 0; j < * (scs_imp->ncord); j++)
3934     {
3935         kf = &scs_imp->cord[j];
3936         C2F(curblk).kfun = *kf;
3937         /* continuous state */
3938         if (scs_imp->blocks[*kf - 1].nx > 0)
3939         {
3940             scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
3941             scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
3942         }
3943         scs_imp->blocks[*kf - 1].nevprt = scs_imp->cord[j + * (scs_imp->ncord)];
3944         if (scs_imp->funtyp[*kf - 1] > -1)
3945         {
3946             flag = 1;
3947             /*callf(told, xd, x, x,x,&flag);*/
3948             callf(told, &scs_imp->blocks[*kf - 1], &flag);
3949             if (flag < 0)
3950             {
3951                 *ierr = 5 - flag;
3952                 return;
3953             }
3954         }
3955
3956         /* Initialize tvec */
3957         if (scs_imp->blocks[*kf - 1].nevout > 0)
3958         {
3959             if (scs_imp->funtyp[*kf - 1] < 0)
3960             {
3961                 i = synchro_nev(scs_imp, *kf, ierr);
3962                 if (*ierr != 0)
3963                 {
3964                     return;
3965                 }
3966                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
3967                 putevs(told, &i2, &ierr1);
3968                 if (ierr1 != 0)
3969                 {
3970                     /* event conflict */
3971                     *ierr = 3;
3972                     return;
3973                 }
3974                 doit(told);
3975                 if (*ierr != 0)
3976                 {
3977                     return;
3978                 }
3979             }
3980         }
3981     }
3982 } /* cdoit_ */
3983 /*--------------------------------------------------------------------------*/
3984 /* Subroutine ddoit */
3985 static void ddoit(double *told)
3986 {
3987     /* update states & event out on discrete activations */
3988     /*     Copyright INRIA */
3989
3990     int i2 = 0, j = 0;
3991     scicos_flag flag = 0;
3992     int kiwa = 0;
3993     int i = 0, i3 = 0, ierr1 = 0;
3994     int ii = 0, keve = 0;
3995
3996     ScicosImport *scs_imp = NULL;
3997     int *kf = NULL;
3998
3999     scs_imp = getscicosimportptr();
4000
4001     /* Function Body */
4002     kiwa = 0;
4003     edoit(told, &kiwa);
4004     if (*ierr != 0)
4005     {
4006         return;
4007     }
4008
4009     /* update continuous and discrete states on event */
4010     if (kiwa == 0)
4011     {
4012         return;
4013     }
4014     for (i = 0; i < kiwa; i++)
4015     {
4016         keve = iwa[i];
4017         if (critev[keve] != 0)
4018         {
4019             hot = 0;
4020         }
4021         i2 = scs_imp->ordptr[keve] - 1;
4022         for (ii = scs_imp->ordptr[keve - 1]; ii <= i2; ii++)
4023         {
4024             kf = &scs_imp->ordclk[ii - 1];
4025             C2F(curblk).kfun = *kf;
4026             /* continuous state */
4027             if (scs_imp->blocks[*kf - 1].nx > 0)
4028             {
4029                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4030                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4031             }
4032
4033             scs_imp->blocks[*kf - 1].nevprt = scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1];
4034
4035             if (scs_imp->blocks[*kf - 1].nevout > 0)
4036             {
4037                 if (scs_imp->funtyp[*kf - 1] >= 0)
4038                 {
4039                     /* initialize evout */
4040                     for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4041                     {
4042                         scs_imp->blocks[*kf - 1].evout[j] = -1;
4043                     }
4044                     flag = 3;
4045
4046                     if (scs_imp->blocks[*kf - 1].nevprt > 0) /* if event has continuous origin don't call*/
4047                     {
4048                         /*callf(told, xd, x, x ,x,&flag);*/
4049                         callf(told, &scs_imp->blocks[*kf - 1], &flag);
4050                         if (flag < 0)
4051                         {
4052                             *ierr = 5 - flag;
4053                             return;
4054                         }
4055                     }
4056
4057                     for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4058                     {
4059                         if (scs_imp->blocks[*kf - 1].evout[j] >= 0.)
4060                         {
4061                             i3 = j + scs_imp->clkptr[*kf - 1] ;
4062                             addevs(scs_imp->blocks[*kf - 1].evout[j] + (*told), &i3, &ierr1);
4063                             if (ierr1 != 0)
4064                             {
4065                                 /* event conflict */
4066                                 *ierr = 3;
4067                                 return;
4068                             }
4069                         }
4070                     }
4071                 }
4072             }
4073
4074             if (scs_imp->blocks[*kf - 1].nevprt > 0)
4075             {
4076                 if (scs_imp->blocks[*kf - 1].nx + scs_imp->blocks[*kf - 1].nz + scs_imp->blocks[*kf - 1].noz > 0 || \
4077                         *scs_imp->blocks[*kf - 1].work != NULL)
4078                 {
4079                     /*  if a hidden state exists, must also call (for new scope eg)  */
4080                     /*  to avoid calling non-real activations */
4081                     flag = 2;
4082                     /*callf(told, xd, x, x,x,&flag);*/
4083                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
4084                     if (flag < 0)
4085                     {
4086                         *ierr = 5 - flag;
4087                         return;
4088                     }
4089                 }
4090             }
4091             else
4092             {
4093                 if (*scs_imp->blocks[*kf - 1].work != NULL)
4094                 {
4095                     flag = 2;
4096                     scs_imp->blocks[*kf - 1].nevprt = 0; /* in case some hidden continuous blocks need updating */
4097                     /*callf(told, xd, x, x,x,&flag);*/
4098                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
4099                     if (flag < 0)
4100                     {
4101                         *ierr = 5 - flag;
4102                         return;
4103                     }
4104                 }
4105             }
4106         }
4107     }
4108 } /* ddoit_ */
4109 /*--------------------------------------------------------------------------*/
4110 /* Subroutine edoit */
4111 static void edoit(double *told, int *kiwa)
4112 {
4113     /* update blocks output on discrete activations */
4114     /*     Copyright INRIA */
4115
4116     int i2 = 0;
4117     scicos_flag flag = 0;
4118     int ierr1 = 0, i = 0;
4119     int kever = 0, ii = 0;
4120
4121     ScicosImport *scs_imp = NULL;
4122     int *kf = NULL;
4123     int nord = 0;
4124
4125     scs_imp = getscicosimportptr();
4126
4127     /* Function Body */
4128     kever = *pointi;
4129
4130     *pointi = evtspt[kever];
4131     evtspt[kever] = -1;
4132
4133     nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
4134     if (nord == 0)
4135     {
4136         return;
4137     }
4138     iwa[*kiwa] = kever;
4139     ++(*kiwa);
4140     for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
4141     {
4142         kf = &scs_imp->ordclk[ii - 1];
4143         C2F(curblk).kfun = *kf;
4144
4145         if (scs_imp->funtyp[*kf - 1] > -1)
4146         {
4147             /* continuous state */
4148             if (scs_imp->blocks[*kf - 1].nx > 0)
4149             {
4150                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4151                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4152             }
4153
4154             scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4155
4156             flag = 1;
4157             /*callf(told, xd, x, x,x,&flag);*/
4158             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4159             if (flag < 0)
4160             {
4161                 *ierr = 5 - flag;
4162                 return;
4163             }
4164         }
4165
4166         /* Initialize tvec */
4167         if (scs_imp->blocks[*kf - 1].nevout > 0)
4168         {
4169             if (scs_imp->funtyp[*kf - 1] < 0)
4170             {
4171                 i = synchro_nev(scs_imp, *kf, ierr);
4172                 if (*ierr != 0)
4173                 {
4174                     return;
4175                 }
4176                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4177                 putevs(told, &i2, &ierr1);
4178                 if (ierr1 != 0)
4179                 {
4180                     /* event conflict */
4181                     *ierr = 3;
4182                     return;
4183                 }
4184                 edoit(told, kiwa);
4185                 if (*ierr != 0)
4186                 {
4187                     return;
4188                 }
4189             }
4190         }
4191     }
4192 } /* edoit_ */
4193 /*--------------------------------------------------------------------------*/
4194 /* Subroutine odoit */
4195 static void odoit(double *told, double *xt, double *xtd, double *residual)
4196 {
4197     /* update blocks derivative of continuous time block */
4198     /*     Copyright INRIA */
4199
4200     int i2 = 0;
4201     scicos_flag flag = 0;
4202     int keve = 0, kiwa = 0;
4203     int ierr1 = 0, i = 0;
4204     int ii = 0, jj = 0;
4205
4206     ScicosImport *scs_imp = NULL;
4207     int *kf = NULL;
4208
4209     scs_imp = getscicosimportptr();
4210
4211     /* Function Body */
4212     kiwa = 0;
4213     for (jj = 0; jj < * (scs_imp->noord); jj++)
4214     {
4215         kf = &scs_imp->oord[jj];
4216         C2F(curblk).kfun = *kf;
4217         /* continuous state */
4218         if (scs_imp->blocks[*kf - 1].nx > 0)
4219         {
4220             scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4221             scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4222             scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
4223         }
4224
4225         scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
4226         if (scs_imp->funtyp[*kf - 1] > -1)
4227         {
4228             flag = 1;
4229             /*callf(told, xtd, xt, residual,x,&flag);*/
4230             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4231             if (flag < 0)
4232             {
4233                 *ierr = 5 - flag;
4234                 return;
4235             }
4236         }
4237
4238         if (scs_imp->blocks[*kf - 1].nevout > 0)
4239         {
4240             if (scs_imp->funtyp[*kf - 1] < 0)
4241             {
4242                 if (scs_imp->blocks[*kf - 1].nmode > 0)
4243                 {
4244                     i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
4245                 }
4246                 else
4247                 {
4248                     i = synchro_nev(scs_imp, *kf, ierr);
4249                     if (*ierr != 0)
4250                     {
4251                         return;
4252                     }
4253                     i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4254                 }
4255                 putevs(told, &i2, &ierr1);
4256                 if (ierr1 != 0)
4257                 {
4258                     /* event conflict */
4259                     *ierr = 3;
4260                     return;
4261                 }
4262                 ozdoit(told, xt, xtd, &kiwa);
4263                 if (*ierr != 0)
4264                 {
4265                     return;
4266                 }
4267             }
4268         }
4269     }
4270
4271     /*  update states derivatives */
4272     for (ii = 0; ii < * (scs_imp->noord); ii++)
4273     {
4274         kf = &scs_imp->oord[ii];
4275         C2F(curblk).kfun = *kf;
4276         if (scs_imp->blocks[*kf - 1].nx > 0 || \
4277                 *scs_imp->blocks[*kf - 1].work != NULL)
4278         {
4279             /* work tests if a hidden state exists, used for delay block */
4280             flag = 0;
4281             /* continuous state */
4282             if (scs_imp->blocks[*kf - 1].nx > 0)
4283             {
4284                 scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4285                 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4286                 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
4287             }
4288             scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
4289             /*callf(told, xtd, xt, residual,xt,&flag);*/
4290             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4291
4292             if (flag < 0)
4293             {
4294                 *ierr = 5 - flag;
4295                 return;
4296             }
4297         }
4298     }
4299
4300     for (i = 0; i < kiwa; i++)
4301     {
4302         keve = iwa[i];
4303         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
4304         {
4305             kf = &scs_imp->ordclk[ii - 1];
4306             C2F(curblk).kfun = *kf;
4307             if (scs_imp->blocks[*kf - 1].nx > 0 || \
4308                     *scs_imp->blocks[*kf - 1].work != NULL)
4309             {
4310                 /* work tests if a hidden state exists */
4311                 flag = 0;
4312                 /* continuous state */
4313                 if (scs_imp->blocks[*kf - 1].nx > 0)
4314                 {
4315                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4316                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4317                     scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
4318                 }
4319                 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4320                 /*callf(told, xtd, xt, residual,xt,&flag);*/
4321                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4322
4323                 if (flag < 0)
4324                 {
4325                     *ierr = 5 - flag;
4326                     return;
4327                 }
4328             }
4329         }
4330     }
4331 } /* odoit_ */
4332 /*--------------------------------------------------------------------------*/
4333 /* Subroutine reinitdoit */
4334 static void reinitdoit(double *told)
4335 {
4336     /* update blocks xproperties of continuous time block */
4337     /*     Copyright INRIA */
4338
4339     int i2 = 0;
4340     scicos_flag flag = 0;
4341     int keve = 0, kiwa = 0;
4342     int ierr1 = 0, i = 0;
4343     int ii = 0, jj = 0;
4344
4345     ScicosImport *scs_imp = NULL;
4346     int *kf = NULL;
4347
4348     scs_imp = getscicosimportptr();
4349
4350     /* Function Body */
4351     kiwa = 0;
4352     for (jj = 0; jj < * (scs_imp->noord); jj++)
4353     {
4354         kf = &scs_imp->oord[jj];
4355         C2F(curblk).kfun = *kf;
4356         /* continuous state */
4357         if (scs_imp->blocks[*kf - 1].nx > 0)
4358         {
4359             scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4360             scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4361         }
4362         scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
4363         if (scs_imp->funtyp[*kf - 1] > -1)
4364         {
4365             flag = 1;
4366             /*callf(told, xd, x, x,x,&flag);*/
4367             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4368             if (flag < 0)
4369             {
4370                 *ierr = 5 - flag;
4371                 return;
4372             }
4373         }
4374
4375         if (scs_imp->blocks[*kf - 1].nevout > 0 && scs_imp->funtyp[*kf - 1] < 0)
4376         {
4377             i = synchro_nev(scs_imp, *kf, ierr);
4378             if (*ierr != 0)
4379             {
4380                 return;
4381             }
4382             if (scs_imp->blocks[*kf - 1].nmode > 0)
4383             {
4384                 scs_imp->blocks[*kf - 1].mode[0] = i;
4385             }
4386             i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4387             putevs(told, &i2, &ierr1);
4388             if (ierr1 != 0)
4389             {
4390                 /* event conflict */
4391                 *ierr = 3;
4392                 return;
4393             }
4394             doit(told);
4395             if (*ierr != 0)
4396             {
4397                 return;
4398             }
4399         }
4400     }
4401
4402     /* re-initialize */
4403     for (ii = 0; ii < * (scs_imp->noord); ii++)
4404     {
4405         kf = &scs_imp->oord[ii];
4406         C2F(curblk).kfun = *kf;
4407         if (scs_imp->blocks[*kf - 1].nx > 0)
4408         {
4409             flag = 7;
4410             scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4411             scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4412             scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4413             scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
4414             scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
4415             /*callf(told, xd, x, xd,x,&flag);*/
4416             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4417             if (flag < 0)
4418             {
4419                 *ierr = 5 - flag;
4420                 return;
4421             }
4422         }
4423     }
4424
4425     for (i = 0; i < kiwa; i++)
4426     {
4427         keve = iwa[i];
4428         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
4429         {
4430             kf = &scs_imp->ordclk[ii - 1];
4431             C2F(curblk).kfun = *kf;
4432             if (scs_imp->blocks[*kf - 1].nx > 0)
4433             {
4434                 flag = 7;
4435                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4436                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4437                 scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4438                 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4439                 scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
4440                 /*callf(told, xd, x, xd,x,&flag);*/
4441                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4442                 if (flag < 0)
4443                 {
4444                     *ierr = 5 - flag;
4445                     return;
4446                 }
4447             }
4448         }
4449     }
4450 } /* reinitdoit_ */
4451 /*--------------------------------------------------------------------------*/
4452 /* Subroutine ozdoit */
4453 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa)
4454 {
4455     /* update blocks output of continuous time block on discrete activations */
4456     /*     Copyright INRIA */
4457
4458     int i2 = 0;
4459     scicos_flag flag = 0;
4460     int nord = 0;
4461     int ierr1 = 0, i = 0;
4462     int ii = 0, kever = 0;
4463
4464     ScicosImport *scs_imp = NULL;
4465     int *kf = NULL;
4466
4467     scs_imp = getscicosimportptr();
4468
4469     /* Function Body */
4470     kever = *pointi;
4471     *pointi = evtspt[kever];
4472     evtspt[kever] = -1;
4473
4474     nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
4475     if (nord == 0)
4476     {
4477         return;
4478     }
4479     iwa[*kiwa] = kever;
4480     ++(*kiwa);
4481
4482     for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
4483     {
4484         kf = &scs_imp->ordclk[ii - 1];
4485         C2F(curblk).kfun = *kf;
4486         if (scs_imp->funtyp[*kf - 1] > -1)
4487         {
4488             /* continuous state */
4489             if (scs_imp->blocks[*kf - 1].nx > 0)
4490             {
4491                 scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4492                 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4493             }
4494             scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4495             flag = 1;
4496             /*callf(told, xtd, xt, xt,x,&flag);*/
4497             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4498             if (flag < 0)
4499             {
4500                 *ierr = 5 - flag;
4501                 return;
4502             }
4503         }
4504
4505         /* Initialize tvec */
4506         if (scs_imp->blocks[*kf - 1].nevout > 0)
4507         {
4508             if (scs_imp->funtyp[*kf - 1] < 0)
4509             {
4510                 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
4511                 {
4512                     i = synchro_nev(scs_imp, *kf, ierr);
4513                     if (*ierr != 0)
4514                     {
4515                         return;
4516                     }
4517                 }
4518                 else
4519                 {
4520                     i = scs_imp->blocks[*kf - 1].mode[0];
4521                 }
4522                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4523                 putevs(told, &i2, &ierr1);
4524                 if (ierr1 != 0)
4525                 {
4526                     /* event conflict */
4527                     *ierr = 3;
4528                     return;
4529                 }
4530                 ozdoit(told, xt, xtd, kiwa);
4531             }
4532         }
4533     }
4534 } /* ozdoit_ */
4535 /*--------------------------------------------------------------------------*/
4536 /* Subroutine zdoit */
4537 static void zdoit(double *told, double *xt, double *xtd, double *g)
4538 {
4539     /* update blocks zcross of continuous time block  */
4540     /*     Copyright INRIA */
4541     int i2 = 0;
4542     scicos_flag flag = 0;
4543     int keve = 0, kiwa = 0;
4544     int ierr1 = 0, i = 0, j = 0;
4545     int ii = 0, jj = 0;
4546
4547     ScicosImport *scs_imp = NULL;
4548     int *kf = NULL;
4549
4550     scs_imp = getscicosimportptr();
4551
4552     /* Function Body */
4553     for (i = 0; i < * (scs_imp->ng); i++)
4554     {
4555         g[i] = 0.;
4556     }
4557
4558     kiwa = 0;
4559     for (jj = 0; jj < * (scs_imp->nzord); jj++)
4560     {
4561         kf = &scs_imp->zord[jj];
4562         C2F(curblk).kfun = *kf;
4563         /* continuous state */
4564         if (scs_imp->blocks[*kf - 1].nx > 0)
4565         {
4566             scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4567             scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4568         }
4569         scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[jj + * (scs_imp->nzord)];
4570
4571         if (scs_imp->funtyp[*kf - 1] > -1)
4572         {
4573             flag = 1;
4574             /*callf(told, xtd, xt, xt,xt,&flag);*/
4575             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4576             if (flag < 0)
4577             {
4578                 *ierr = 5 - flag;
4579                 return;
4580             }
4581         }
4582
4583         /* Initialize tvec */
4584         if (scs_imp->blocks[*kf - 1].nevout > 0)
4585         {
4586             if (scs_imp->funtyp[*kf - 1] < 0)
4587             {
4588                 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
4589                 {
4590                     i = synchro_nev(scs_imp, *kf, ierr);
4591                     if (*ierr != 0)
4592                     {
4593                         return;
4594                     }
4595                 }
4596                 else
4597                 {
4598                     i = scs_imp->blocks[*kf - 1].mode[0];
4599                 }
4600                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4601                 putevs(told, &i2, &ierr1);
4602                 if (ierr1 != 0)
4603                 {
4604                     /* event conflict */
4605                     *ierr = 3;
4606                     return;
4607                 }
4608                 ozdoit(told, xt, xtd, &kiwa);
4609                 if (*ierr != 0)
4610                 {
4611                     return;
4612                 }
4613             }
4614         }
4615     }
4616
4617     /* update zero crossing surfaces */
4618     for (ii = 0; ii < * (scs_imp->nzord); ii++)
4619     {
4620         kf = &scs_imp->zord[ii];
4621         C2F(curblk).kfun = *kf;
4622         if (scs_imp->blocks[*kf - 1].ng > 0)
4623         {
4624             /* update g array ptr */
4625             scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
4626             if (scs_imp->funtyp[*kf - 1] > 0)
4627             {
4628                 flag = 9;
4629                 /* continuous state */
4630                 if (scs_imp->blocks[*kf - 1].nx > 0)
4631                 {
4632                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4633                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4634                 }
4635                 scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[ii + * (scs_imp->nzord)];
4636                 /*callf(told, xtd, xt, xtd,g,&flag);*/
4637                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4638                 if (flag < 0)
4639                 {
4640                     *ierr = 5 - flag;
4641                     return;
4642                 }
4643             }
4644             else
4645             {
4646                 j = synchro_g_nev(scs_imp, g, *kf, ierr);
4647                 if (*ierr != 0)
4648                 {
4649                     return;
4650                 }
4651                 if ( (phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0) )
4652                 {
4653                     scs_imp->blocks[*kf - 1].mode[0] = j;
4654                 }
4655             }
4656
4657             // scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
4658
4659         }
4660     }
4661
4662     for (i = 0; i < kiwa; i++)
4663     {
4664         keve = iwa[i];
4665         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
4666         {
4667             kf = &scs_imp->ordclk[ii - 1];
4668             C2F(curblk).kfun = *kf;
4669             if (scs_imp->blocks[*kf - 1].ng > 0)
4670             {
4671                 /* update g array ptr */
4672                 scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
4673                 if (scs_imp->funtyp[*kf - 1] > 0)
4674                 {
4675                     flag = 9;
4676                     /* continuous state */
4677                     if (scs_imp->blocks[*kf - 1].nx > 0)
4678                     {
4679                         scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4680                         scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4681                     }
4682                     scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4683                     /*callf(told, xtd, xt, xtd,g,&flag);*/
4684                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
4685                     if (flag < 0)
4686                     {
4687                         *ierr = 5 - flag;
4688                         return;
4689                     }
4690                 }
4691                 else
4692                 {
4693                     j = synchro_g_nev(scs_imp, g, *kf, ierr);
4694                     if (*ierr != 0)
4695                     {
4696                         return;
4697                     }
4698                     if ((phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0))
4699                     {
4700                         scs_imp->blocks[*kf - 1].mode[0] = j;
4701                     }
4702                 }
4703
4704                 //scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
4705             }
4706         }
4707     }
4708 } /* zdoit_ */
4709 /*--------------------------------------------------------------------------*/
4710 /* Subroutine Jdoit */
4711 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job)
4712 {
4713     /* update blocks jacobian of continuous time block  */
4714     /*     Copyright INRIA */
4715
4716     int i2 = 0;
4717     scicos_flag flag = 0;
4718     int keve = 0, kiwa = 0;
4719     int ierr1 = 0, i = 0;
4720     int ii = 0, jj = 0;
4721
4722     ScicosImport *scs_imp = NULL;
4723     int *kf = NULL;
4724
4725     scs_imp = getscicosimportptr();
4726
4727     /* Function Body */
4728     kiwa = 0;
4729     for (jj = 0; jj < * (scs_imp->noord); jj++)
4730     {
4731         kf = &scs_imp->oord[jj];
4732         C2F(curblk).kfun = *kf;
4733         scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
4734         if (scs_imp->funtyp[*kf - 1] > -1)
4735         {
4736             flag = 1;
4737             /* applying desired output */
4738             if ((*job == 2) && (scs_imp->oord[jj] == AJacobian_block))
4739             {
4740             }
4741             else
4742                 /* continuous state */
4743                 if (scs_imp->blocks[*kf - 1].nx > 0)
4744                 {
4745                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4746                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4747                     scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
4748                 }
4749
4750             /*callf(told, xtd, xt, residual,x,&flag);*/
4751             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4752             if (flag < 0)
4753             {
4754                 *ierr = 5 - flag;
4755                 return;
4756             }
4757         }
4758
4759         if (scs_imp->blocks[*kf - 1].nevout > 0)
4760         {
4761             if (scs_imp->funtyp[*kf - 1] < 0)
4762             {
4763                 if (scs_imp->blocks[*kf - 1].nmode > 0)
4764                 {
4765                     i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
4766                 }
4767                 else
4768                 {
4769                     i = synchro_nev(scs_imp, *kf, ierr);
4770                     if (*ierr != 0)
4771                     {
4772                         return;
4773                     }
4774                     i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4775                 }
4776                 putevs(told, &i2, &ierr1);
4777                 if (ierr1 != 0)
4778                 {
4779                     /* event conflict */
4780                     *ierr = 3;
4781                     return;
4782                 }
4783                 ozdoit(told, xt, xtd, &kiwa);
4784             }
4785         }
4786     }
4787
4788     /* update states derivatives */
4789     for (ii = 0; ii < * (scs_imp->noord); ii++)
4790     {
4791         kf = &scs_imp->oord[ii];
4792         C2F(curblk).kfun = *kf;
4793         if (scs_imp->blocks[*kf - 1].nx > 0 || \
4794                 *scs_imp->blocks[*kf - 1].work != NULL)
4795         {
4796             /* work tests if a hidden state exists, used for delay block */
4797             flag = 0;
4798             if (((*job == 1) && (scs_imp->oord[ii] == AJacobian_block)) || (*job != 1))
4799             {
4800                 if (*job == 1)  flag = 10;
4801                 /* continuous state */
4802                 if (scs_imp->blocks[*kf - 1].nx > 0)
4803                 {
4804                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4805                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4806                     scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
4807                 }
4808                 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
4809                 /*callf(told, xtd, xt, residual,xt,&flag);*/
4810                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4811             }
4812             if (flag < 0)
4813             {
4814                 *ierr = 5 - flag;
4815                 return;
4816             }
4817         }
4818     }
4819
4820     for (i = 0; i < kiwa; i++)
4821     {
4822         keve = iwa[i];
4823         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1;  ii++)
4824         {
4825             kf = &scs_imp->ordclk[ii - 1];
4826             C2F(curblk).kfun = *kf;
4827             if (scs_imp->blocks[*kf - 1].nx > 0 || \
4828                     *scs_imp->blocks[*kf - 1].work != NULL)
4829             {
4830                 /* work tests if a hidden state exists */
4831                 flag = 0;
4832                 if (((*job == 1) && (scs_imp->oord[ii - 1] == AJacobian_block)) || (*job != 1))
4833                 {
4834                     if (*job == 1)  flag = 10;
4835                     /* continuous state */
4836                     if (scs_imp->blocks[*kf - 1].nx > 0)
4837                     {
4838                         scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
4839                         scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4840                         scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
4841                     }
4842                     scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4843                     /*callf(told, xtd, xt, residual,xt,&flag);*/
4844                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
4845                 }
4846                 if (flag < 0)
4847                 {
4848                     *ierr = 5 - flag;
4849                     return;
4850                 }
4851             }
4852         }
4853     }
4854 } /* Jdoit_ */
4855 /*--------------------------------------------------------------------------*/
4856 /* Subroutine synchro_nev */
4857 static int synchro_nev(ScicosImport *scs_imp, int kf, int *ierr)
4858 {
4859     /* synchro blocks computation  */
4860     /*     Copyright INRIA */
4861     SCSREAL_COP *outtbdptr = NULL;     /*to store double of outtb*/
4862     SCSINT8_COP *outtbcptr = NULL;     /*to store int8 of outtb*/
4863     SCSINT16_COP *outtbsptr = NULL;    /*to store int16 of outtb*/
4864     SCSINT32_COP *outtblptr = NULL;    /*to store int32 of outtb*/
4865     SCSUINT8_COP *outtbucptr = NULL;   /*to store unsigned int8 of outtb */
4866     SCSUINT16_COP *outtbusptr = NULL;  /*to store unsigned int16 of outtb */
4867     SCSUINT32_COP *outtbulptr = NULL;  /*to store unsigned int32 of outtb */
4868
4869     int cond = 0;
4870     int i = 0; /* return 0 by default */
4871
4872     /* variable for param */
4873     int *outtbtyp = 0;
4874     void **outtbptr = NULL;
4875     int *funtyp = 0;
4876     int *inplnk = 0;
4877     int *inpptr = 0;
4878
4879     /* get param ptr */
4880     outtbtyp = scs_imp->outtbtyp;
4881     outtbptr = scs_imp->outtbptr;
4882     funtyp = scs_imp->funtyp;
4883     inplnk = scs_imp->inplnk;
4884     inpptr = scs_imp->inpptr;
4885
4886     /* if-then-else blk */
4887     if (funtyp[kf - 1] == -1)
4888     {
4889         switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
4890         {
4891             case SCSREAL_N    :
4892                 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4893                 cond = (*outtbdptr <= 0.);
4894                 break;
4895
4896             case SCSCOMPLEX_N :
4897                 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4898                 cond = (*outtbdptr <= 0.);
4899                 break;
4900
4901             case SCSINT8_N    :
4902                 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4903                 cond = (*outtbcptr <= 0);
4904                 break;
4905
4906             case SCSINT16_N   :
4907                 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4908                 cond = (*outtbsptr <= 0);
4909                 break;
4910
4911             case SCSINT32_N   :
4912                 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4913                 cond = (*outtblptr <= 0);
4914                 break;
4915
4916             case SCSUINT8_N   :
4917                 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4918                 cond = (*outtbucptr <= 0);
4919                 break;
4920
4921             case SCSUINT16_N  :
4922                 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4923                 cond = (*outtbusptr <= 0);
4924                 break;
4925
4926             case SCSUINT32_N  :
4927                 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4928                 cond = (*outtbulptr <= 0);
4929                 break;
4930
4931             default  : /* Add a message here */
4932                 *ierr = 25;
4933                 return 0;
4934         }
4935         if (cond)
4936         {
4937             i = 2;
4938         }
4939         else
4940         {
4941             i = 1;
4942         }
4943     }
4944     /* eselect blk */
4945     else if (funtyp[kf - 1] == -2)
4946     {
4947         switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
4948         {
4949             case SCSREAL_N    :
4950                 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4951                 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
4952                 break;
4953
4954             case SCSCOMPLEX_N :
4955                 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4956                 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
4957                 break;
4958
4959             case SCSINT8_N    :
4960                 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4961                 i = Max(Min((int) * outtbcptr,
4962                             scs_imp->blocks[kf - 1].nevout), 1);
4963                 break;
4964
4965             case SCSINT16_N   :
4966                 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4967                 i = Max(Min((int) * outtbsptr, scs_imp->blocks[kf - 1].nevout), 1);
4968                 break;
4969
4970             case SCSINT32_N   :
4971                 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4972                 i = Max(Min((int) * outtblptr, scs_imp->blocks[kf - 1].nevout), 1);
4973                 break;
4974
4975             case SCSUINT8_N   :
4976                 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4977                 i = Max(Min((int) * outtbucptr, scs_imp->blocks[kf - 1].nevout), 1);
4978                 break;
4979
4980             case SCSUINT16_N  :
4981                 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4982                 i = Max(Min((int) * outtbusptr, scs_imp->blocks[kf - 1].nevout), 1);
4983                 break;
4984
4985             case SCSUINT32_N  :
4986                 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
4987                 i = Max(Min((int) * outtbulptr, scs_imp->blocks[kf - 1].nevout), 1);
4988                 break;
4989
4990             default  : /* Add a message here */
4991                 *ierr = 25;
4992                 return 0;
4993         }
4994     }
4995     return i;
4996 } /* synchro_nev */
4997 /*--------------------------------------------------------------------------*/
4998 /* Subroutine synchro_g_nev */
4999 static int synchro_g_nev(ScicosImport *scs_imp, double *g, int kf, int *ierr)
5000 {
5001     /* synchro blocks with zcross computation  */
5002     /*     Copyright INRIA */
5003     SCSREAL_COP *outtbdptr = NULL;     /*to store double of outtb*/
5004     SCSINT8_COP *outtbcptr = NULL;     /*to store int8 of outtb*/
5005     SCSINT16_COP *outtbsptr = NULL;    /*to store int16 of outtb*/
5006     SCSINT32_COP *outtblptr = NULL;    /*to store int32 of outtb*/
5007     SCSUINT8_COP *outtbucptr = NULL;   /*to store unsigned int8 of outtb */
5008     SCSUINT16_COP *outtbusptr = NULL;  /*to store unsigned int16 of outtb */
5009     SCSUINT32_COP *outtbulptr = NULL;  /*to store unsigned int32 of outtb */
5010
5011     int cond = 0;
5012     int i = 0; /* return 0 by default */
5013     int jj = 0;
5014
5015     /* variable for param */
5016     int *outtbtyp = NULL;
5017     void **outtbptr = NULL;
5018     int *funtyp = NULL;
5019     int *inplnk = NULL;
5020     int *inpptr = NULL;
5021     int *zcptr = NULL;
5022
5023     /* get param ptr */
5024     outtbtyp = scs_imp->outtbtyp;
5025     outtbptr = scs_imp->outtbptr;
5026     funtyp = scs_imp->funtyp;
5027     inplnk = scs_imp->inplnk;
5028     inpptr = scs_imp->inpptr;
5029     zcptr = scs_imp->zcptr;
5030
5031     /* if-then-else blk */
5032     if (funtyp[kf - 1] == -1)
5033     {
5034         switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5035         {
5036             case SCSREAL_N    :
5037                 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5038                 g[zcptr[kf - 1] - 1] = *outtbdptr;
5039                 cond = (*outtbdptr <= 0.);
5040                 break;
5041
5042             case SCSCOMPLEX_N :
5043                 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5044                 g[zcptr[kf - 1] - 1] = *outtbdptr;
5045                 cond = (*outtbdptr <= 0.);
5046                 break;
5047
5048             case SCSINT8_N    :
5049                 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5050                 g[zcptr[kf - 1] - 1] = (double) * outtbcptr;
5051                 cond = (*outtbcptr <= 0);
5052                 break;
5053
5054             case SCSINT16_N   :
5055                 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5056                 g[zcptr[kf - 1] - 1] = (double) * outtbsptr;
5057                 cond = (*outtbsptr <= 0);
5058                 break;
5059
5060             case SCSINT32_N   :
5061                 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5062                 g[zcptr[kf - 1] - 1] = (double) * outtblptr;
5063                 cond = (*outtblptr <= 0);
5064                 break;
5065
5066             case SCSUINT8_N   :
5067                 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5068                 g[zcptr[kf - 1] - 1] = (double) * outtbucptr;
5069                 cond = (*outtbucptr <= 0);
5070                 break;
5071
5072             case SCSUINT16_N  :
5073                 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5074                 g[zcptr[kf - 1] - 1] = (double) * outtbusptr;
5075                 cond = (*outtbusptr <= 0);
5076                 break;
5077
5078             case SCSUINT32_N  :
5079                 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5080                 g[zcptr[kf - 1] - 1] = (double) * outtbulptr;
5081                 cond = (*outtbulptr <= 0);
5082                 break;
5083
5084             default  : /* Add a message here */
5085                 *ierr = 25;
5086                 return 0;
5087         }
5088         if (cond)
5089         {
5090             i = 2;
5091         }
5092         else
5093         {
5094             i = 1;
5095         }
5096     }
5097     /* eselect blk */
5098     else if (funtyp[kf - 1] == -2)
5099     {
5100         switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5101         {
5102             case SCSREAL_N    :
5103                 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5104                 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5105                 {
5106                     g[zcptr[kf - 1] - 1 + jj] = *outtbdptr - (double)(jj + 2);
5107                 }
5108                 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5109                 break;
5110
5111             case SCSCOMPLEX_N :
5112                 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5113                 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5114                 {
5115                     g[zcptr[kf - 1] - 1 + jj] = *outtbdptr - (double)(jj + 2);
5116                 }
5117                 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5118                 break;
5119
5120             case SCSINT8_N    :
5121                 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5122                 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5123                 {
5124                     g[zcptr[kf - 1] - 1 + jj] = (double) * outtbcptr - (double)(jj + 2);
5125                 }
5126                 i = Max(Min((int) * outtbcptr, scs_imp->blocks[kf - 1].nevout), 1);
5127                 break;
5128
5129             case SCSINT16_N   :
5130                 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5131                 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5132                 {
5133                     g[zcptr[kf - 1] - 1 + jj] = (double) * outtbsptr - (double)(jj + 2);
5134                 }
5135                 i = Max(Min((int) * outtbsptr, scs_imp->blocks[kf - 1].nevout), 1);
5136                 break;
5137
5138             case SCSINT32_N   :
5139                 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5140                 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5141                 {
5142                     g[zcptr[kf - 1] - 1 + jj] = (double) * outtblptr - (double)(jj + 2);
5143                 }
5144                 i = Max(Min((int) * outtblptr, scs_imp->blocks[kf - 1].nevout), 1);
5145                 break;
5146
5147             case SCSUINT8_N   :
5148                 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5149                 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5150                 {
5151                     g[zcptr[kf - 1] - 1 + jj] = (double) * outtbucptr - (double)(jj + 2);
5152                 }
5153                 i = Max(Min((int) * outtbucptr, scs_imp->blocks[kf - 1].nevout), 1);
5154                 break;
5155
5156             case SCSUINT16_N  :
5157                 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5158                 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5159                 {
5160                     g[zcptr[kf - 1] - 1 + jj] = (double) * outtbusptr - (double)(jj + 2);
5161                 }
5162                 i = Max(Min((int) * outtbusptr, scs_imp->blocks[kf - 1].nevout), 1);
5163                 break;
5164
5165             case SCSUINT32_N  :
5166                 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5167                 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5168                 {
5169                     g[zcptr[kf - 1] - 1 + jj] = (double) * outtbulptr - (double)(jj + 2);
5170                 }
5171                 i = Max(Min((int) * outtbulptr, scs_imp->blocks[kf - 1].nevout), 1);
5172                 break;
5173
5174             default  : /* Add a message here */
5175                 *ierr = 25;
5176                 return 0;
5177         }
5178     }
5179     return i;
5180 } /* synchro_g_nev */
5181 /*--------------------------------------------------------------------------*/
5182 /* FREE_blocks */
5183 static void FREE_blocks()
5184 {
5185     int kf = 0;
5186     for (kf = 0; kf < nblk; ++kf)
5187     {
5188         if (Blocks[kf].insz != NULL)
5189         {
5190             FREE(Blocks[kf].insz);
5191         }
5192         else
5193         {
5194             break;
5195         }
5196         if (Blocks[kf].inptr != NULL)
5197         {
5198             FREE(Blocks[kf].inptr);
5199         }
5200         else
5201         {
5202             break;
5203         }
5204         if (Blocks[kf].outsz != NULL)
5205         {
5206             FREE(Blocks[kf].outsz);
5207         }
5208         else
5209         {
5210             break;
5211         }
5212         if (Blocks[kf].outptr != NULL)
5213         {
5214             FREE(Blocks[kf].outptr);
5215         }
5216         else
5217         {
5218             break;
5219         }
5220         if (Blocks[kf].oparsz != NULL)
5221         {
5222             FREE(Blocks[kf].oparsz);
5223         }
5224         else
5225         {
5226             break;
5227         }
5228         if (Blocks[kf].ozsz != NULL)
5229         {
5230             FREE(Blocks[kf].ozsz);
5231         }
5232         else
5233         {
5234             break;
5235         }
5236         if (Blocks[kf].label != NULL)
5237         {
5238             FREE(Blocks[kf].label);
5239         }
5240         else
5241         {
5242             break;
5243         }
5244         if (Blocks[kf].evout != NULL)
5245         {
5246             FREE(Blocks[kf].evout);
5247         }
5248         else
5249         {
5250             break;
5251         }
5252     }
5253     FREE(Blocks);
5254
5255     if (nx > 0) FREE(xprop);
5256
5257     if (nmod > 0) FREE(mod);
5258
5259     if (ng > 0) FREE(g);
5260
5261     return;
5262 } /* FREE_blocks */
5263 /*--------------------------------------------------------------------------*/
5264 /* Subroutine funnum */
5265 int C2F(funnum)(char * fname)
5266 {
5267     int i = 0, ln = 0;
5268     int loc = -1;
5269     while ( tabsim[i].name != (char *) NULL)
5270     {
5271         if ( strcmp(fname, tabsim[i].name) == 0 ) return(i + 1);
5272         i++;
5273     }
5274     ln = (int)strlen(fname);
5275     C2F(iislink)(fname, &loc);
5276     C2F(iislink)(fname, &loc);
5277     if (loc >= 0) return(ntabsim + (int)loc + 1);
5278     return(0);
5279 }/* funnum */
5280 /*--------------------------------------------------------------------------*/
5281 int get_phase_simulation(void)
5282 {
5283     return phase;
5284 }
5285 /*--------------------------------------------------------------------------*/
5286 void do_cold_restart(void)
5287 {
5288     hot = 0;
5289     return;
5290 }
5291 /*--------------------------------------------------------------------------*/
5292 /* get_scicos_time : return the current
5293 * simulation time
5294 */
5295 double get_scicos_time(void)
5296 {
5297     return scicos_time;
5298 }
5299 /*--------------------------------------------------------------------------*/
5300 /* get_block_number : return the current
5301 * block number
5302 */
5303 int get_block_number(void)
5304 {
5305     return C2F(curblk).kfun;
5306 }
5307 /*--------------------------------------------------------------------------*/
5308 /* set_block_error : set an error number
5309 * for block_error
5310 */
5311 void set_block_error(int err)
5312 {
5313     *block_error = err;
5314 }
5315 /*--------------------------------------------------------------------------*/
5316 /* Coserror : copy an error message
5317 * in coserr.buf an set block_error to
5318 * -16
5319 */
5320 #if _MSC_VER
5321 #ifndef vsnprintf
5322 #define vsnprintf _vsnprintf
5323 #endif
5324 #endif
5325
5326
5327 void Coserror(const char *fmt, ...)
5328 {
5329     int retval;
5330     va_list ap;
5331
5332     va_start(ap, fmt);
5333
5334 #ifdef vsnprintf
5335     retval = vsnprintf(coserr.buf, 4095, fmt, ap);
5336 #else
5337     retval = vsprintf(coserr.buf, fmt, ap);
5338 #endif
5339
5340     if (retval == -1)
5341     {
5342         coserr.buf[0] = '\0';
5343     }
5344
5345     va_end(ap);
5346
5347     /* coserror use error number 10 */
5348     *block_error = -5;
5349 }
5350 /*--------------------------------------------------------------------------*/
5351 /* get_block_error : get the block error
5352 * number
5353 */
5354 int get_block_error()
5355 {
5356     return *block_error;
5357 }
5358 /*--------------------------------------------------------------------------*/
5359 void end_scicos_sim()
5360 {
5361     C2F(coshlt).halt = 2;
5362     return;
5363 }
5364 /*--------------------------------------------------------------------------*/
5365 /* get_pointer_xproperty */
5366 int* get_pointer_xproperty()
5367 {
5368     return &xprop[-1 + xptr[C2F(curblk).kfun]];
5369 }
5370 /*--------------------------------------------------------------------------*/
5371 /* get_Npointer_xproperty */
5372 int get_npointer_xproperty()
5373 {
5374     return Blocks[C2F(curblk).kfun - 1].nx;
5375 }
5376 /*--------------------------------------------------------------------------*/
5377 /* set_pointer_xproperty */
5378 void set_pointer_xproperty(int* pointer)
5379 {
5380     int i;
5381     for (i = 0; i < Blocks[C2F(curblk).kfun - 1].nx; i++)
5382     {
5383         Blocks[C2F(curblk).kfun - 1].xprop[i] = pointer[i];
5384     }
5385 }
5386 /*--------------------------------------------------------------------------*/
5387 /* Jacobian */
5388 void Set_Jacobian_flag(int flag)
5389 {
5390     Jacobian_Flag = flag;
5391     return;
5392 }
5393 /*--------------------------------------------------------------------------*/
5394 /* Get_Jacobian_ci et Get_Jacobian_cj were called by the C file only produced
5395 by Modelicac v 1.11.2 */
5396 /* double Get_Jacobian_ci(void)
5397 {
5398 return CI;
5399 } */
5400 /*--------------------------------------------------------------------------*/
5401 /* double Get_Jacobian_cj(void)
5402 {
5403 return CJ;
5404 } */
5405 /*--------------------------------------------------------------------------*/
5406 /* Fonction called by the C file produced by Modelicac 1.7.3 and 1.12.1 */
5407 double Get_Jacobian_parameter(void)
5408 {
5409     return CJJ;
5410 }
5411 /*--------------------------------------------------------------------------*/
5412 double Get_Scicos_SQUR(void)
5413 {
5414     return  SQuround;
5415 }
5416 /*--------------------------------------------------------------------------*/
5417 static int Jacobians(long int Neq, realtype tt, N_Vector yy, N_Vector yp,
5418                      N_Vector resvec, realtype cj, void *jdata, DenseMat Jacque,
5419                      N_Vector tempv1, N_Vector tempv2, N_Vector tempv3)
5420 {
5421     double  ttx = 0;
5422     double *xc = NULL, *xcdot = NULL, *residual = NULL;
5423     /*  char chr;*/
5424     int i = 0, j = 0, n = 0, nx = 0, ni = 0, no = 0, nb = 0, m = 0, flag = 0;
5425     double *RX = NULL, *Fx = NULL, *Fu = NULL, *Gx = NULL, *Gu = NULL, *ERR1 = NULL, *ERR2 = NULL;
5426     double *Hx = NULL, *Hu = NULL, *Kx = NULL, *Ku = NULL, *HuGx = NULL, *FuKx = NULL, *FuKuGx = NULL, *HuGuKx = NULL;
5427     double ysave = 0;
5428     int job = 0;
5429     double **y = NULL;
5430     double **u = NULL;
5431     /*  taill1= 3*n+(n+ni)*(n+no)+nx(2*nx+ni+2*m+no)+m*(2*m+no+ni)+2*ni*no*/
5432     double inc = 0., inc_inv = 0., xi = 0., xpi = 0., srur = 0.;
5433     realtype *Jacque_col = NULL;
5434
5435     UserData data;
5436     realtype hh = 0.;
5437     N_Vector ewt;
5438     double *ewt_data = NULL;
5439
5440     *ierr = 0;
5441
5442     data = (UserData) jdata;
5443     ewt = data->ewt;
5444
5445     flag = IDAGetCurrentStep(data->ida_mem, &hh);
5446     if (flag < 0)
5447     {
5448         *ierr = 200 + (-flag);
5449         return (*ierr);
5450     };
5451
5452     flag = IDAGetErrWeights(data->ida_mem, ewt);
5453     if (flag < 0)
5454     {
5455         *ierr = 200 + (-flag);
5456         return (*ierr);
5457     };
5458
5459     ewt_data = NV_DATA_S(ewt);
5460     xc   = (double *) N_VGetArrayPointer(yy);
5461     xcdot  = (double *) N_VGetArrayPointer(yp);
5462     /*residual=(double *) NV_DATA_S(resvec);*/
5463     ttx = (double)tt;
5464     // CJ=(double)cj;  // for fonction Get_Jacobian_cj
5465     CJJ = (double)cj;  // returned by Get_Jacobian_parameter
5466
5467     srur = (double) RSqrt(UNIT_ROUNDOFF);
5468
5469     if (AJacobian_block > 0)
5470     {
5471         nx = Blocks[AJacobian_block - 1].nx; /* quant on est là cela signifie que AJacobian_block>0 */
5472         no = Blocks[AJacobian_block - 1].nout;
5473         ni = Blocks[AJacobian_block - 1].nin;
5474         y = (double **)Blocks[AJacobian_block - 1].outptr; /*for compatibility */
5475         u = (double **)Blocks[AJacobian_block - 1].inptr; /*warning pointer of y and u have changed to void ***/
5476     }
5477     else
5478     {
5479         nx = 0;
5480         no = 0;
5481         ni = 0;
5482     }
5483     n = Neq;
5484     nb = nblk;
5485     m = n - nx;
5486
5487     residual = (double *)data->rwork;
5488     ERR1 = residual + n;
5489     ERR2 = ERR1 + n;
5490     RX = ERR2 + n;
5491     Fx = RX + (n + ni) * (n + no); /* car (nx+ni)*(nx+no) peut etre > `a n*n*/
5492     Fu = Fx + nx * nx;
5493     Gx = Fu + nx * ni;
5494     Gu = Gx + no * nx;
5495     Hx = Gu + no * ni;
5496     Hu = Hx + m * m;
5497     Kx = Hu + m * no;
5498     Ku = Kx + ni * m;
5499     HuGx = Ku + ni * no;
5500     FuKx = HuGx + m * nx;
5501     FuKuGx = FuKx + nx * m;
5502     HuGuKx = FuKuGx + nx * nx;
5503     /* HuGuKx+m*m; =>  m*m=size of HuGuKx */
5504     /* ------------------ Numerical Jacobian--->> Hx,Kx */
5505
5506     /* read residuals;*/
5507     job = 0;
5508     Jdoit(&ttx, xc, xcdot, residual, &job);
5509     if (*ierr < 0) return -1;
5510
5511     /* "residual" already contains the current residual,
5512     so the first call to Jdoit can be remoevd*/
5513
5514     for (i = 0; i < m; i++)
5515         for (j = 0; j < ni; j++)
5516             Kx[j + i * ni] = u[j][0];
5517
5518     for (i = 0; i < m; i++)
5519     {
5520         xi = xc[i];
5521         xpi = xcdot[i];
5522         inc = MAX( srur * MAX( ABS(xi), ABS(hh * xpi)), ONE / ewt_data[i] );
5523         if (hh * xpi < ZERO) inc = -inc;
5524         inc = (xi + inc) - xi;
5525
5526         /* if (CI==0) {
5527         inc = MAX( srur * ABS(hh*xpi),ONE );
5528         if (hh*xpi < ZERO) inc = -inc;
5529         inc = (xpi + inc) - xi;
5530         } */
5531         // xc[i] += CI*inc;
5532         // xcdot[i] += CJ*inc;
5533         xc[i] += inc;
5534         xcdot[i] += CJJ * inc;
5535         /*a= Max(abs(H[0]*xcdot[i]),abs(1.0/Ewt[i]));
5536         b= Max(1.0,abs(xc[i]));
5537         del=SQUR[0]*Max(a,b);    */
5538         job = 0; /* read residuals */
5539         Jdoit(&ttx, xc, xcdot, ERR2, &job);
5540         if (*ierr < 0) return -1;
5541         inc_inv = ONE / inc;
5542         for (j = 0; j < m; j++)
5543             Hx[m * i + j] = (ERR2[j] - residual[j]) * inc_inv;
5544         for (j = 0; j < ni; j++)
5545             Kx[j + i * ni] = (u[j][0] - Kx[j + i * ni]) * inc_inv;
5546         xc[i] = xi;
5547         xcdot[i] = xpi;
5548     }
5549     /*----- Numerical Jacobian--->> Hu,Ku */
5550
5551     if (AJacobian_block == 0)
5552     {
5553         for (j = 0; j < m; j++)
5554         {
5555             Jacque_col = DENSE_COL(Jacque, j);
5556             for (i = 0; i < m; i++)
5557             {
5558                 Jacque_col[i] = Hx[i + j * m];
5559             }
5560         }
5561         C2F(ierode).iero = *ierr;
5562         return 0;
5563     }
5564     /****------------------***/
5565     job = 0;
5566     Jdoit(&ttx, xc, xcdot, ERR1, &job);
5567     for (i = 0; i < no; i++)
5568         for (j = 0; j < ni; j++)
5569             Ku[j + i * ni] = u[j][0];
5570
5571     for (i = 0; i < no; i++)
5572     {
5573         ysave = y[i][0];
5574         inc = srur * MAX( ABS(ysave), 1);
5575         inc = (ysave + inc) - ysave;
5576         /*del=SQUR[0]* Max(1.0,abs(y[i][0]));
5577         del=(y[i][0]+del)-y[i][0];*/
5578         y[i][0] += inc;
5579         job = 2; /* applying y[i][0] to the output of imp block*/
5580         Jdoit(&ttx, xc, xcdot, ERR2, &job);
5581         if (*ierr < 0) return -1;
5582         inc_inv = ONE / inc;
5583         for (j = 0; j < m; j++)
5584             Hu[m * i + j] = (ERR2[j] - ERR1[j]) * inc_inv;
5585         for (j = 0; j < ni; j++)
5586             Ku[j + i * ni] = (u[j][0] - Ku[j + i * ni]) * inc_inv;
5587         y[i][0] = ysave;
5588     }
5589     /*----------------------------------------------*/
5590     job = 1; /* read jacobian through flag=10; */
5591     Jdoit(&ttx, xc, xcdot, &Fx[-m], &job);/* Filling up the FX:Fu:Gx:Gu*/
5592     if (*block_error != 0) sciprint(_("\n error in Jacobian"));
5593     /*-------------------------------------------------*/
5594
5595     Multp(Fu, Ku, RX, nx, ni, ni, no);
5596     Multp(RX, Gx, FuKuGx, nx, no, no, nx);
5597
5598     for (j = 0; j < nx; j++)
5599     {
5600         Jacque_col = DENSE_COL(Jacque, j + m);
5601         for (i = 0; i < nx; i++)
5602         {
5603             Jacque_col[i + m] = Fx[i + j * nx] + FuKuGx[i + j * nx];
5604         }
5605     }
5606
5607     Multp(Hu, Gx, HuGx, m, no, no, nx);
5608
5609     for (i = 0; i < nx; i++)
5610     {
5611         Jacque_col = DENSE_COL(Jacque, i + m);
5612         for (j = 0; j < m; j++)
5613         {
5614             Jacque_col[j] = HuGx[j + i * m];
5615         }
5616     }
5617
5618     Multp(Fu, Kx, FuKx, nx, ni, ni, m);
5619
5620     for (i = 0; i < m; i++)
5621     {
5622         Jacque_col = DENSE_COL(Jacque, i);
5623         for (j = 0; j < nx; j++)
5624         {
5625             Jacque_col[j + m] = FuKx[j + i * nx];
5626         }
5627     }
5628
5629
5630     Multp(Hu, Gu, RX, m, no, no, ni);
5631     Multp(RX, Kx, HuGuKx, m, ni, ni, m);
5632
5633     for (j = 0; j < m; j++)
5634     {
5635         Jacque_col = DENSE_COL(Jacque, j);
5636         for (i = 0; i < m; i++)
5637         {
5638             Jacque_col[i] = Hx[i + j * m] + HuGuKx[i + j * m];
5639         }
5640     }
5641
5642     /*  chr='Z';   printf("\n t=%g",ttx); DISP(Z,n,n,chr);*/
5643     C2F(ierode).iero = *ierr;
5644     return 0;
5645
5646 }
5647 /*----------------------------------------------------*/
5648 static void Multp(double *A, double *B, double *R, int ra, int rb, int ca, int cb)
5649 {
5650     int i = 0, j = 0, k = 0;
5651     /*if (ca!=rb) sciprint(_("\n Error in matrix multiplication"));*/
5652     for (i = 0; i < ra; i++)
5653         for (j = 0; j < cb; j++)
5654         {
5655             R[i + ra * j] = 0.0;
5656             for (k = 0; k < ca; k++)
5657                 R[i + ra * j] += A[i + k * ra] * B[k + j * rb];
5658         }
5659     return;
5660 }
5661 /*--------------------------------------------------------------------------*/
5662 int read_xml_initial_states(int nvar, const char * xmlfile, char **ids, double *svars)
5663 {
5664     ezxml_t model, elements;
5665     int result = 0, i = 0;
5666     double vr = 0.;
5667
5668     if (nvar == 0) return 0;
5669     result = 0;
5670     for (i = 0; i < nvar; i++)
5671     {
5672         if (strcmp(ids[i], "") != 0)
5673         {
5674             result = 1;
5675             break;
5676         }
5677     }
5678     if (result == 0) return 0;
5679
5680     model = ezxml_parse_file(xmlfile);
5681
5682     if (model == NULL)
5683     {
5684         sciprint(_("Error: Cannot find file '%s'.\n"), xmlfile);
5685         return -1;/* file does not exist*/
5686     }
5687
5688     elements = ezxml_child(model, "elements");
5689     for (i = 0; i < nvar; i++)
5690     {
5691         vr = 0.0;
5692         result = read_id(&elements, ids[i], &vr);
5693         if (result == 1) svars[i] = vr;
5694     }
5695     ezxml_free(model);
5696     return 0;
5697 }
5698 /*--------------------------------------------------------------------------*/
5699 static int read_id(ezxml_t *elements, char *id, double *value)
5700 {
5701     char V1[100], V2[100];
5702     int ok = 0, i = 0, ln = 0;
5703
5704     if (strcmp(id, "") == 0) return 0;
5705     ok = search_in_child(elements, id, V1);
5706     if (ok == 0 )
5707     {
5708         /*sciprint(_("Cannot find: %s=%s  \n"),id,V1);      */
5709         return 0;
5710     }
5711     else
5712     {
5713         if (Convert_number(V1, value) != 0)
5714         {
5715             ln = (int)(strlen(V1));
5716             if (ln > 2)
5717             {
5718                 for (i = 1; i <= ln - 2; i++) V2[i - 1] = V1[i];
5719                 V2[ln - 2] = '\0';
5720                 ok = read_id(elements, V2, value);
5721                 return ok;
5722             }
5723             else return 0;
5724         }
5725         else
5726         {
5727             /*      printf("\n ---->>>%s= %g",V1,*value);*/
5728             return 1;
5729         }
5730     }
5731 }
5732 /*--------------------------------------------------------------------------*/
5733 int Convert_number(char *s, double *out)
5734 {
5735     char *endp = NULL;
5736     double d = 0.;
5737     long int l = 0;
5738     d = strtod(s, &endp);
5739     if (s != endp && *endp == '\0')
5740     {
5741         /*    printf("  It's a float with value %g ", d); */
5742         *out = d;
5743         return 0;
5744     }
5745     else
5746     {
5747         l = strtol(s, &endp, 0);
5748         if (s != endp && *endp == '\0')
5749         {
5750             /*printf("  It's an int with value %ld ", 1); */
5751             *out = (double)l;
5752             return 0;
5753         }
5754         else
5755         {
5756             /*printf("  string "); */
5757             return -1;
5758         }
5759     }
5760 }
5761 /*--------------------------------------------------------------------------*/
5762 int write_xml_states(int nvar, const char * xmlfile, char **ids, double *x)
5763 {
5764     ezxml_t model, elements;
5765     int result = 0, i = 0, err = 0;
5766     FILE *fd = NULL;
5767     char *s = NULL;
5768     char **xv = NULL;
5769
5770     if (nvar == 0) return 0;
5771     result = 0;
5772     for (i = 0; i < nvar; i++)
5773     {
5774         if (strcmp(ids[i], "") != 0)
5775         {
5776             result = 1;
5777             break;
5778         }
5779     }
5780     if (result == 0) return 0;
5781
5782     xv = MALLOC(nvar * sizeof(char*));
5783     for (i = 0; i < nvar; i++)
5784     {
5785         xv[i] = MALLOC(nvar * 100 * sizeof(char));
5786         sprintf(xv[i], "%g", x[i]);
5787     }
5788
5789     model = ezxml_parse_file(xmlfile);
5790     if (model == NULL)
5791     {
5792         sciprint(_("Error: Cannot find file '%s'.\n"), xmlfile);
5793         for (i = 0; i < nvar; i++)
5794         {
5795             FREE(xv[i]);
5796         }
5797         FREE(xv);
5798         return -1;/* file does not exist */
5799     }
5800
5801     elements = ezxml_child(model, "elements");
5802
5803     for (i = 0; i < nvar; i++)
5804     {
5805         if (strcmp(ids[i], "") == 0) continue;
5806         result = write_in_child(&elements, ids[i], xv[i]);
5807         if (result == 0 )
5808         {
5809             /* sciprint(_("cannot find %s in '%s' \n"),ids[i],xmlfile);      */
5810             /* err= -1;*/ /* Variable does not exist*/
5811         }
5812     }
5813
5814     s = ezxml_toxml(model);
5815     ezxml_free(model);
5816
5817
5818     wcfopen(fd, (char*)xmlfile, "wb");
5819     if (fd < 0)
5820     {
5821         sciprint(_("Error: cannot write to  '%s'  \n"), xmlfile);
5822         return -3;/* cannot write to file*/
5823     }
5824
5825     fputs (s, fd);
5826     fclose(fd);
5827
5828     return err;
5829 }
5830 /*--------------------------------------------------------------------------*/
5831 int C2F(fx)(double *x, double *residual) /* used for homotopy*/
5832 {
5833     double  *xdot = NULL, t = 0;
5834     xdot = x + *neq;
5835     t = 0;
5836     *ierr = 0;
5837     C2F(ierode).iero = 0;
5838     odoit(&t, x, xdot, residual);
5839     C2F(ierode).iero = *ierr;
5840     return (*ierr);
5841 }
5842 /*--------------------------------------------------------------------------*/
5843 int rho_(double *a, double *L, double *x, double *rho, double *rpar, int *ipar) /* used for homotopy*/
5844 {
5845     int i = 0, N = 0;
5846     N = *neq;
5847
5848     fx_(x, rho);
5849     for (i = 0; i < N; i++)
5850         rho[i] += (-1 + *L) * a[i];
5851     return 0;
5852 }
5853 /*--------------------------------------------------------------------------*/
5854 int rhojac_(double *a, double *lambda, double  *x, double  *jac, int *col, double *rpar, int *ipar) /* used for homotopy*/
5855 {
5856     /* MATRIX [d_RHO/d_LAMBDA, d_RHO/d_X_col] */
5857     int j = 0, N = 0;
5858     double *work = NULL;
5859     int job = 0;
5860     double inc = 0., inc_inv = 0., xi = 0., srur = 0.;
5861     N = *neq;
5862     if (*col == 1)
5863     {
5864         for (j = 0; j < N; j++)
5865             jac[j] = a[j];
5866     }
5867     else
5868     {
5869         if ((work = (double *) MALLOC(N * sizeof(double))) == NULL)
5870         {
5871             *ierr = 10000;
5872             return *ierr;
5873         }
5874         rho_(a, lambda, x, work, rpar, ipar);
5875         srur = 1e-10;
5876         xi = x[*col - 2];
5877         inc = srur * Max(fabs(xi), 1);
5878         inc = (xi + inc) - xi;
5879         x[*col - 2] += inc;
5880
5881         job = 0;
5882         rho_(a, lambda, x, jac, rpar, ipar);
5883         inc_inv = 1.0 / inc;
5884
5885         for (j = 0; j < N; j++)
5886             jac[j] = (jac[j] - work[j]) * inc_inv;
5887
5888         x[*col - 2] = xi;
5889         FREE(work);
5890     }
5891     return 0;
5892 }
5893 /*--------------------------------------------------------------------------*/
5894 int C2F(hfjac)(double *x, double *jac, int *col)
5895 {
5896     int N = 0, j = 0;
5897     double *work = NULL;
5898     double  *xdot = NULL;
5899     int job = 0;
5900     double inc = 0., inc_inv = 0., xi = 0., srur = 0.;
5901
5902     N = *neq;
5903     if ((work = (double *) MALLOC(N * sizeof(double))) == NULL)
5904     {
5905         *ierr = 10000;
5906         return&