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