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