Clean old static structs from dynamic_link
[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., 675 Mass Ave, Cambridge, MA 02139, 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 #include "scicos_math.h"
73 #include "sciblk2.h"
74 #include "sciblk4.h"
75 #include "dynlib_scicos.h"
76
77 #if defined(linux) && defined(__i386__)
78 #include "setPrecisionFPU.h"
79 #endif
80
81 #include "localization.h"
82 #include "charEncoding.h"
83 /*--------------------------------------------------------------------------*/
84 typedef struct {
85         void *ida_mem;
86         N_Vector ewt;
87         double *rwork;
88         int *iwork;
89         double *gwork; /* just added for a very special use: a
90                                    space passing to grblkdakr for zero crossing surfaces
91                                    when updating mode variables during initialization */
92 } *UserData;
93
94 SCICOS_IMPEXP SCSPTR_struct C2F(scsptr);
95 /*--------------------------------------------------------------------------*/
96
97 #define freeall                                 \
98         if (*neq>0) CVodeFree(&cvode_mem);              \
99         if (*neq>0) N_VDestroy_Serial(y);               \
100         if ( ng>0 ) FREE(jroot);                        \
101         if ( ng>0 ) FREE(zcros);
102
103
104 /* TJacque allocates by sundials */
105 #define freeallx                                \
106         if (*neq>0) free(TJacque);      \
107         if (*neq>0) FREE(data->rwork);          \
108         if (( ng>0 )&& (*neq>0)) FREE(data->gwork);     \
109         if (*neq>0) N_VDestroy_Serial(data->ewt);       \
110         if (*neq>0) FREE(data);                 \
111         if (*neq>0) IDAFree(&ida_mem);          \
112         if (*neq>0) N_VDestroy_Serial(IDx);             \
113         if (*neq>0) N_VDestroy_Serial(yp);              \
114         if (*neq>0) N_VDestroy_Serial(yy);              \
115         if ( ng>0 ) FREE(jroot);                        \
116         if ( ng>0 ) FREE(zcros);                        \
117         if (nmod>0) FREE(Mode_save);
118
119 #define freeouttbptr                            \
120         FREE(outtbd);                                   \
121         FREE(outtbc);                                   \
122         FREE(outtbs);                                   \
123         FREE(outtbl);                                   \
124         FREE(outtbuc);                          \
125         FREE(outtbus);                          \
126         FREE(outtbul);
127
128 #define freekinsol                              \
129         FREE(Mode_save);                                \
130         N_VDestroy_Serial(y);                           \
131         N_VDestroy_Serial(fscale);                      \
132         N_VDestroy_Serial(yscale);                      \
133         KINFree(&kin_mem); 
134
135
136 #define ONE   RCONST(1.0)
137 #define ZERO  RCONST(0.0)
138 #define T0    RCONST(0.0)   
139 /*--------------------------------------------------------------------------*/
140 /* Table of constant values */
141 static int c__90 = 90;
142 static int c__91 = 91;
143 static int c__0 = 0;
144 static double c_b14 = 0.;
145 static int c__1 = 1;
146
147 int TCritWarning = 0;
148
149 static double *t0 = NULL, *tf = NULL;
150 static double *x = NULL,*tevts = NULL,*xd = NULL,*g = NULL;
151 static  double Atol = 0., rtol = 0., ttol = 0., deltat = 0.,hmax = 0.;
152 static int *ierr = NULL;
153 static int *pointi = NULL;
154 static int *xptr = NULL,*modptr = NULL, *evtspt = NULL;
155 static int *funtyp = NULL, *inpptr = NULL, *outptr = NULL, *inplnk = NULL, *outlnk = NULL;
156 static int *clkptr = NULL, *ordptr = NULL, *ordclk = NULL, *cord = NULL, *iord = NULL, *oord = NULL,  *zord = NULL,  *critev = NULL,  *zcptr = NULL;
157 static int nblk = 0, nordptr = 0, nlnk = 0, nx = 0, ng = 0, ncord = 0, noord = 0, nzord = 0;
158 static int nordclk = 0,niord = 0,nmod = 0;
159 static int nelem = 0;
160 static int *mod = NULL;
161 static int *neq = NULL;
162 static int *xprop = NULL; /* xproperty */
163 static int debug_block = 0;
164 static int *iwa = NULL;
165 static int hot;
166
167 /* declaration of ptr for typed port */
168 static void **outtbptr = NULL; /*pointer array of object of outtb*/
169 static int *outtbsz = NULL;    /*size of object of outtb*/
170 static int *outtbtyp = NULL;   /*type of object of outtb*/
171
172 SCSREAL_COP *outtbdptr = NULL;     /*to store double of outtb*/
173 SCSINT8_COP *outtbcptr = NULL;     /*to store int8 of outtb*/
174 SCSINT16_COP *outtbsptr = NULL;    /*to store int16 of outtb*/
175 SCSINT32_COP *outtblptr = NULL;    /*to store int32 of outtb*/
176 SCSUINT8_COP *outtbucptr = NULL;   /*to store unsigned int8 of outtb */
177 SCSUINT16_COP *outtbusptr = NULL;  /*to store unsigned int16 of outtb */
178 SCSUINT32_COP *outtbulptr = NULL;  /*to store unsigned int32 of outtb */
179
180 static outtb_el *outtb_elem = NULL;
181
182 static scicos_block *Blocks = NULL;
183
184 /* pass to external variable for code generation */
185 /* reserved variable name */
186 int *block_error = NULL;
187 double scicos_time = 0.;
188 int phase = 0;
189 int Jacobian_Flag = 0;
190 // double CI = 0., CJ = 0.;  // doubles returned by Get_Jacobian_ci and Get_Jacobian_cj respectively 
191 double  CJJ = 0.;            // returned by Get_Jacobian_parameter
192 double SQuround = 0.;
193 /* Jacobian*/
194 static int AJacobian_block = 0;
195
196
197 /* Variable declaration moved to scicos.c because it was in the scicos-def.h therefore
198 * multiple declaration of the variable and linkers were complaining about duplicate
199 * symbols
200 */
201 SCICOS_IMPEXP COSDEBUGCOUNTER_struct C2F(cosdebugcounter);
202 SCICOS_IMPEXP RTFACTOR_struct C2F(rtfactor);
203 SCICOS_IMPEXP SOLVER_struct C2F(cmsolver);
204 SCICOS_IMPEXP CURBLK_struct C2F(curblk);
205 SCICOS_IMPEXP COSDEBUG_struct C2F(cosdebug);
206 SCICOS_IMPEXP COSHLT_struct C2F(coshlt);
207 SCICOS_IMPEXP DBCOS_struct C2F(dbcos);
208 SCICOS_IMPEXP COSTOL_struct C2F(costol);
209 SCICOS_IMPEXP COSERR_struct coserr;
210 /*--------------------------------------------------------------------------*/
211 static void FREE_blocks();
212 static void cosini(double *told);
213 static void cosend(double *told);
214 static void cossim(double *told);
215 static void cossimdaskr(double *told);
216 static void doit(double *told);
217 static void idoit(double *told);
218 static void zdoit(double *told, double *xt, double *xtd, double *g);
219 static void cdoit(double *told);
220 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa);
221 static void odoit(double *told, double *xt, double *xtd, double *residual);
222 static void ddoit(double *told);
223 static void edoit(double *told, int *kiwa);
224 static void reinitdoit(double *told);
225 static int CallKinsol(double *told);
226 static int simblk(realtype t,N_Vector yy,N_Vector yp, void *f_data);
227 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);
228 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data);
229 static void addevs(double t, int *evtnb, int *ierr1);
230 static int synchro_g_nev(ScicosImport *scs_imp,double *g,int kf,int *ierr);
231 static void Multp(double *A,double *B,double *R,int ra,int rb,int ca,int cb);
232 static int read_id(ezxml_t *elements,char *id,double *value);
233 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata);
234 static int Jacobians(long int Neq, realtype tt, N_Vector yy, N_Vector yp,
235                                          N_Vector resvec, realtype cj, void *jdata, DenseMat Jacque,
236                                          N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
237 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk);
238 static int synchro_nev(ScicosImport *scs_imp,int kf,int *ierr);
239 /*--------------------------------------------------------------------------*/
240 extern int C2F(dset)(int *n, double *dx, double *dy, int *incy);
241 extern int C2F(dcopy)(int *,double *,int *,double *,int *);
242 /*--------------------------------------------------------------------------*/
243 void putevs(double *t, int *evtnb, int *ierr1);
244 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job);
245 int simblkKinsol(N_Vector yy, N_Vector resval, void *rdata);
246 /*--------------------------------------------------------------------------*/
247 int C2F(scicos)(double *x_in, int *xptr_in, double *z__,
248                                 void **work,int *zptr,int *modptr_in,
249                                 void **oz,int *ozsz,int *oztyp,int *ozptr,
250                                 int *iz,int *izptr,double *t0_in,
251                                 double *tf_in,double *tevts_in,int *evtspt_in,
252                                 int *nevts,int *pointi_in,void **outtbptr_in,
253                                 int *outtbsz_in,int *outtbtyp_in,
254                                 outtb_el *outtb_elem_in,int *nelem1,int *nlnk1,
255                                 int *funptr,int *funtyp_in,int *inpptr_in,
256                                 int *outptr_in, int *inplnk_in,int *outlnk_in,
257                                 double *rpar,int *rpptr,int *ipar,int *ipptr,
258                                 void **opar,int *oparsz,int *opartyp,int *opptr,
259                                 int *clkptr_in,int *ordptr_in,int *nordptr1,
260                                 int *ordclk_in,int *cord_in,int *ncord1,
261                                 int *iord_in,int *niord1,int *oord_in,
262                                 int *noord1,int *zord_in,int *nzord1,
263                                 int *critev_in,int *nblk1,int *ztyp,
264                                 int *zcptr_in,int *subscr,int *nsubs,
265                                 double *simpar,int *flag__,int *ierr_out)
266 {
267         int i1,kf,lprt,in,out,job=1;
268
269         extern int C2F(msgs)();
270         extern int C2F(getscsmax)();
271
272         extern int C2F(makescicosimport)();
273         extern int C2F(clearscicosimport)();
274
275         static int mxtb = 0, ierr0 = 0, kfun0 = 0, i = 0, j = 0, k = 0, jj = 0;
276         static int ni = 0, no = 0;
277         static int nz = 0, noz = 0, nopar = 0;
278         double *W = NULL;
279
280         // Set FPU Flag to Extended for scicos simulation
281         // in order to override Java setting it to Double.
282 #if defined(linux) && defined(__i386__)
283         setFPUToExtended();
284 #endif
285         /*     Copyright INRIA */
286         /* iz,izptr are used to pass block labels */
287         TCritWarning = 0;
288
289         t0=t0_in;
290         tf=tf_in;
291         ierr=ierr_out;
292
293         /* Parameter adjustments */
294         pointi=pointi_in;
295         x=x_in;
296         xptr=xptr_in-1;
297         modptr=modptr_in-1;
298         --zptr;
299         --izptr;
300         --ozptr;
301         evtspt=evtspt_in-1;
302         tevts=tevts_in-1;
303         outtbptr=outtbptr_in;
304         outtbsz=outtbsz_in;
305         outtbtyp=outtbtyp_in;
306         outtb_elem=outtb_elem_in;
307         funtyp=funtyp_in-1;
308         inpptr=inpptr_in-1;
309         outptr=outptr_in-1;
310         inplnk=inplnk_in-1;
311         outlnk=outlnk_in-1;
312         --rpptr;
313         --ipptr;
314         --opptr;
315         clkptr=clkptr_in-1;
316         ordptr=ordptr_in-1;
317         ordclk=ordclk_in-1;
318         cord=cord_in-1;
319         iord=iord_in-1;
320         oord=oord_in-1;
321         zord=zord_in-1;
322
323         critev=critev_in-1;
324         --ztyp;
325         zcptr=zcptr_in-1;
326         --simpar;
327
328         /* Function Body */
329         Atol = simpar[1];
330         rtol = simpar[2];
331         ttol = simpar[3];
332         deltat = simpar[4];
333         C2F(rtfactor).scale = simpar[5];
334         C2F(cmsolver).solver = (int) simpar[6];
335         hmax=simpar[7];
336
337         nordptr = *nordptr1;
338         nblk  = *nblk1;
339         ncord = *ncord1;
340         noord = *noord1;
341         nzord = *nzord1;
342         niord = *niord1;
343         nlnk  = *nlnk1;
344         nelem = *nelem1;
345         *ierr = 0;
346
347         nordclk = ordptr[nordptr]-1;    /* number of rows in ordclk is ordptr(nclkp1)-1 */
348         ng      = zcptr[nblk + 1] - 1;  /* computes number of zero crossing surfaces */
349         nmod    = modptr[nblk + 1] - 1; /* computes number of modes */
350         nz      = zptr[nblk + 1] - 1;   /* number of discrete real states */
351         noz     = ozptr[nblk + 1] - 1;  /* number of discrete object states */
352         nopar   = opptr[nblk + 1] - 1;  /* number of object parameters */
353         nx      = xptr[nblk +1] - 1;    /* number of object parameters */
354         neq     = &nx;
355
356         xd = &x[xptr[nblk+1]-1];
357
358         /* check for hard coded maxsize */
359         for (i = 1; i <= nblk; ++i) {
360                 if (funtyp[i] < 10000) {
361                         funtyp[i] %= 1000;
362                 } else {
363                         funtyp[i] = funtyp[i] % 1000 + 10000;
364                 }
365                 ni = inpptr[i + 1] - inpptr[i];
366                 no = outptr[i + 1] - outptr[i];
367                 if (funtyp[i] == 1) {
368                         if (ni + no > 11) {
369                                 /*     hard coded maxsize in callf.c */
370                                 C2F(msgs)(&c__90, &c__0);
371                                 C2F(curblk).kfun = i;
372                                 *ierr = i + 1005;
373                                 return 0;
374                         }
375                 } else if (funtyp[i] == 2 || funtyp[i] == 3) {
376                         /*     hard coded maxsize in scicos.h */
377                         if (ni + no > SZ_SIZE) {
378                                 C2F(msgs)(&c__90, &c__0);
379                                 C2F(curblk).kfun = i;
380                                 *ierr = i + 1005;
381                                 return 0;
382                         }
383                 }
384                 mxtb = 0;
385                 if (funtyp[i] == 0) {
386                         if (ni > 1) {
387                                 for (j = 1; j <= ni; ++j) {
388                                         k = inplnk[inpptr[i] - 1 + j];
389                                         mxtb = mxtb + (outtbsz[k-1]*outtbsz[(k-1)+nlnk]);
390                                 }
391                         }
392                         if (no > 1) {
393                                 for (j = 1; j <= no; ++j) {
394                                         k = outlnk[outptr[i] - 1 + j];
395                                         mxtb = mxtb + (outtbsz[k-1]*outtbsz[(k-1)+nlnk]);
396                                 }
397                         }
398                         if (mxtb > TB_SIZE) {
399                                 C2F(msgs)(&c__91, &c__0);
400                                 C2F(curblk).kfun = i;
401                                 *ierr = i + 1005;
402                                 return 0;
403                         }
404                 }
405         }
406
407         if(nx>0) { /* xprop */
408                 if((xprop=MALLOC(sizeof(int)*nx))== NULL ) {
409                         *ierr =5;
410                         return 0;
411                 }
412         }
413         for(i=0;i<nx;i++) { /* initialize */
414                 xprop[i]=1;
415         }
416         if(nmod>0) { /* mod */
417                 if((mod=MALLOC(sizeof(int)*nmod))== NULL ) {
418                         *ierr =5;
419                         if(nx>0) FREE(xprop);
420                         return 0;
421                 }
422         }
423         if(ng>0) { /* g becomes global */
424                 if((g=MALLOC(sizeof(double)*ng))== NULL ) {
425                         *ierr =5;
426                         if(nmod>0) FREE(mod);
427                         if(nx>0) FREE(xprop);
428                         return 0;
429                 }
430         }
431
432         debug_block=-1; /* no debug block for start */
433         C2F(cosdebugcounter).counter=0;
434
435         /** Create Block's array **/
436         if((Blocks=MALLOC(sizeof(scicos_block)*nblk))== NULL ) {
437                 *ierr =5;
438                 if(nx>0) FREE(xprop);
439                 if(nmod>0) FREE(mod);
440                 if(ng>0) FREE(g);
441                 return 0;
442         }
443
444         /** Setting blocks properties for each entry in Block's array **/
445
446         /* 1 : type and pointer on simulation function */
447         for (kf = 0; kf < nblk; ++kf) { /*for each block */
448                 C2F(curblk).kfun = kf+1;
449                 i=funptr[kf];
450                 Blocks[kf].type=funtyp[kf+1];
451                 if (i<0) {
452                         switch (funtyp[kf+1]) {
453           case 0:
454                   Blocks[kf].funpt=F2C(sciblk);
455                   break;
456           case 1:
457                   sciprint(_("type 1 function not allowed for scilab blocks\n"));
458                   *ierr =1000+kf+1;
459                   FREE_blocks();
460                   return 0;
461           case 2:
462                   sciprint(_("type 2 function not allowed for scilab blocks\n"));
463                   *ierr =1000+kf+1;
464                   FREE_blocks();
465                   return 0;
466           case 3:
467                   Blocks[kf].funpt=sciblk2;
468                   Blocks[kf].type=2;
469                   break;
470           case 5:
471                   Blocks[kf].funpt=sciblk4;
472                   Blocks[kf].type=4;
473                   break;
474           case 99: /* debugging block */
475                   Blocks[kf].funpt=sciblk4;
476                   /*Blocks[kf].type=4;*/
477                   debug_block=kf;
478                   break;
479
480           case 10005:
481                   Blocks[kf].funpt=sciblk4;
482                   Blocks[kf].type=10004;
483                   break;
484           default :
485                   sciprint(_("Undefined Function type\n"));
486                   *ierr =1000+kf+1;
487                   FREE_blocks();
488                   return 0;
489                         }
490                         Blocks[kf].scsptr=-i; /* set scilab function adress for sciblk */
491                 }
492                 else if (i<=ntabsim){
493                         Blocks[kf].funpt=*(tabsim[i-1].fonc);
494                         Blocks[kf].scsptr=0;     /* this is done for being able to test if a block
495                                                                          is a scilab block in the debugging phase when 
496                                                                          sciblk4 is called */
497                 }
498                 else {
499                         i -= (ntabsim+1);
500                         //GetDynFunc(i,&Blocks[kf].funpt);  // TODO: old Scilab 5 function to be replaced
501                         if ( Blocks[kf].funpt == (voidf) 0) {
502                                 sciprint(_("Function not found\n"));
503                                 *ierr =1000+kf+1;
504                                 FREE_blocks();
505                                 return 0;
506                         }
507                         Blocks[kf].scsptr=0;   /* this is done for being able to test if a block
508                                                                    is a scilab block in the debugging phase when 
509                                                                    sciblk4 is called */
510                 }
511
512                 /* 2 : Dimension properties */
513                 Blocks[kf].ztyp=ztyp[kf+1];
514                 Blocks[kf].nx=xptr[kf+2]-xptr[kf+1]; /* continuuous state dimension*/
515                 Blocks[kf].ng=zcptr[kf+2]-zcptr[kf+1];/* number of zero crossing surface*/
516                 Blocks[kf].nz=zptr[kf+2]-zptr[kf+1];/* number of double discrete state*/
517                 Blocks[kf].noz=ozptr[kf+2]-ozptr[kf+1];/* number of other discrete state*/
518                 Blocks[kf].nrpar=rpptr[kf+2]-rpptr[kf+1];/* size of double precision parameter vector*/
519                 Blocks[kf].nipar=ipptr[kf+2]-ipptr[kf+1];/* size of integer precision parameter vector*/
520                 Blocks[kf].nopar=opptr[kf+2]-opptr[kf+1];/* number of other parameters (matrix, data structure,..)*/
521                 Blocks[kf].nin=inpptr[kf+2]-inpptr[kf+1]; /* number of input ports */
522                 Blocks[kf].nout=outptr[kf+2]-outptr[kf+1];/* number of output ports */
523
524                 /* 3 : input port properties */
525                 /* in insz, we store :
526                 *  - insz[0..nin-1] : first dimension of input ports
527                 *  - insz[nin..2*nin-1] : second dimension of input ports
528                 *  - insz[2*nin..3*nin-1] : type of data of input ports
529                 */
530                 /* allocate size and pointer arrays (number of input ports)*/
531                 Blocks[kf].insz=NULL;
532                 Blocks[kf].inptr=NULL;
533                 if (Blocks[kf].nin!=0) {
534                         if ((Blocks[kf].insz=MALLOC(Blocks[kf].nin*3*sizeof(int)))== NULL ){
535                                 FREE_blocks();
536                                 *ierr =5;
537                                 return 0;
538                         }
539                         if ((Blocks[kf].inptr=MALLOC(Blocks[kf].nin*sizeof(double*)))== NULL ){
540                                 FREE_blocks();
541                                 *ierr =5;
542                                 return 0;
543                         }
544                 }
545                 for(in=0;in<Blocks[kf].nin;in++) {
546                         lprt=inplnk[inpptr[kf+1]+in];
547                         Blocks[kf].inptr[in]=outtbptr[lprt-1];/* pointer on the data*/
548                         Blocks[kf].insz[in]=outtbsz[lprt-1];/* row dimension of the input port*/
549                         Blocks[kf].insz[Blocks[kf].nin+in]=outtbsz[(lprt-1)+nlnk];/* column dimension of the input port*/
550                         Blocks[kf].insz[2*Blocks[kf].nin+in]=outtbtyp[lprt-1];/*type of data of the input port*/
551                 }
552                 /* 4 : output port properties */
553                 /* in outsz, we store :
554                 *  - outsz[0..nout-1] : first dimension of output ports
555                 *  - outsz[nout..2*nout-1] : second dimension of output ports
556                 *  - outsz[2*nout..3*nout-1] : type of data of output ports
557                 */
558                 /* allocate size and pointer arrays (number of output ports)*/
559                 Blocks[kf].outsz=NULL;
560                 Blocks[kf].outptr=NULL;
561                 if (Blocks[kf].nout!=0) {
562                         if ((Blocks[kf].outsz=MALLOC(Blocks[kf].nout*3*sizeof(int)))== NULL ){
563                                 FREE_blocks();
564                                 *ierr =5;
565                                 return 0;
566                         }
567                         if ((Blocks[kf].outptr=MALLOC(Blocks[kf].nout*sizeof(double*)))== NULL ){
568                                 FREE_blocks();
569                                 *ierr =5;
570                                 return 0;
571                         }
572                 }
573                 /* set the values */
574                 for(out=0;out<Blocks[kf].nout;out++) { /*for each output port */
575                         lprt=outlnk[outptr[kf+1]+out];
576                         Blocks[kf].outptr[out]=outtbptr[lprt-1]; /*pointer on data */
577                         Blocks[kf].outsz[out]=outtbsz[lprt-1]; /*row dimension of output port*/
578                         Blocks[kf].outsz[Blocks[kf].nout+out]=outtbsz[(lprt-1)+nlnk]; /*column dimension of output ports*/
579                         Blocks[kf].outsz[2*Blocks[kf].nout+out]=outtbtyp[lprt-1];/*type of data of output port */
580                 }
581
582                 /* 5 : event output port properties */
583                 Blocks[kf].evout=NULL;
584                 Blocks[kf].nevout=clkptr[kf+2] - clkptr[kf+1];
585                 if (Blocks[kf].nevout!=0) {
586                         if ((Blocks[kf].evout=CALLOC(Blocks[kf].nevout,sizeof(double)))== NULL ){
587                                 FREE_blocks();
588                                 *ierr =5;
589                                 return 0;
590                         }
591                 }
592
593                 /* 6 : pointer on the begining of the double discrete state array ( z) */
594                 Blocks[kf].z=&(z__[zptr[kf+1]-1]);
595
596                 /* 7 : type, size and pointer on the other discrete states  data structures (oz) */
597                 Blocks[kf].ozsz=NULL;
598                 if (Blocks[kf].noz==0) {
599                         Blocks[kf].ozptr=NULL;
600                         Blocks[kf].oztyp=NULL;
601                 }
602                 else {
603                         Blocks[kf].ozptr=&(oz[ozptr[kf+1]-1]);
604                         if ((Blocks[kf].ozsz=MALLOC(Blocks[kf].noz*2*sizeof(int)))== NULL ) {
605                                 FREE_blocks();
606                                 *ierr =5;
607                                 return 0;
608                         }
609                         for (i=0;i<Blocks[kf].noz;i++) {
610                                 Blocks[kf].ozsz[i]=ozsz[(ozptr[kf+1]-1)+i];
611                                 Blocks[kf].ozsz[i+Blocks[kf].noz]=ozsz[(ozptr[kf+1]-1+noz)+i];
612                         }
613                         Blocks[kf].oztyp=&(oztyp[ozptr[kf+1]-1]);
614                 }
615
616                 /* 8 : pointer on the begining of the double parameter array ( rpar ) */
617                 Blocks[kf].rpar=&(rpar[rpptr[kf+1]-1]);
618
619                 /* 9 : pointer on the begining of the integer parameter array ( ipar ) */
620                 Blocks[kf].ipar=&(ipar[ipptr[kf+1]-1]);
621
622                 /* 10 : type, size and pointer on the other parameters  data structures (opar) */
623                 Blocks[kf].oparsz=NULL;
624                 if (Blocks[kf].nopar==0) {
625                         Blocks[kf].oparptr=NULL;
626                         Blocks[kf].opartyp=NULL;
627                 }
628                 else {
629                         Blocks[kf].oparptr=&(opar[opptr[kf+1]-1]);
630                         if ((Blocks[kf].oparsz=MALLOC(Blocks[kf].nopar*2*sizeof(int)))==NULL) {
631                                 FREE_blocks();
632                                 *ierr =5;
633                                 return 0;
634                         }
635                         for (i=0;i<Blocks[kf].nopar;i++) {
636                                 Blocks[kf].oparsz[i]=oparsz[(opptr[kf+1]-1)+i];
637                                 Blocks[kf].oparsz[i+Blocks[kf].nopar]=oparsz[(opptr[kf+1]-1+nopar)+i];
638                         }
639                         Blocks[kf].opartyp=&(opartyp[opptr[kf+1]-1]);
640                 }
641
642                 /* 10 : pointer on the beginning of the residual array (res) */
643                 Blocks[kf].res=NULL;
644                 if (Blocks[kf].nx!=0) {
645                         if ((Blocks[kf].res=MALLOC(Blocks[kf].nx*sizeof(double)))==NULL){
646                                 FREE_blocks();
647                                 *ierr =5;
648                                 return 0;
649                         }
650                 }
651
652                 /* 11 : block label (label) */
653                 i1=izptr[kf+2]-izptr[kf+1];
654                 if ((Blocks[kf].label=MALLOC(sizeof(char)*(i1+1)))==NULL){
655                         FREE_blocks();
656                         *ierr =5;
657                         return 0;
658                 }
659                 Blocks[kf].label[i1]='\0';
660                 C2F(cvstr)(&i1,&(iz[izptr[kf+1]-1]),Blocks[kf].label,&job,i1);
661
662                 /* 12 : block array of crossed surfaces (jroot) */
663                 Blocks[kf].jroot=NULL;
664                 if (Blocks[kf].ng!=0) {
665                         if ((Blocks[kf].jroot=CALLOC(Blocks[kf].ng,sizeof(int)))==NULL){
666                                 FREE_blocks();
667                                 *ierr =5;
668                                 return 0;
669                         }
670                 }
671
672                 /* 13 : block work array  (work) */
673                 Blocks[kf].work=(void **)(((double *)work)+kf);
674
675                 /* 14 : block modes  array  (mode) */
676                 Blocks[kf].nmode=modptr[kf+2]-modptr[kf+1];
677                 if (Blocks[kf].nmode!=0) {
678                         Blocks[kf].mode=&(mod[modptr[kf+1]-1]);
679                 }
680
681                 /* 15 : block xprop  array  (xprop) */  
682                 Blocks[kf].xprop=NULL;
683                 if (Blocks[kf].nx!=0) {
684                         Blocks[kf].xprop=&(xprop[xptr[kf+1]-1]);
685                 }
686
687                 /* 16 : pointer on the zero crossing surface computation function of the block (g) */
688                 Blocks[kf].g=NULL;
689                 if (Blocks[kf].ng!=0) {
690                         Blocks[kf].g=&(g[zcptr[kf+1]-1]);
691                 }
692         }
693         /** all block properties are stored in the Blocks array **/
694
695         /* iwa */
696         iwa=NULL;
697         if ((*nevts)!=0) {
698                 if((iwa=MALLOC(sizeof(int)*(*nevts)))==NULL) {
699                         FREE_blocks();
700                         *ierr =5;
701                         return 0;
702                 }
703         }
704
705         /* save ptr of scicos in import structure */
706         C2F(makescicosimport)(x,&nx,&xptr[1],&zcptr[1],z__,&nz,&zptr[1],
707                 &noz,oz,ozsz,oztyp,&ozptr[1],
708                 g,&ng,mod,&nmod,&modptr[1],iz,&izptr[1],
709                 &inpptr[1],&inplnk[1],&outptr[1],&outlnk[1],
710                 outtbptr,outtbsz,outtbtyp,
711                 outtb_elem,&nelem,
712                 &nlnk,rpar,&rpptr[1],ipar,&ipptr[1],
713                 opar,oparsz,opartyp,&opptr[1],
714                 &nblk,subscr,nsubs,
715                 &tevts[1],&evtspt[1],nevts,pointi,
716                 &iord[1],&niord,&oord[1],&noord,&zord[1],&nzord,
717                 funptr,&funtyp[1],&ztyp[1],
718                 &cord[1],&ncord,&ordclk[1],&nordclk,&clkptr[1],
719                 &ordptr[1],&nordptr,&critev[1],iwa,Blocks,
720                 t0,tf,&Atol,&rtol,&ttol,&deltat,&hmax,
721                 xprop,xd);
722
723         if (*flag__ == 1) { /*start*/
724                 /*     initialisation des blocks */
725                 for (kf = 0; kf < nblk; ++kf) {
726                         *(Blocks[kf].work)=NULL;
727                 }
728                 cosini(t0);
729                 if (*ierr != 0) {
730                         ierr0=*ierr;
731                         kfun0 = C2F(curblk).kfun;
732                         cosend(t0);
733                         *ierr=ierr0;
734                         C2F(curblk).kfun = kfun0;
735                 }
736
737         } else if (*flag__ == 2) { /*run*/
738
739                 /*     integration */
740                 if (C2F(cmsolver).solver == 0) {      /*  CVODE: Method: BDF,   Nonlinear solver= NEWTON     */
741                         cossim(t0);
742                 }else if (C2F(cmsolver).solver == 1) {/*  CVODE: Method: BDF,   Nonlinear solver= FUNCTIONAL */
743                         cossim(t0);
744                 }else if (C2F(cmsolver).solver == 2) {/*  CVODE: Method: ADAMS, Nonlinear solver= NEWTON     */
745                         cossim(t0);
746                 }else if (C2F(cmsolver).solver == 3) {/*  CVODE: Method: ADAMS, Nonlinear solver= FUNCTIONAL */
747                         cossim(t0);
748                 } else if(C2F(cmsolver).solver == 100){/* IDA  : Method:       , Nonlinear solver=  */
749                         cossimdaskr(t0);
750                 } else {
751                         /*     add a warning message please */
752                 }
753                 if (*ierr != 0) {
754                         ierr0=*ierr;
755                         kfun0 = C2F(curblk).kfun;
756                         cosend(t0);
757                         *ierr=ierr0;
758                         C2F(curblk).kfun = kfun0;
759                 }
760
761         } else if (*flag__ == 3) { /*finish*/
762                 /*     fermeture des blocks */
763                 cosend(t0);
764         } else if (*flag__ == 4) { /*linear*/
765                 phase=1;
766                 idoit(t0);
767                 if (*ierr == 0) {
768                         if((W=MALLOC(sizeof(double)*(max(nx,ng))))== NULL ){
769                                 FREE(iwa);
770                                 FREE_blocks();
771                                 *ierr =5;
772                                 return 0;
773                         }
774
775                         /*---------à la place de old simblk--------*/
776                         /*  C2F(simblk)(&nx, t0, x, W);  */
777
778                         if (ng>0&&nmod>0){
779                                 zdoit(t0, x, x+nx, W); /* updating modes as a function of state values; this was necessary in iGUI*/
780                         }
781                         for(jj=0;jj<nx;jj++) W[jj]=0.0;
782                         C2F(ierode).iero = 0;   *ierr= 0;
783                         if(C2F(cmsolver).solver < 100){
784                                 odoit(t0, x, W, W);
785                         }
786                         else {
787                                 odoit(t0, x, x+nx, W);
788                         }
789                         C2F(ierode).iero = *ierr;
790                         /*-----------------------------------------*/
791                         for (i = 0; i < nx; ++i) {
792                                 x[i] = W[i];
793                         }
794                         FREE(W);
795                 }
796         } else if (*flag__ == 5) { /* initial_KINSOL= "Kinsol" */
797                 C2F(ierode).iero = 0;   *ierr= 0;
798                 idoit(t0);
799                 CallKinsol(t0);
800                 *ierr=C2F(ierode).iero;
801         }
802
803
804         FREE(iwa);
805         FREE_blocks();
806
807         C2F(clearscicosimport)();
808         return 0;
809 } /* scicos_ */
810 /*--------------------------------------------------------------------------*/
811 /* check_flag */
812 static int check_flag(void *flagvalue, char *funcname, int opt)
813 {
814         int *errflag = NULL;
815
816         /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
817         if (opt == 0 && flagvalue == NULL) {
818                 sciprint(_("\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n"), funcname);
819                 return(1); }
820         /* Check if flag < 0 */
821         else if (opt == 1) {
822                 errflag = (int *) flagvalue;
823                 if (*errflag < 0) {
824                         sciprint(_("\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n"),
825                                 funcname, *errflag);
826                         return(1); }}
827         /* Check if function returned NULL pointer - no memory allocated */
828         else if (opt == 2 && flagvalue == NULL) {
829                 sciprint(_("\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n"),funcname);
830                 return(1); }
831
832         return(0);
833 } /* check_flag */
834
835 /*--------------------------------------------------------------------------*/
836 static void cosini(double *told)
837 {
838         static scicos_flag flag__ = 0;
839         static int i = 0;
840
841         static int kfune = 0;
842         static int jj = 0;
843
844         SCSREAL_COP *outtbd = NULL;    /*to save double of outtb*/
845         SCSINT8_COP *outtbc = NULL;    /*to save int8 of outtb*/
846         SCSINT16_COP *outtbs = NULL;   /*to save int16 of outtb*/
847         SCSINT32_COP *outtbl = NULL;   /*to save int32 of outtb*/
848         SCSUINT8_COP *outtbuc = NULL;  /*to save unsigned int8 of outtb*/
849         SCSUINT16_COP *outtbus = NULL; /*to save unsigned int16 of outtb*/
850         SCSUINT32_COP *outtbul=NULL; /*to save unsigned int32 of outtb*/
851         int szouttbd = 0;  /*size of arrays*/
852         int szouttbc = 0,  szouttbs = 0,  szouttbl = 0;
853         int szouttbuc = 0, szouttbus = 0, szouttbul = 0;
854         int curouttbd = 0; /*current position in arrays*/
855         int curouttbc = 0,  curouttbs = 0,  curouttbl = 0;
856         int curouttbuc = 0, curouttbus = 0, curouttbul = 0;
857
858         int ii = 0, kk = 0; /*local counters*/
859         int sszz = 0;  /*local size of element of outtb*/
860         /*Allocation of arrays for outtb*/
861         for (ii=0;ii<nlnk;ii++)
862         {
863                 switch (outtbtyp[ii])
864                 {
865                 case SCSREAL_N    : szouttbd+=outtbsz[ii]*outtbsz[ii+nlnk]; /*double real matrix*/
866                         outtbd=(SCSREAL_COP *) REALLOC (outtbd,szouttbd*sizeof(SCSREAL_COP));
867                         break;
868
869                 case SCSCOMPLEX_N : szouttbd+=2*outtbsz[ii]*outtbsz[ii+nlnk]; /*double complex matrix*/
870                         outtbd=(SCSCOMPLEX_COP *) REALLOC (outtbd,szouttbd*sizeof(SCSCOMPLEX_COP));
871                         break;
872
873                 case SCSINT8_N    : szouttbc+=outtbsz[ii]*outtbsz[ii+nlnk]; /*int8*/
874                         outtbc=(SCSINT8_COP *) REALLOC (outtbc,szouttbc*sizeof(SCSINT8_COP));
875                         break;
876
877                 case SCSINT16_N   : szouttbs+=outtbsz[ii]*outtbsz[ii+nlnk]; /*int16*/
878                         outtbs=(SCSINT16_COP *) REALLOC (outtbs,szouttbs*sizeof(SCSINT16_COP));
879                         break;
880
881                 case SCSINT32_N   : szouttbl+=outtbsz[ii]*outtbsz[ii+nlnk]; /*int32*/
882                         outtbl=(SCSINT32_COP *) REALLOC (outtbl,szouttbl*sizeof(SCSINT32_COP));
883                         break;
884
885                 case SCSUINT8_N   : szouttbuc+=outtbsz[ii]*outtbsz[ii+nlnk]; /*uint8*/
886                         outtbuc=(SCSUINT8_COP *) REALLOC (outtbuc,szouttbuc*sizeof(SCSUINT8_COP));
887                         break;
888
889                 case SCSUINT16_N  : szouttbus+=outtbsz[ii]*outtbsz[ii+nlnk]; /*uint16*/
890                         outtbus=(SCSUINT16_COP *) REALLOC (outtbus,szouttbus*sizeof(SCSUINT16_COP));
891                         break;
892
893                 case SCSUINT32_N  : szouttbul+=outtbsz[ii]*outtbsz[ii+nlnk]; /*uint32*/
894                         outtbul=(SCSUINT32_COP *) REALLOC (outtbul,szouttbul*sizeof(SCSUINT32_COP));
895                         break;
896
897                 default  : /* Add a message here */
898                         break;
899                 }
900         }
901
902         /* Jacobian*/
903         AJacobian_block = 0;
904
905         /* Function Body */
906         *ierr = 0;
907
908         /*     initialization (flag 4) */
909         /*     loop on blocks */
910         C2F(dset)(&ng, &c_b14, g, &c__1);
911
912         for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun) {
913                 flag__ = 4;
914                 if(Blocks[C2F(curblk).kfun-1].nx > 0) {
915                         Blocks[C2F(curblk).kfun-1].x  = &x[xptr[C2F(curblk).kfun]-1];
916                         Blocks[C2F(curblk).kfun-1].xd = &xd[xptr[C2F(curblk).kfun]-1];
917                 }
918                 Blocks[C2F(curblk).kfun-1].nevprt = 0;
919                 if (funtyp[C2F(curblk).kfun] >= 0) { /* debug_block is not called here */
920                         /*callf(told, xd, x, x,g,&flag__);*/
921                         Jacobian_Flag=0;
922                         callf(told, &Blocks[C2F(curblk).kfun-1], &flag__);
923                         if (flag__ < 0 && *ierr == 0) {
924                                 *ierr = 5 - flag__;
925                                 kfune = C2F(curblk).kfun;
926                         }
927                         if ((Jacobian_Flag==1)&&(AJacobian_block==0)) AJacobian_block=C2F(curblk).kfun;
928                 }
929         }
930         if (*ierr != 0) {
931                 C2F(curblk).kfun = kfune;
932                 freeouttbptr;
933                 return;
934         }
935
936         /*     initialization (flag 6) */
937         flag__ = 6;
938         for (jj = 1; jj <= ncord; ++jj) {
939                 C2F(curblk).kfun = cord[jj];
940                 Blocks[C2F(curblk).kfun-1].nevprt = 0;
941                 if (funtyp[C2F(curblk).kfun] >= 0) {
942                         /*callf(told, xd, x, x,g,&flag__);*/
943                         callf(told, &Blocks[C2F(curblk).kfun-1], &flag__);
944                         if (flag__ < 0) {
945                                 *ierr = 5 - flag__;
946                                 freeouttbptr;
947                                 return;
948                         }
949                 }
950         }
951
952         /*     point-fix iterations */
953         flag__ = 6;
954         for (i = 1; i <= nblk + 1; ++i) { /*for each block*/
955                 /*     loop on blocks */
956                 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun) {
957                         Blocks[C2F(curblk).kfun-1].nevprt = 0;
958                         if (funtyp[C2F(curblk).kfun] >= 0) {
959                                 /*callf(told, xd, x, x,g,&flag__);*/
960                                 callf(told, &Blocks[C2F(curblk).kfun-1], &flag__);
961                                 if (flag__ < 0) {
962                                         *ierr = 5 - flag__;
963                                         freeouttbptr;
964                                         return;
965                                 }
966                         }
967                 }
968
969                 flag__ = 6;
970                 for (jj = 1; jj <= ncord; ++jj) { /*for each continous block*/
971                         C2F(curblk).kfun = cord[jj];
972                         if (funtyp[C2F(curblk).kfun] >= 0) {
973                                 /*callf(told, xd, x, x,g,&flag__);*/
974                                 callf(told, &Blocks[C2F(curblk).kfun-1], &flag__);
975                                 if (flag__ < 0) {
976                                         *ierr = 5 - flag__;
977                                         freeouttbptr;
978                                         return;
979                                 }
980                         }
981                 }
982
983                 /*comparison between outtb and arrays*/
984                 curouttbd = 0;  curouttbc = 0;  curouttbs = 0; curouttbl = 0;
985                 curouttbuc = 0; curouttbus = 0; curouttbul = 0;
986                 for (jj = 0; jj<nlnk; jj++)
987                 {
988                         switch (outtbtyp[jj]) /*for each type of ports*/
989                         {
990                         case SCSREAL_N    : outtbdptr=(SCSREAL_COP *)outtbptr[jj]; /*double real matrix*/
991                                 sszz=outtbsz[jj]*outtbsz[jj+nlnk];
992                                 for(kk=0;kk<sszz;kk++)
993                                 {
994                                         if(outtbdptr[kk]!=(SCSREAL_COP)outtbd[curouttbd+kk]) goto L30;
995                                 }
996                                 curouttbd+=sszz;
997                                 break;
998
999                         case SCSCOMPLEX_N : outtbdptr=(SCSCOMPLEX_COP *)outtbptr[jj]; /*double complex matrix*/
1000                                 sszz=2*outtbsz[jj]*outtbsz[jj+nlnk];
1001                                 for(kk=0;kk<sszz;kk++)
1002                                 {
1003                                         if(outtbdptr[kk]!=(SCSCOMPLEX_COP)outtbd[curouttbd+kk]) goto L30;
1004                                 }
1005                                 curouttbd+=sszz;
1006                                 break;
1007
1008                         case SCSINT8_N    : outtbcptr=(SCSINT8_COP *)outtbptr[jj]; /*int8*/
1009                                 sszz=outtbsz[jj]*outtbsz[jj+nlnk];
1010                                 for(kk=0;kk<sszz;kk++)
1011                                 {
1012                                         if(outtbcptr[kk]!=(SCSINT8_COP)outtbc[curouttbc+kk]) goto L30;
1013                                 }
1014                                 curouttbc+=sszz;
1015                                 break;
1016
1017                         case SCSINT16_N   : outtbsptr=(SCSINT16_COP *)outtbptr[jj]; /*int16*/
1018                                 sszz=outtbsz[jj]*outtbsz[jj+nlnk];
1019                                 for (kk=0;kk<sszz;kk++)
1020                                 {
1021                                         if(outtbsptr[kk]!=(SCSINT16_COP)outtbs[curouttbs+kk]) goto L30;
1022                                 }
1023                                 curouttbs+=sszz;
1024                                 break;
1025
1026                         case SCSINT32_N   : outtblptr=(SCSINT32_COP *)outtbptr[jj]; /*int32*/
1027                                 sszz=outtbsz[jj]*outtbsz[jj+nlnk];
1028                                 for (kk=0;kk<sszz;kk++)
1029                                 {
1030                                         if(outtblptr[kk]!=(SCSINT32_COP)outtbl[curouttbl+kk]) goto L30;
1031                                 }
1032                                 curouttbl+=sszz;
1033                                 break;
1034
1035                         case SCSUINT8_N   : outtbucptr=(SCSUINT8_COP *)outtbptr[jj]; /*uint8*/
1036                                 sszz=outtbsz[jj]*outtbsz[jj+nlnk];
1037                                 for (kk=0;kk<sszz;kk++)
1038                                 {
1039                                         if(outtbucptr[kk]!=(SCSUINT8_COP)outtbuc[curouttbuc+kk]) goto L30;
1040                                 }
1041                                 curouttbuc+=sszz;
1042                                 break;
1043
1044                         case SCSUINT16_N  : outtbusptr=(SCSUINT16_COP *)outtbptr[jj]; /*uint16*/
1045                                 sszz=outtbsz[jj]*outtbsz[jj+nlnk];
1046                                 for (kk=0;kk<sszz;kk++)
1047                                 {
1048                                         if(outtbusptr[kk]!=(SCSUINT16_COP)outtbus[curouttbus+kk]) goto L30;
1049                                 }
1050                                 curouttbus+=sszz;
1051                                 break;
1052
1053                         case SCSUINT32_N  : outtbulptr=(SCSUINT32_COP *)outtbptr[jj]; /*uint32*/
1054                                 sszz=outtbsz[jj]*outtbsz[jj+nlnk];
1055                                 for (kk=0;kk<sszz;kk++)
1056                                 {
1057                                         if(outtbulptr[kk]!=(SCSUINT32_COP)outtbul[curouttbul+kk]) goto L30;
1058                                 }
1059                                 curouttbul+=sszz;
1060                                 break;
1061
1062                         default  : /* Add a message here */
1063                                 break;
1064                         }
1065                 }
1066                 freeouttbptr;
1067                 return;
1068
1069 L30:
1070                 /*Save data of outtb in arrays*/
1071                 curouttbd = 0;
1072                 curouttbc = 0;  curouttbs = 0;  curouttbl = 0;
1073                 curouttbuc = 0; curouttbus = 0; curouttbul = 0;
1074                 for (ii=0;ii<nlnk;ii++) /*for each link*/
1075                 {
1076                         switch (outtbtyp[ii])  /*switch to type of outtb object*/
1077                         {
1078                         case SCSREAL_N    : outtbdptr=(SCSREAL_COP *)outtbptr[ii];  /*double real matrix*/
1079                                 sszz=outtbsz[ii]*outtbsz[ii+nlnk];
1080                                 C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1081                                 curouttbd+=sszz;
1082                                 break;
1083
1084                         case SCSCOMPLEX_N : outtbdptr=(SCSCOMPLEX_COP *)outtbptr[ii];  /*double complex matrix*/
1085                                 sszz=2*outtbsz[ii]*outtbsz[ii+nlnk];
1086                                 C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1087                                 curouttbd+=sszz;
1088                                 break;
1089
1090                         case SCSINT8_N    : outtbcptr=(SCSINT8_COP *)outtbptr[ii];    /*int8*/
1091                                 sszz=outtbsz[ii]*outtbsz[ii+nlnk];
1092                                 for (kk=0;kk<sszz;kk++) outtbc[curouttbc+kk]=(SCSINT8_COP)outtbcptr[kk];
1093                                 curouttbc+=sszz;
1094                                 break;
1095
1096                         case SCSINT16_N   : outtbsptr=(SCSINT16_COP *)outtbptr[ii];   /*int16*/
1097                                 sszz=outtbsz[ii]*outtbsz[ii+nlnk];
1098                                 for (kk=0;kk<sszz;kk++) outtbs[curouttbs+kk]=(SCSINT16_COP)outtbsptr[kk];
1099                                 curouttbs+=sszz;
1100                                 break;
1101
1102                         case SCSINT32_N   : outtblptr=(SCSINT32_COP *)outtbptr[ii];    /*int32*/
1103                                 sszz=outtbsz[ii]*outtbsz[ii+nlnk];
1104                                 for (kk=0;kk<sszz;kk++) outtbl[curouttbl+kk]=(SCSINT32_COP)outtblptr[kk];
1105                                 curouttbl+=sszz;
1106                                 break;
1107
1108                         case SCSUINT8_N   : outtbucptr=(SCSUINT8_COP *)outtbptr[ii];  /*uint8*/
1109                                 sszz=outtbsz[ii]*outtbsz[ii+nlnk];
1110                                 for (kk=0;kk<sszz;kk++) outtbuc[curouttbuc+kk]=(SCSUINT8_COP)outtbucptr[kk];
1111                                 curouttbuc+=sszz;
1112                                 break;
1113
1114                         case SCSUINT16_N  : outtbusptr=(SCSUINT16_COP *)outtbptr[ii]; /*uint16*/
1115                                 sszz=outtbsz[ii]*outtbsz[ii+nlnk];
1116                                 for (kk=0;kk<sszz;kk++) outtbus[curouttbus+kk]=(SCSUINT16_COP)outtbusptr[kk];
1117                                 curouttbus+=sszz;
1118                                 break;
1119
1120                         case SCSUINT32_N  : outtbulptr=(SCSUINT32_COP *)outtbptr[ii];  /*uint32*/
1121                                 sszz=outtbsz[ii]*outtbsz[ii+nlnk];
1122                                 for (kk=0;kk<sszz;kk++) outtbul[curouttbul+kk]=(SCSUINT32_COP)outtbulptr[kk];
1123                                 curouttbul+=sszz;
1124                                 break;
1125
1126                         default  : /* Add a message here */
1127                                 break;
1128                         }
1129                 }
1130         }
1131         *ierr = 20;
1132         freeouttbptr;
1133 } /* cosini_ */
1134
1135 /*--------------------------------------------------------------------------*/
1136 static void cossim(double *told)
1137 {
1138         /* System generated locals */
1139         int i3 = 0;
1140
1141         //** used for the [stop] button
1142         static char CommandToUnstack[1024];
1143         static int CommandLength = 0;
1144         static int SeqSync = 0;
1145         static int zero = 0;
1146         static int one = 1;
1147
1148         /* Local variables */
1149         static scicos_flag flag__ = 0;
1150         static int ierr1 = 0;
1151         static int j = 0, k = 0;
1152         static double t = 0.;
1153         static int jj = 0;
1154         static double rhotmp = 0., tstop = 0.;
1155         static int inxsci = 0;
1156         static int kpo = 0, kev = 0;
1157         int Discrete_Jump = 0;
1158         int *jroot = NULL, *zcros = NULL;
1159         realtype reltol = 0., abstol = 0.;
1160         N_Vector y=NULL;
1161         void *cvode_mem = NULL;
1162         int flag = 0, flagr = 0;
1163         int cnt=0;
1164         jroot=NULL;
1165         if (ng!=0) {
1166                 if((jroot=MALLOC(sizeof(int)*ng))== NULL ){
1167                         *ierr =10000;
1168                         return;
1169                 }
1170         }
1171
1172         for ( jj = 0 ; jj < ng ; jj++ )
1173                 jroot[jj] = 0 ;
1174
1175         zcros=NULL;
1176         if (ng!=0) {
1177                 if((zcros=MALLOC(sizeof(int)*ng))== NULL ){
1178                         *ierr=10000;
1179                         if (ng>0) FREE(jroot);
1180                         return;
1181                 }
1182         }
1183
1184         reltol = (realtype) rtol;
1185         abstol = (realtype) Atol;  /* Ith(abstol,1) = realtype) Atol;*/
1186
1187         if (*neq>0) { /* Unfortunately CVODE does not work with NEQ==0 */
1188                 y = N_VNewEmpty_Serial(*neq); 
1189                 if (check_flag((void *)y, "N_VNewEmpty_Serial", 0)) {
1190                         *ierr=10000;
1191                         if (ng>0) FREE(jroot);
1192                         if (ng>0) FREE(zcros);
1193                         return;
1194                 }
1195
1196                 NV_DATA_S(y)=x;
1197
1198                 cvode_mem = NULL;
1199
1200                 /* Set extension of Sundials for scicos */
1201                 set_sundials_with_extension(TRUE);
1202
1203                 switch (C2F(cmsolver).solver)
1204                 {
1205                 case 0:   cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);break;
1206                 case 1:   cvode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);break;
1207                 case 2:   cvode_mem = CVodeCreate(CV_ADAMS, CV_NEWTON);break;
1208                 case 3:   cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);break;
1209                 }
1210
1211                 /*    cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);*/
1212
1213                 if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) {
1214                         *ierr=10000;
1215                         N_VDestroy_Serial(y); FREE(jroot); FREE(zcros);
1216                         return;
1217                 }
1218
1219                 flag = CVodeMalloc(cvode_mem, simblk, T0, y, CV_SS, reltol, &abstol);
1220                 if (check_flag(&flag, "CVodeMalloc", 1)) {
1221                         *ierr=300+(-flag);
1222                         freeall
1223                                 return;
1224                 }
1225
1226                 flag = CVodeRootInit(cvode_mem, ng, grblk, NULL);
1227                 if (check_flag(&flag, "CVodeRootInit", 1)) {
1228                         *ierr=300+(-flag);
1229                         freeall
1230                                 return;
1231                 }
1232
1233                 /* Call CVDense to specify the CVDENSE dense linear solver */
1234                 flag = CVDense(cvode_mem, *neq);
1235                 if (check_flag(&flag, "CVDense", 1)) {
1236                         *ierr=300+(-flag);
1237                         freeall
1238                                 return;
1239                 }
1240
1241                 if(hmax>0){
1242                         flag=CVodeSetMaxStep(cvode_mem, (realtype) hmax);
1243                         if (check_flag(&flag, "CVodeSetMaxStep", 1)) {
1244                                 *ierr=300+(-flag);
1245                                 freeall;
1246                                 return;
1247                         }
1248                 }
1249                 /* Set the Jacobian routine to Jac (user-supplied) 
1250                 flag = CVDenseSetJacFn(cvode_mem, Jac, NULL);
1251                 if (check_flag(&flag, "CVDenseSetJacFn", 1)) return(1);  */
1252
1253         }/* testing if neq>0 */
1254
1255         /* Function Body */
1256         C2F(coshlt).halt = 0;
1257         *ierr = 0;
1258
1259         C2F(xscion)(&inxsci);
1260         /*     initialization */
1261         C2F(realtimeinit)(told, &C2F(rtfactor).scale);
1262
1263         phase=1;
1264         hot = 0;
1265
1266         jj = 0;
1267         for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun) {
1268                 if (Blocks[C2F(curblk).kfun-1].ng >= 1) {
1269                         zcros[jj] = C2F(curblk).kfun;
1270                         ++jj;
1271                 }
1272         }
1273         /*     . Il faut:  ng >= jj */
1274         if (jj != ng) {
1275                 zcros[jj] = -1;
1276         }
1277         /*     initialisation (propagation of constant blocks outputs) */
1278         idoit(told);
1279         if (*ierr != 0) {
1280                 freeall;
1281                 return;
1282         }
1283         /*--discrete zero crossings----dzero--------------------*/
1284         if (ng>0){ /* storing ZC signs just after a solver call*/
1285                 /*zdoit(told, g, x, x);*/
1286                 zdoit(told, x, x, g);
1287                 if (*ierr != 0) {freeall;return;  }
1288                 for (jj = 0; jj < ng; ++jj)
1289                         if(g[jj]>=0) jroot[jj]=5;else jroot[jj]=-5;
1290         }
1291         /*--discrete zero crossings----dzero--------------------*/
1292
1293         /*     main loop on time */
1294         while(*told < *tf) {
1295                 while (ismenu()) //** if the user has done something, do the actions
1296                 {
1297                         int ierr2=0;
1298                         SeqSync = GetCommand(CommandToUnstack); //** get at the action
1299                         CommandLength = (int)strlen(CommandToUnstack);
1300                         syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
1301                 }
1302                 if (C2F(coshlt).halt != 0) {
1303                         if (C2F(coshlt).halt ==2) *told=*tf; /* end simulation */
1304                         C2F(coshlt).halt = 0;
1305                         freeall;
1306                         return;
1307                 }
1308                 if (*pointi == 0) {
1309                         t = *tf;
1310                 } else {
1311                         t = tevts[*pointi];
1312                 }
1313                 if (abs(t - *told) < ttol) {
1314                         t = *told;
1315                         /*     update output part */
1316                 }
1317                 if (*told > t) {
1318                         /*     !  scheduling problem */
1319                         *ierr = 1;
1320                         freeall;
1321                         return;
1322                 }
1323                 if (*told != t) {
1324                         if (xptr[nblk+1] == 1) {
1325                                 /*     .     no continuous state */
1326                                 if (*told + deltat + ttol > t) {
1327                                         *told = t;
1328                                 } else {
1329                                         *told += deltat;
1330                                 }
1331                                 /*     .     update outputs of 'c' type blocks with no continuous state */
1332                                 if (*told >= *tf) {
1333                                         /*     .     we are at the end, update continuous part before leaving */
1334                                         if (ncord > 0) {
1335                                                 cdoit(told);
1336                                                 freeall;
1337                                                 return;
1338                                         }
1339                                 }
1340                         } else {
1341                                 /*     integrate */
1342                                 rhotmp = *tf + ttol;
1343                                 if (*pointi!=0){
1344                                         kpo = *pointi;
1345 L20:
1346                                         if (critev[kpo] == 1) {
1347                                                 rhotmp = tevts[kpo];
1348                                                 goto L30;
1349                                         }
1350                                         kpo = evtspt[kpo];
1351                                         if (kpo != 0) {
1352                                                 goto L20;
1353                                         }
1354 L30:
1355                                         if (rhotmp < tstop) {
1356                                                 hot = 0;
1357                                         }
1358                                 }
1359                                 tstop = rhotmp;
1360                                 t = min(*told + deltat,min(t,*tf + ttol));
1361
1362                                 if (ng>0 &&  hot == 0 && nmod>0) {
1363                                         zdoit(told,x,x,g);
1364                                         if (*ierr != 0){
1365                                                 freeall;
1366                                                 return;
1367                                         }
1368                                 }
1369
1370                                 if (hot==0){ /* hot==0 : cold restart*/
1371                                         flag = CVodeSetStopTime(cvode_mem, (realtype)tstop);  /* Setting the stop time*/
1372                                         if (check_flag(&flag, "CVodeSetStopTime", 1)) {    
1373                                                 *ierr=300+(-flag);
1374                                                 freeall;
1375                                                 return;
1376                                         }
1377
1378                                         flag =CVodeReInit(cvode_mem, simblk, (realtype)(*told), y, CV_SS, reltol, &abstol);
1379                                         if (check_flag(&flag, "CVodeReInit", 1)) {
1380                                                 *ierr=300+(-flag);
1381                                                 freeall;
1382                                                 return;
1383                                         }
1384                                 }
1385
1386                                 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3)) {
1387                                         sciprint(_("****SUNDIALS.Cvode from: %f to %f hot= %d  \n"), *told,t,hot);
1388                                 }
1389
1390                                 /*--discrete zero crossings----dzero--------------------*/
1391                                 /*--check for Dzeros after Mode settings or ddoit()----*/
1392                                 Discrete_Jump=0;
1393
1394                                 if (ng>0 && hot==0){
1395                                         zdoit(told, x, x, g);
1396                                         if (*ierr != 0) {freeall;return;}
1397                                         for (jj = 0; jj < ng; ++jj) {
1398                                                 if((g[jj]>=0.0)&&(jroot[jj]==-5)) {Discrete_Jump=1;jroot[jj]=1;}
1399                                                 else if((g[jj]<0.0)&&(jroot[jj]==5)) {Discrete_Jump=1;jroot[jj]=-1;}
1400                                                 else jroot[jj]=0;
1401                                         }
1402                                 }
1403                                 /*--discrete zero crossings----dzero--------------------*/
1404
1405                                 if (Discrete_Jump==0){/* if there was a dzero, its event should be activated*/
1406                                         phase=2;
1407                                         flag = CVode(cvode_mem, t, y, told, CV_NORMAL_TSTOP);
1408                                         if (*ierr != 0) {freeall;return;}
1409                                         phase=1;
1410                                 }else{
1411                                         flag = CV_ROOT_RETURN; /* in order to handle discrete jumps */
1412                                 }
1413
1414                                 /*     .     update outputs of 'c' type  blocks if we are at the end*/
1415                                 if (*told >= *tf) {
1416                                         if (ncord > 0) {
1417                                                 cdoit(told);
1418                                                 freeall;
1419                                                 return;
1420                                         }
1421                                 }
1422
1423                                 if (flag>=0){
1424                                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1425                                                 sciprint(_("****SUNDIALS.Cvode reached: %f\n"),*told);
1426                                         hot = 1;
1427                                         cnt = 0;
1428                                 } else if ( flag==CV_TOO_MUCH_WORK ||  flag == CV_CONV_FAILURE || flag==CV_ERR_FAILURE) {  
1429                                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1430                                                 sciprint(_("****SUNDIALS.Cvode: too much work at time=%g (stiff region, change RTOL and ATOL)\n"),*told);         
1431                                         hot = 0;
1432                                         cnt++;
1433                                         if (cnt>5) {
1434                                                 *ierr=300+(-flag);
1435                                                 freeall;
1436                                                 return;
1437                                         }
1438                                 }else{
1439                                         if (flag<0) *ierr=300+(-flag);    /* raising errors due to internal errors, other wise erros due to flagr*/
1440                                         freeall;
1441                                         return;
1442                                 }
1443
1444                                 if (flag == CV_ZERO_DETACH_RETURN){hot=0;};  /* new feature of sundials, detects zero-detaching */
1445
1446                                 if (flag == CV_ROOT_RETURN) {
1447                                         /*     .        at a least one root has been found */
1448                                         hot = 0;
1449                                         if (Discrete_Jump==0){
1450                                                 flagr = CVodeGetRootInfo(cvode_mem, jroot);
1451                                                 if (check_flag(&flagr, "CVodeGetRootInfo", 1)) {
1452                                                         *ierr=300+(-flagr);
1453                                                         freeall;
1454                                                         return;
1455                                                 }
1456                                         }
1457                                         /*     .        at a least one root has been found */
1458                                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1459                                         {
1460                                                 sciprint(_("root found at t=: %f\n"),*told);
1461                                         }
1462                                         /*     .        update outputs affecting ztyp blocks ONLY FOR OLD BLOCKS */
1463                                         zdoit(told, x, xd, g);
1464                                         if (*ierr != 0) {
1465                                                 freeall;
1466                                                 return;
1467                                         }
1468                                         for (jj = 0; jj < ng; ++jj) {
1469                                                 C2F(curblk).kfun = zcros[ jj];
1470                                                 if (C2F(curblk).kfun == -1) {
1471                                                         break; 
1472                                                 }
1473                                                 kev = 0;
1474
1475                                                 for (j = zcptr[C2F(curblk).kfun]-1 ; 
1476                                                         j <zcptr[C2F(curblk).kfun+1]-1 ; ++j) {
1477                                                                 if(jroot[j]!=0){
1478                                                                         kev=1;
1479                                                                         break;
1480                                                                 }
1481                                                 }
1482                                                 /*   */
1483                                                 if (kev != 0) {
1484                                                         Blocks[C2F(curblk).kfun-1].jroot=&jroot[zcptr[C2F(curblk).kfun]-1];
1485                                                         if (funtyp[C2F(curblk).kfun] > 0) {
1486
1487                                                                 if (Blocks[C2F(curblk).kfun-1].nevout > 0) {
1488                                                                         flag__ = 3;
1489                                                                         if(Blocks[C2F(curblk).kfun-1].nx > 0) {
1490                                                                                 Blocks[C2F(curblk).kfun-1].x  = &x[xptr[C2F(curblk).kfun]-1];
1491                                                                                 Blocks[C2F(curblk).kfun-1].xd = &xd[xptr[C2F(curblk).kfun]-1];
1492                                                                         }
1493                                                                         /* call corresponding block to determine output event (kev) */
1494                                                                         Blocks[C2F(curblk).kfun-1].nevprt = -kev;
1495                                                                         /*callf(told, xd, x, x,g,&flag__);*/
1496                                                                         callf(told, &Blocks[C2F(curblk).kfun-1], &flag__);
1497                                                                         if (flag__ < 0) {
1498                                                                                 *ierr = 5 - flag__;
1499                                                                                 freeall;
1500                                                                                 return;
1501                                                                         }
1502                                                                         /*     .              update event agenda */
1503                                                                         for (k = 0; k < Blocks[C2F(curblk).kfun-1].nevout; ++k) {
1504                                                                                 if (Blocks[C2F(curblk).kfun-1].evout[k] >= 0.) {
1505                                                                                         i3 = k + clkptr[C2F(curblk).kfun] ;
1506                                                                                         addevs(Blocks[C2F(curblk).kfun-1].evout[k]+(*told), &i3, &ierr1);
1507                                                                                         if (ierr1 != 0) {
1508                                                                                                 /*     .                       nevts too small */
1509                                                                                                 *ierr = 3;
1510                                                                                                 freeall;
1511                                                                                                 return;
1512                                                                                         }
1513                                                                                 }
1514                                                                         }
1515                                                                 }
1516                                                                 /*     .              update state */
1517                                                                 if (Blocks[C2F(curblk).kfun-1].nx > 0) {
1518                                                                         /*     .              call corresponding block to update state */
1519                                                                         flag__ = 2;
1520                                                                         Blocks[C2F(curblk).kfun-1].x  = &x[xptr[C2F(curblk).kfun]-1];
1521                                                                         Blocks[C2F(curblk).kfun-1].xd = &xd[xptr[C2F(curblk).kfun]-1];
1522                                                                         Blocks[C2F(curblk).kfun-1].nevprt = -kev;
1523                                                                         /*callf(told, xd, x, x,g,&flag__);*/
1524                                                                         callf(told, &Blocks[C2F(curblk).kfun-1], &flag__);
1525
1526                                                                         if (flag__ < 0) {
1527                                                                                 *ierr = 5 - flag__;
1528                                                                                 freeall;
1529                                                                                 return;
1530                                                                         }
1531                                                                 }
1532                                                         }
1533                                                 }
1534                                         }
1535                                 }
1536                         }
1537                         /*--discrete zero crossings----dzero--------------------*/
1538                         if (ng>0) { /* storing ZC signs just after a sundials call*/
1539                                 zdoit(told, x, x, g);
1540                                 if (*ierr != 0) {
1541                                         freeall;
1542                                         return;
1543                                 }
1544                                 for (jj = 0; jj < ng; ++jj) {
1545                                         if(g[jj]>=0) {
1546                                                 jroot[jj]=5;
1547                                         }
1548                                         else {
1549                                                 jroot[jj]=-5;
1550                                         }
1551                                 }
1552                         }
1553                         /*--discrete zero crossings----dzero--------------------*/
1554
1555                         C2F(realtime)(told);
1556                 } else {
1557                         /*     .  t==told */
1558                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1559                         {
1560                                 sciprint(_("Event: %d activated at t=%f\n"),*pointi,*told);
1561                                 for(kev=0;kev<nblk;kev++){
1562                                         if (Blocks[kev].nmode>0){
1563                                                 sciprint(_("mode of block %d=%d, "),kev,Blocks[kev].mode[0]);
1564                                         }
1565                                 }
1566                                 sciprint(_("**mod**\n"));
1567                         }
1568
1569                         ddoit(told);
1570                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1571                         {
1572                                 sciprint(_("End of activation\n"));
1573                         }
1574                         if (*ierr != 0) {
1575                                 freeall;
1576                                 return;
1577                         }
1578
1579                 }
1580                 /*     end of main loop on time */
1581         }
1582         freeall;
1583 } /* cossim_ */
1584
1585 /*--------------------------------------------------------------------------*/
1586 static void cossimdaskr(double *told)
1587 {
1588         /* Initialized data */
1589         static int otimer = 0;
1590         /* System generated locals */
1591         int i3;
1592         //** used for the [stop] button
1593         static char CommandToUnstack[1024];
1594         static int CommandLength = 0;
1595         static int SeqSync = 0;
1596         static int zero = 0;
1597         static int one = 1;
1598
1599         /* Local variables */
1600         static scicos_flag flag__ = 0;
1601         static int ierr1 = 0;
1602         static int j = 0, k = 0;
1603         static double t = 0.;
1604         static int kk = 0,jj = 0, jt = 0;
1605         static int  ntimer = 0;
1606         static double rhotmp = 0.,tstop = 0.;
1607         static int inxsci = 0;
1608         static int kpo = 0, kev = 0;
1609
1610         int *jroot = NULL,*zcros = NULL;
1611         int maxord = 0;
1612         int *Mode_save = NULL;
1613         int Mode_change=0;
1614
1615         int flag = 0, flagr = 0; 
1616         N_Vector   yy = NULL, yp = NULL;
1617         realtype reltol = 0., abstol = 0.; 
1618         int Discrete_Jump = 0;
1619         N_Vector IDx = NULL;
1620         realtype *scicos_xproperty=NULL;
1621         N_Vector bidon = NULL, tempv1 = NULL, tempv2 = NULL, tempv3 = NULL;
1622         DenseMat TJacque = NULL;
1623         realtype *Jacque_col = NULL;
1624
1625         void *ida_mem = NULL;
1626         UserData data = NULL;
1627         IDAMem copy_IDA_mem = NULL;
1628         int maxnj = 0,maxnit = 0;
1629         /*-------------------- Analytical Jacobian memory allocation ----------*/
1630         int  Jn = 0, Jnx = 0, Jno = 0, Jni = 0, Jactaille = 0;
1631         double uround = 0.;
1632         int cnt = 0, N_iters = 0;
1633         maxord = 5;
1634
1635
1636         /* Set extension of Sundials for scicos */
1637         set_sundials_with_extension(TRUE);
1638
1639         // CI=1.0;   /* for function Get_Jacobian_ci */  
1640         jroot=NULL;
1641         if (ng!=0) {
1642                 if((jroot=MALLOC(sizeof(int)*ng))== NULL ){
1643                         *ierr =10000;
1644                         return;
1645                 }
1646         }
1647         for ( jj = 0 ; jj < ng ; jj++ )  jroot[jj] = 0 ;
1648
1649         zcros=NULL;
1650         if (ng!=0) {
1651                 if((zcros=MALLOC(sizeof(int)*ng))== NULL ){
1652                         *ierr =10000;
1653                         if (ng!=0) FREE(jroot);
1654                         return;
1655                 }
1656         }
1657
1658         Mode_save=NULL;
1659         if (nmod!=0) {
1660                 if((Mode_save=MALLOC(sizeof(int)*nmod))== NULL ){
1661                         *ierr =10000;
1662                         if (ng!=0) FREE(jroot);
1663                         if (ng!=0) FREE(zcros);
1664                         return;
1665                 }
1666         }
1667
1668         reltol = (realtype) rtol;
1669         abstol = (realtype) Atol;  /*  Ith(abstol,1) = realtype) Atol;*/
1670
1671         if (*neq>0) {
1672                 yy=NULL;
1673                 yy = N_VNewEmpty_Serial(*neq);
1674                 if(check_flag((void *)yy, "N_VNew_Serial", 0)){
1675                         if (ng!=0) FREE(jroot);
1676                         if (ng!=0) FREE(zcros);
1677                         if (nmod!=0) FREE(Mode_save);
1678                 }
1679                 NV_DATA_S(yy)=x;
1680
1681                 yp=NULL;
1682                 yp = N_VNewEmpty_Serial(*neq);
1683                 if(check_flag((void *)yp, "N_VNew_Serial", 0)){
1684                         if (*neq>0) N_VDestroy_Serial(yy);
1685                         if (ng!=0)  FREE(jroot);
1686                         if (ng!=0) FREE(zcros);
1687                         if (nmod!=0) FREE(Mode_save);
1688                         return;
1689                 }
1690                 NV_DATA_S(yp)=xd;
1691
1692                 IDx = NULL;
1693                 IDx = N_VNew_Serial(*neq); 
1694                 if (check_flag((void *)IDx, "N_VNew_Serial", 0)) {
1695                         *ierr=10000;
1696                         if (*neq>0) N_VDestroy_Serial(yp);
1697                         if (*neq>0) N_VDestroy_Serial(yy);
1698                         if (ng!=0)  FREE(jroot);
1699                         if (ng!=0) FREE(zcros);
1700                         if (nmod!=0) FREE(Mode_save);
1701                         return;
1702                 }
1703
1704                 /* Call IDACreate and IDAMalloc to initialize IDA memory */
1705                 ida_mem = NULL;
1706                 ida_mem = IDACreate();
1707                 if(check_flag((void *)ida_mem, "IDACreate", 0)) {      
1708                         if (*neq>0) N_VDestroy_Serial(IDx);
1709                         if (*neq>0) N_VDestroy_Serial(yp);
1710                         if (*neq>0) N_VDestroy_Serial(yy);
1711                         if (ng!=0) FREE(jroot);
1712                         if (ng!=0) FREE(zcros);
1713                         if (nmod!=0)  FREE(Mode_save);
1714                         return;
1715                 }
1716                 copy_IDA_mem= (IDAMem) ida_mem;
1717
1718                 flag = IDAMalloc(ida_mem, simblkdaskr, T0, yy, yp, IDA_SS, reltol, &abstol);
1719                 if(check_flag(&flag, "IDAMalloc", 1)){
1720                         *ierr=200+(-flag);  
1721                         if (*neq>0)IDAFree(&ida_mem);
1722                         if (*neq>0)N_VDestroy_Serial(IDx);
1723                         if (*neq>0) N_VDestroy_Serial(yp);
1724                         if (*neq>0)N_VDestroy_Serial(yy);
1725                         if (ng!=0) FREE(jroot);
1726                         if (ng!=0) FREE(zcros);
1727                         if (nmod!=0) FREE(Mode_save);
1728                         return;
1729                 }
1730
1731
1732                 flag = IDARootInit(ida_mem, ng, grblkdaskr, NULL);
1733                 if (check_flag(&flag, "IDARootInit", 1)) {
1734                         *ierr=200+(-flag);  
1735                         if (*neq>0)IDAFree(&ida_mem);
1736                         if (*neq>0)N_VDestroy_Serial(IDx);
1737                         if (*neq>0)N_VDestroy_Serial(yp);
1738                         if (*neq>0)N_VDestroy_Serial(yy);
1739                         if (ng!=0) FREE(jroot);
1740                         if (ng!=0) FREE(zcros);
1741                         if (nmod!=0) FREE(Mode_save);
1742                         return;
1743                 }
1744
1745
1746                 flag = IDADense(ida_mem, *neq);
1747                 if (check_flag(&flag, "IDADense", 1)) {
1748                         *ierr=200+(-flag);  
1749                         if (*neq>0)IDAFree(&ida_mem);
1750                         if (*neq>0)N_VDestroy_Serial(IDx);
1751                         if (*neq>0)N_VDestroy_Serial(yp);
1752                         if (*neq>0)N_VDestroy_Serial(yy);
1753                         if (ng!=0) FREE(jroot);
1754                         if (ng!=0) FREE(zcros);
1755                         if (nmod!=0) FREE(Mode_save);
1756                         return;
1757                 }
1758
1759                 data=NULL;
1760                 if ((data = (UserData) MALLOC(sizeof(*data)))==NULL){
1761                         *ierr=10000;
1762                         if (*neq>0)IDAFree(&ida_mem);
1763                         if (*neq>0)N_VDestroy_Serial(IDx);
1764                         if (*neq>0)N_VDestroy_Serial(yp);
1765                         if (*neq>0)N_VDestroy_Serial(yy);
1766                         if (ng!=0) FREE(jroot);
1767                         if (ng!=0) FREE(zcros);
1768                         if (nmod!=0) FREE(Mode_save);
1769                         return;
1770                 }
1771                 data->ida_mem = ida_mem;
1772                 data->ewt   = NULL;
1773                 data->iwork = NULL;
1774                 data->rwork = NULL;
1775                 data->gwork = NULL;
1776
1777                 data->ewt   = N_VNew_Serial(*neq);
1778                 if (check_flag((void *)data->ewt, "N_VNew_Serial", 0)) {
1779                         *ierr=200+(-flag);  
1780                         if (*neq>0)FREE(data);
1781                         if (*neq>0)IDAFree(&ida_mem);
1782                         if (*neq>0) N_VDestroy_Serial(IDx);
1783                         if (*neq>0) N_VDestroy_Serial(yp);
1784                         if (*neq>0) N_VDestroy_Serial(yy);
1785                         if (ng!=0) FREE(jroot);
1786                         if (ng!=0) FREE(zcros);
1787                         if (nmod!=0) FREE(Mode_save);
1788                         return;
1789                 }
1790                 if ( ng>0 ){
1791                         if ((data->gwork = (double *) MALLOC(ng * sizeof(double)))==NULL){
1792                                 if (*neq>0) N_VDestroy_Serial(data->ewt);
1793                                 if (*neq>0)FREE(data);
1794                                 if (*neq>0)IDAFree(&ida_mem);
1795                                 if (*neq>0)N_VDestroy_Serial(IDx);
1796                                 if (*neq>0)N_VDestroy_Serial(yp);
1797                                 if (*neq>0)N_VDestroy_Serial(yy);
1798                                 if (ng!=0) FREE(jroot);
1799                                 if (ng!=0) FREE(zcros);
1800                                 if (nmod!=0) FREE(Mode_save);
1801                                 return;
1802                         }
1803                 }
1804                 /*Jacobian_Flag=0; */
1805                 if (AJacobian_block>0){/* set by the block with A-Jac in flag-4 using Set_Jacobian_flag(1); */
1806                         Jn=*neq;
1807                         Jnx=Blocks[AJacobian_block-1].nx;
1808                         Jno=Blocks[AJacobian_block-1].nout;
1809                         Jni=Blocks[AJacobian_block-1].nin;
1810                 }else {
1811                         Jn=*neq;
1812                         Jnx=0;
1813                         Jno=0;
1814                         Jni=0;
1815                 }
1816                 Jactaille= 3*Jn+(Jn+Jni)*(Jn+Jno)+Jnx*(Jni+2*Jn+Jno)+(Jn-Jnx)*(2*(Jn-Jnx)+Jno+Jni)+2*Jni*Jno;    
1817
1818                 if ((data->rwork = (double *) MALLOC(Jactaille * sizeof(double)))==NULL){
1819                         if ( ng>0 ) FREE(data->gwork);
1820                         if (*neq>0)N_VDestroy_Serial(data->ewt);
1821                         if (*neq>0)FREE(data);
1822                         if (*neq>0)IDAFree(&ida_mem);
1823                         if (*neq>0)N_VDestroy_Serial(IDx);
1824                         if (*neq>0)N_VDestroy_Serial(yp);
1825                         if (*neq>0)N_VDestroy_Serial(yy);
1826                         if (ng!=0) FREE(jroot);
1827                         if (ng!=0) FREE(zcros);
1828                         if (nmod!=0) FREE(Mode_save);
1829                         *ierr =10000;
1830                         return;
1831                 }
1832
1833                 flag = IDADenseSetJacFn(ida_mem, Jacobians, data);
1834                 if(check_flag(&flag, "IDADenseSetJacFn", 1)) {
1835                         *ierr=200+(-flag);
1836                         freeallx
1837                                 return;
1838                 }
1839
1840                 TJacque= (DenseMat) DenseAllocMat(*neq, *neq);
1841
1842                 flag = IDASetRdata(ida_mem, data);
1843                 if (check_flag(&flag, "IDASetRdata", 1)) {
1844                         *ierr=200+(-flag);  
1845                         freeallx
1846                                 return;
1847                 }
1848
1849                 if(hmax>0){
1850                         flag=IDASetMaxStep(ida_mem, (realtype) hmax);
1851                         if (check_flag(&flag, "IDASetMaxStep", 1)) {
1852                                 *ierr=200+(-flag);
1853                                 freeallx
1854                                         return;
1855                         }
1856                 }
1857
1858                 maxnj=100; /* setting the maximum number of Jacobian evaluation during a Newton step */
1859                 flag=IDASetMaxNumJacsIC(ida_mem, maxnj);
1860                 if (check_flag(&flag, "IDASetMaxNumJacsIC", 1)) {
1861                         *ierr=200+(-flag);
1862                         freeallx
1863                                 return;
1864                 }
1865
1866                 maxnit=10; /* setting the maximum number of Newton iterations in any one attemp to solve CIC */
1867                 flag=IDASetMaxNumItersIC(ida_mem, maxnit);
1868                 if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
1869                         *ierr=200+(-flag);
1870                         freeallx
1871                                 return;
1872                 }
1873
1874                 /* setting the maximum number of steps in an integration interval */
1875                 flag=IDASetMaxNumSteps(ida_mem, 2000);
1876                 if (check_flag(&flag, "IDASetMaxNumSteps", 1)) {
1877                         *ierr=200+(-flag);
1878                         freeallx
1879                                 return;
1880                 }
1881
1882         } /* testing if neq>0 */
1883
1884         uround = 1.0;
1885         do{
1886                 uround = uround*0.5;
1887         }while ( 1.0 + uround != 1.0);
1888         uround = uround*2.0;
1889         SQuround=sqrt(uround);
1890         /* Function Body */
1891
1892         C2F(coshlt).halt = 0;
1893         *ierr = 0;
1894         /*     hot = .false. */
1895         phase=1;
1896         hot = 0;
1897
1898         jt = 2;
1899
1900         /*      stuck=.false. */
1901         C2F(xscion)(&inxsci);
1902         /*     initialization */
1903         C2F(realtimeinit)(told, &C2F(rtfactor).scale);
1904         /*     ATOL and RTOL are scalars */
1905
1906         jj = 0;
1907         for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun) {
1908                 if (Blocks[C2F(curblk).kfun-1].ng >= 1) {
1909                         zcros[jj] = C2F(curblk).kfun;
1910                         ++jj;
1911                 }
1912         }
1913         /*     . Il faut:  ng >= jj */
1914         if (jj != ng) {
1915                 zcros[jj] = -1;
1916         }
1917         /*     initialisation (propagation of constant blocks outputs) */
1918         idoit(told);
1919         if (*ierr != 0) {
1920                 freeallx;
1921                 return;
1922         }
1923
1924         /*--discrete zero crossings----dzero--------------------*/
1925         if (ng>0){ /* storing ZC signs just after a solver call*/
1926                 zdoit(told, x, x, g);
1927                 if (*ierr != 0) {freeallx;return;  }
1928                 for (jj = 0; jj < ng; ++jj)
1929                         if(g[jj]>=0) jroot[jj]=5;else jroot[jj]=-5;
1930         }
1931         /*     main loop on time */
1932         while (*told < *tf) {
1933                 while (ismenu()) //** if the user has done something, do the actions
1934                 {
1935                         int ierr2=0;
1936                         SeqSync = GetCommand(CommandToUnstack); //** get at the action
1937                         CommandLength = (int)strlen(CommandToUnstack);
1938                         syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
1939                 }
1940                 if (C2F(coshlt).halt != 0) {
1941                         if (C2F(coshlt).halt ==2) *told=*tf; /* end simulation */
1942                         C2F(coshlt).halt = 0;
1943                         freeallx;
1944                         return;
1945                 }
1946                 if (*pointi == 0) {
1947                         t = *tf;
1948                 } else {
1949                         t = tevts[*pointi];
1950                 }
1951                 if (abs(t - *told) < ttol) {
1952                         t = *told;
1953                         /*     update output part */
1954                 }
1955                 if (*told > t) {
1956                         /*     !  scheduling problem */
1957                         *ierr = 1;
1958                         freeallx;
1959                         return;
1960                 }
1961                 if (*told != t) {
1962                         if (xptr[nblk+1] == 1) {
1963                                 /*     .     no continuous state */
1964                                 if (*told + deltat + ttol > t) {
1965                                         *told = t;
1966                                 } else {
1967                                         *told += deltat;
1968                                 }
1969                                 /*     .     update outputs of 'c' type blocks with no continuous state */
1970                                 if (*told >= *tf) {
1971                                         /*     .     we are at the end, update continuous part before leaving */
1972                                         cdoit(told);
1973                                         freeallx;
1974                                         return;
1975                                 }
1976                         } else {
1977                                 rhotmp = *tf + ttol;
1978                                 if (*pointi!=0){
1979                                         kpo = *pointi;
1980 L20:
1981                                         if (critev[kpo] == 1) {
1982                                                 rhotmp = tevts[kpo];
1983                                                 goto L30;
1984                                         }
1985                                         kpo = evtspt[kpo];
1986                                         if (kpo != 0) {
1987                                                 goto L20;
1988                                         }
1989 L30:
1990                                         if (rhotmp < tstop) {
1991                                                 hot = 0;/* Do cold-restart the solver:if the new TSTOP isn't beyong the previous one*/ 
1992                                         }
1993                                 }
1994                                 tstop = rhotmp;
1995                                 t = min(*told + deltat,min(t,*tf + ttol));
1996
1997                                 if (hot == 0){ /* CIC calculation when hot==0 */
1998
1999                                         /* Setting the stop time*/
2000                                         flag = IDASetStopTime(ida_mem, (realtype)tstop);
2001                                         if (check_flag(&flag, "IDASetStopTime", 1)) {
2002                                                 *ierr=200+(-flag);
2003                                                 freeallx;
2004                                                 return;
2005                                         }
2006
2007                                         if (ng>0&&nmod>0){
2008                                                 phase=1;
2009                                                 zdoit(told, x, xd, g);
2010                                                 if (*ierr != 0) {
2011                                                         freeallx;
2012                                                         return;
2013                                                 }
2014                                         }
2015
2016                                         /*----------ID setting/checking------------*/
2017                                         N_VConst(ONE, IDx); /* Initialize id to 1's. */
2018                                         scicos_xproperty=NV_DATA_S(IDx);
2019                                         reinitdoit(told);   
2020                                         if(*ierr >0) {
2021                                                 freeallx;
2022                                                 return;
2023                                         }
2024                                         for(jj=0;jj<*neq;jj++) {
2025                                                 if (xprop[jj] ==  1) scicos_xproperty[jj] = ONE;
2026                                                 if (xprop[jj] == -1) scicos_xproperty[jj] = ZERO;
2027                                         }         
2028                                         /* CI=0.0;CJ=100.0; // for functions Get_Jacobian_ci and Get_Jacobian_cj
2029                                         Jacobians(*neq, (realtype) (*told), yy, yp,     bidon, (realtype) CJ, data, TJacque, tempv1, tempv2, tempv3);       
2030                                         for (jj=0;jj<*neq;jj++){
2031                                         Jacque_col=DENSE_COL(TJacque,jj);
2032                                         CI=ZERO;
2033                                         for (kk=0;kk<*neq;kk++){
2034                                         if ((Jacque_col[kk]-Jacque_col[kk]!=0)) {
2035                                         CI=-ONE;
2036                                         break;
2037                                         }else{
2038                                         if (Jacque_col[kk]!=0){
2039                                         CI=ONE;
2040                                         break;
2041                                         }
2042                                         }
2043                                         }           
2044                                         if (CI>=ZERO){  scicos_xproperty[jj]=CI;}else{fprintf(stderr,"\nWarinng! Xproperties are not match for i=%d!",jj);}
2045                                         } */
2046                                         /* printf("\n"); for(jj=0;jj<*neq;jj++) { printf("x%d=%g ",jj,scicos_xproperty[jj]); }*/
2047
2048                                         flag=IDASetId(ida_mem,IDx);
2049                                         if (check_flag(&flag, "IDASetId", 1)) {
2050                                                 *ierr=200+(-flag); 
2051                                                 freeallx
2052                                                         return;
2053                                         }
2054                                         // CI=1.0;  // for function Get_Jacobian_ci
2055                                         /*--------------------------------------------*/
2056                                         // maxnj=100; /* setting the maximum number of Jacobian evaluation during a Newton step */
2057                                         // flag=IDASetMaxNumJacsIC(ida_mem, maxnj);
2058                                         // if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
2059                                         //   *ierr=200+(-flag);
2060                                         //   freeallx;
2061                                         //   return;
2062                                         // };
2063                                         // flag=IDASetLineSearchOffIC(ida_mem,FALSE);  /* (def=false)  */
2064                                         // if (check_flag(&flag, "IDASetLineSearchOffIC", 1)) { 
2065                                         //   *ierr=200+(-flag);
2066                                         //   freeallx;
2067                                         //   return;
2068                                         // };
2069                                         // flag=IDASetMaxNumItersIC(ida_mem, 10);/* (def=10) setting the maximum number of Newton iterations in any one attemp to solve CIC */
2070                                         // if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) { 
2071                                         //   *ierr=200+(-flag);
2072                                         //   freeallx;
2073                                         //   return;
2074                                         // };
2075
2076                                         N_iters=4+nmod*4;
2077                                         for(j=0;j<=N_iters;j++){/* counter to reevaluate the
2078                                                                                         modes in  mode->CIC->mode->CIC-> loop
2079                                                                                         do it once in the absence of mode (nmod=0) */
2080                                                 /* updating the modes through Flag==9, Phase==1 */
2081
2082                                                 /* Serge Steer 29/06/2009 */
2083                                                 while (ismenu()) //** if the user has done something, do the actions
2084                                                 {
2085                                                         int ierr2=0;
2086                                                         SeqSync = GetCommand(CommandToUnstack); //** get at the action
2087                                                         CommandLength = (int)strlen(CommandToUnstack);
2088                                                         syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2089                                                 }
2090                                                 if (C2F(coshlt).halt != 0) {
2091                                                         C2F(coshlt).halt = 0;
2092                                                         freeallx;
2093                                                         return;
2094                                                 }
2095
2096                                                 /* yy->PH */
2097                                                 flag =IDAReInit(ida_mem, simblkdaskr, (realtype)(*told), yy, yp, IDA_SS, reltol, &abstol);
2098                                                 if (check_flag(&flag, "CVodeReInit", 1)) {
2099                                                         *ierr=200+(-flag);
2100                                                         freeallx;
2101                                                         return;
2102                                                 }
2103
2104                                                 phase=2; /* IDACalcIC: PHI-> yy0: if (ok) yy0_cic-> PHI*/
2105                                                 copy_IDA_mem->ida_kk=1;
2106
2107                                                 // the initial conditons y0 and yp0 do not satisfy the DAE
2108                                                 flagr=IDACalcIC(ida_mem, IDA_YA_YDP_INIT, (realtype)(t));
2109
2110                                                 phase=1;
2111                                                 flag = IDAGetConsistentIC(ida_mem, yy, yp); /* PHI->YY */
2112
2113                                                 if (*ierr > 5) {  /* *ierr>5 => singularity in block */
2114                                                         freeallx;
2115                                                         return;
2116                                                 }
2117
2118                                                 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2119                                                 {
2120                                                         if (flagr>=0) {
2121                                                                 sciprint(_("**** SUNDIALS.IDA successfully initialized *****\n") );
2122                                                         }
2123                                                         else{
2124                                                                 sciprint(_("**** SUNDIALS.IDA failed to initialize ->try again *****\n") );
2125                                                         }
2126                                                 }
2127                                                 /*-------------------------------------*/
2128                                                 /* saving the previous modes*/
2129                                                 for (jj = 0; jj < nmod; ++jj) {
2130                                                         Mode_save[jj] = mod[jj];
2131                                                 }
2132                                                 if (ng>0&&nmod>0){
2133                                                         phase=1;
2134                                                         zdoit(told, x, xd, g);
2135                                                         if (*ierr != 0) {
2136                                                                 freeallx;
2137                                                                 return; 
2138                                                         }
2139                                                 }
2140                                                 /*------------------------------------*/
2141                                                 Mode_change=0;
2142                                                 for (jj = 0; jj < nmod; ++jj) {
2143                                                         if(Mode_save[jj] != mod[jj])
2144                                                         {Mode_change=1;
2145                                                         break;
2146                                                         }
2147                                                 }
2148                                                 if (Mode_change==0) {
2149                                                         if (flagr >=0 ){
2150                                                                 break;  /*   if (flagr>=0) break;  else{ *ierr=200+(-flagr); freeallx;  return; }*/
2151                                                         }else if(j>=(int)( N_iters/2)) {
2152                                                                 /* IDASetMaxNumStepsIC(mem,10); */     /* maxnh (def=5) */
2153                                                                 IDASetMaxNumJacsIC(ida_mem,10);       /* maxnj 100 (def=4)*/
2154                                                                 /* IDASetMaxNumItersIC(mem,100000); */ /* maxnit in IDANewtonIC (def=10) */
2155                                                                 IDASetLineSearchOffIC(ida_mem,TRUE);  /* (def=false)  */
2156                                                                 /* IDASetNonlinConvCoefIC(mem,1.01);*/ /* (def=0.01-0.33*/
2157                                                                 flag=IDASetMaxNumItersIC(ida_mem, 1000);
2158                                                                 if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
2159                                                                         *ierr=200+(-flag); 
2160                                                                         freeallx; 
2161                                                                         return;
2162                                                                 };
2163                                                         }
2164                                                 }
2165                                         }/* mode-CIC  counter*/
2166                                         if(Mode_change==1){ 
2167                                                 /* In tghis case, we try again by relaxing all modes and calling IDA_calc again 
2168                                                 /Masoud */
2169                                                 phase=1;
2170                                                 copy_IDA_mem->ida_kk=1;
2171                                                 flagr=IDACalcIC(ida_mem, IDA_YA_YDP_INIT, (realtype)(t));
2172                                                 phase=1;
2173                                                 flag = IDAGetConsistentIC(ida_mem, yy, yp); /* PHI->YY */
2174                                                 if ((flagr<0)||(*ierr>5)) {  /* *ierr>5 => singularity in block */
2175                                                         *ierr = 23;
2176                                                         freeallx;
2177                                                         return;
2178                                                 }
2179                                         }
2180                                         /*-----If flagr<0 the initialization solver has not converged-----*/
2181                                         if (flagr<0) {
2182                                                 *ierr = 237;
2183                                                 freeallx;
2184                                                 return;
2185                                         }
2186
2187                                 } /* CIC calculation when hot==0 */
2188
2189                                 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2190                                 {
2191                                         sciprint(_("****daskr from: %f to %f hot= %d  \n"), *told,t,hot);
2192                                 }
2193
2194                                 /*--discrete zero crossings----dzero--------------------*/
2195                                 /*--check for Dzeros after Mode settings or ddoit()----*/
2196                                 Discrete_Jump=0;
2197                                 if (ng>0 && hot==0){
2198                                         zdoit(told, x, xd, g);
2199                                         if (*ierr != 0) {freeallx;return;  }
2200                                         for (jj = 0; jj < ng; ++jj) {
2201                                                 if((g[jj]>=0.0)&&( jroot[jj]==-5)) {Discrete_Jump=1;jroot[jj]=1;}
2202                                                 else if((g[jj]<0.0)&&( jroot[jj]==5)) {Discrete_Jump=1;jroot[jj]=-1;}
2203                                                 else jroot[jj]=0;
2204                                         }
2205                                 }
2206
2207                                 /*--discrete zero crossings----dzero--------------------*/
2208                                 if (Discrete_Jump==0){/* if there was a dzero, its event should be activated*/
2209                                         phase=2;
2210                                         flagr = IDASolve(ida_mem, t, told, yy, yp, IDA_NORMAL_TSTOP);
2211                                         phase=1;
2212                                         if (*ierr != 0) {freeallx;return;}
2213                                 }else{
2214                                         flagr = IDA_ROOT_RETURN; /* in order to handle discrete jumps */
2215                                 }
2216                                 if (flagr>=0){
2217                                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2218                                                 sciprint(_("****SUNDIALS.Ida reached: %f\n"),*told);
2219                                         hot = 1;
2220                                         cnt = 0;
2221                                 } else if ( flagr==IDA_TOO_MUCH_WORK ||  flagr == IDA_CONV_FAIL || flagr==IDA_ERR_FAIL) {  
2222                                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2223                                                 sciprint(_("**** SUNDIALS.Ida: too much work at time=%g (stiff region, change RTOL and ATOL)\n"),*told);
2224                                         hot = 0;
2225                                         cnt++;
2226                                         if (cnt>5) {
2227                                                 *ierr=200+(-flagr);
2228                                                 freeallx;
2229                                                 return;
2230                                         }
2231                                 }else{
2232                                         if (flagr<0) *ierr=200+(-flagr);    /* raising errors due to internal errors, other wise erros due to flagr*/
2233                                         freeallx;
2234                                         return;
2235                                 }
2236
2237                                 /*     update outputs of 'c' type  blocks if we are at the end*/
2238                                 if (*told >= *tf) {
2239                                         cdoit(told);
2240                                         freeallx;
2241                                         return;
2242                                 }
2243
2244                                 if (flagr == IDA_ZERO_DETACH_RETURN){hot=0;}; /* new feature of sundials, detects unmasking */
2245                                 if (flagr == IDA_ROOT_RETURN) {
2246                                         /*     .        at a least one root has been found */
2247                                         hot = 0;
2248                                         if (Discrete_Jump==0){
2249                                                 flagr = IDAGetRootInfo(ida_mem, jroot);
2250                                                 if (check_flag(&flagr, "IDAGetRootInfo", 1)) {
2251                                                         *ierr=200+(-flagr);    
2252                                                         freeallx;
2253                                                         return;
2254                                                 }
2255                                         }
2256
2257                                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2258                                         {
2259                                                 sciprint(_("root found at t=: %f\n"),*told);
2260                                         }
2261                                         /*     .        update outputs affecting ztyp blocks  ONLY FOR OLD BLOCKS*/
2262                                         zdoit(told, x, xd, g);
2263                                         if (*ierr != 0) {
2264                                                 freeallx;
2265                                                 return;
2266                                         }
2267                                         for (jj = 0; jj < ng; ++jj) {
2268                                                 C2F(curblk).kfun = zcros[jj];
2269                                                 if (C2F(curblk).kfun == -1) {
2270                                                         break; 
2271                                                 }
2272                                                 kev = 0;
2273                                                 for (j = zcptr[C2F(curblk).kfun]-1 ; 
2274                                                         j <zcptr[C2F(curblk).kfun+1]-1 ; ++j) {
2275                                                                 if(jroot[j]!=0){
2276                                                                         kev=1;
2277                                                                         break;
2278                                                                 }
2279                                                 }
2280                                                 if (kev != 0) {
2281                                                         Blocks[C2F(curblk).kfun-1].jroot=&jroot[zcptr[C2F(curblk).kfun]-1];
2282                                                         if (funtyp[C2F(curblk).kfun] > 0) {
2283                                                                 if (Blocks[C2F(curblk).kfun-1].nevout > 0) {
2284                                                                         flag__ = 3;
2285                                                                         if(Blocks[C2F(curblk).kfun-1].nx > 0) {
2286                                                                                 Blocks[C2F(curblk).kfun-1].x  = &x[xptr[C2F(curblk).kfun]-1];
2287                                                                                 Blocks[C2F(curblk).kfun-1].xd = &xd[xptr[C2F(curblk).kfun]-1];
2288                                                                         }
2289                                                                         /*     call corresponding block to determine output event (kev) */
2290                                                                         Blocks[C2F(curblk).kfun-1].nevprt = -kev;
2291                                                                         /*callf(told, xd, x, x,g,&flag__);*/
2292                                                                         callf(told, &Blocks[C2F(curblk).kfun-1], &flag__);
2293                                                                         if (flag__ < 0) {
2294                                                                                 *ierr = 5 - flag__;
2295                                                                                 freeallx;
2296                                                                                 return;
2297                                                                         }
2298                                                                         /*     update event agenda */
2299                                                                         for (k = 0; k < Blocks[C2F(curblk).kfun-1].nevout; ++k) {
2300                                                                                 if (Blocks[C2F(curblk).kfun-1].evout[k] >= 0) {
2301                                                                                         i3 = k + clkptr[C2F(curblk).kfun] ;
2302                                                                                         addevs(Blocks[C2F(curblk).kfun-1].evout[k]+(*told), &i3, &ierr1);
2303                                                                                         if (ierr1 != 0) {
2304                                                                                                 /*     .                       nevts too small */
2305                                                                                                 *ierr = 3;
2306                                                                                                 freeallx;
2307                                                                                                 return;
2308                                                                                         }
2309                                                                                 }
2310                                                                         }
2311                                                                 }
2312                                                                 /* update state */
2313                                                                 if ((Blocks[C2F(curblk).kfun-1].nx > 0) || (*Blocks[C2F(curblk).kfun-1].work != NULL) ) {
2314                                                                         /* call corresponding block to update state */
2315                                                                         flag__ = 2;
2316                                                                         if(Blocks[C2F(curblk).kfun-1].nx > 0) {
2317                                                                                 Blocks[C2F(curblk).kfun-1].x  = &x[xptr[C2F(curblk).kfun]-1];
2318                                                                                 Blocks[C2F(curblk).kfun-1].xd = &xd[xptr[C2F(curblk).kfun]-1];
2319                                                                         }
2320                                                                         Blocks[C2F(curblk).kfun-1].nevprt = -kev;
2321
2322                                                                         Blocks[C2F(curblk).kfun-1].xprop = &xprop[-1+xptr[C2F(curblk).kfun]];
2323                                                                         /*callf(told, xd, x, x,g,&flag__);*/
2324                                                                         callf(told, &Blocks[C2F(curblk).kfun-1], &flag__);
2325
2326                                                                         if (flag__ < 0) {
2327                                                                                 *ierr = 5 - flag__;
2328                                                                                 freeallx;
2329                                                                                 return;
2330                                                                         }
2331                                                                         for(j=0;j<*neq;j++) { /* Adjust xprop for IDx */
2332                                                                                 if (xprop[j] ==  1) scicos_xproperty[j] = ONE;
2333                                                                                 if (xprop[j] == -1) scicos_xproperty[j] = ZERO;
2334                                                                         }
2335                                                                 }
2336                                                         }
2337                                                 }
2338                                         }
2339                                 }
2340                                 /* Serge Steer 29/06/2009 */
2341                                 while (ismenu()) //** if the user has done something, do the actions
2342                                 {
2343                                         int ierr2=0;
2344                                         SeqSync = GetCommand(CommandToUnstack); //** get at the action
2345                                         CommandLength = (int)strlen(CommandToUnstack);
2346                                         syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2347                                 }
2348
2349                                 if (C2F(coshlt).halt != 0) {
2350                                         C2F(coshlt).halt = 0;
2351                                         freeallx;
2352                                         return;
2353                                 }
2354                                 /* if(*pointi!=0){
2355                                 t=tevts[*pointi];
2356                                 if(*told<t-ttol){
2357                                 cdoit(told);
2358                                 goto L15;
2359                                 }
2360                                 }else{
2361                                 if(*told<*tf){
2362                                 cdoit(told);
2363                                 goto L15;
2364                                 }
2365                                 }*/
2366
2367                                 /*--discrete zero crossings----dzero--------------------*/
2368                                 if (ng>0) { /* storing ZC signs just after a ddaskr call*/
2369                                         zdoit(told, x, xd, g);
2370                                         if (*ierr != 0) {
2371                                                 freeallx;
2372                                                 return;
2373                                         }
2374                                         for (jj = 0; jj < ng; ++jj) {
2375                                                 if(g[jj]>=0) {
2376                                                         jroot[jj]=5;
2377                                                 }
2378                                                 else {
2379                                                         jroot[jj]=-5;
2380                                                 }
2381                                         }
2382                                 }
2383                                 /*--discrete zero crossings----dzero--------------------*/
2384                         }
2385                         C2F(realtime)(told);
2386                 } else {
2387                         /*     .  t==told */
2388                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3)) {
2389                                 sciprint(_("Event: %d activated at t=%f\n"),*pointi,*told);
2390                         }
2391
2392                         ddoit(told);
2393                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3)) {
2394                                 sciprint(_("End of activation"));
2395                         }
2396                         if (*ierr != 0) {
2397                                 freeallx;
2398                                 return;
2399                         }
2400                 }
2401                 /*     end of main loop on time */
2402         }
2403         freeallx;
2404 } /* cossimdaskr_ */
2405 /*--------------------------------------------------------------------------*/
2406 /* Subroutine cosend */
2407 static void cosend(double *told)
2408 {
2409         /* Local variables */
2410         static scicos_flag flag__ = 0;
2411
2412         static int kfune = 0;
2413
2414         /* Function Body */
2415         *ierr = 0;
2416         /*     loop on blocks */
2417         for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun) {
2418                 flag__ = 5;
2419                 Blocks[C2F(curblk).kfun-1].nevprt = 0;
2420                 if (funtyp[C2F(curblk).kfun] >= 0) {
2421                         if(Blocks[C2F(curblk).kfun-1].nx > 0) {
2422                                 Blocks[C2F(curblk).kfun-1].x  = &x[xptr[C2F(curblk).kfun]-1];
2423                                 Blocks[C2F(curblk).kfun-1].xd = &xd[xptr[C2F(curblk).kfun]-1];
2424                         }
2425                         /*callf(told, xd, x, x,x,&flag__);*/
2426                         callf(told, &Blocks[C2F(curblk).kfun-1], &flag__);
2427                         if (flag__ < 0 && *ierr == 0) {
2428                                 *ierr = 5 - flag__;
2429                                 kfune = C2F(curblk).kfun;
2430                         }
2431                 }
2432         }
2433         if (*ierr != 0) {
2434                 C2F(curblk).kfun = kfune;
2435                 return;
2436         }
2437 } /* cosend_ */
2438 /*--------------------------------------------------------------------------*/
2439 /* callf */
2440 void callf(double *t, scicos_block *block, scicos_flag *flag)
2441 {
2442         double* args[SZ_SIZE];
2443         int sz[SZ_SIZE];
2444         double intabl[TB_SIZE];
2445         double outabl[TB_SIZE];
2446
2447         int ii = 0,in = 0,out = 0,ki = 0,ko = 0,no = 0,ni = 0,k = 0,j = 0;
2448         int szi = 0,flagi = 0;
2449         double *ptr_d=NULL;
2450
2451         /* function pointers type def */
2452         voidf loc ;
2453         ScicosF0 loc0;
2454         ScicosF loc1;
2455         /*  ScicosFm1 loc3;*/
2456         ScicosF2 loc2;
2457         ScicosF2z loc2z;
2458         ScicosFi loci1;
2459         ScicosFi2 loci2;
2460         ScicosFi2z loci2z;
2461         ScicosF4 loc4;
2462
2463         int solver = C2F(cmsolver).solver;
2464         int cosd   = C2F(cosdebug).cosd;
2465         /*int kf     = C2F(curblk).kfun;*/
2466         scicos_time     = *t;
2467         block_error     = (int*) flag;
2468
2469         /* debug block is never called */
2470         /*if (kf==(debug_block+1)) return;*/
2471         if (block->type==99) return;
2472
2473         /* flag 7 implicit initialization */
2474         flagi = (int) *flag;
2475         /* change flag to zero if flagi==7 for explicit block */
2476         if(flagi==7 && block->type<10000) {
2477                 *flag=0;
2478         }
2479
2480         /* display information for debugging mode */
2481         if (cosd > 1) {
2482                 if (cosd != 3) {
2483                         sciprint(_("block %d is called "),C2F(curblk).kfun);
2484                         sciprint(_("with flag %d "),*flag);
2485                         sciprint(_("at time %f \n"),*t);
2486                 }
2487                 if(debug_block>-1) {
2488                         if (cosd != 3) sciprint(_("Entering the block \n"));
2489                         call_debug_scicos(block,flag,flagi,debug_block);
2490                         if (*flag<0) return;  /* error in debug block */
2491                 }
2492         }
2493
2494         C2F(scsptr).ptr = block->scsptr;
2495
2496         /* get pointer of the function */
2497         loc=block->funpt;
2498
2499         /* continuous state */
2500         if(solver==100 && block->type<10000 && *flag==0) {
2501                 ptr_d = block->xd;
2502                 block->xd  = block->res;
2503         }
2504
2505         /* switch loop */
2506         //sciprint("callf type=%d flag=%d\n",block->type,flagi);
2507         switch (block->type) {
2508                 /*******************/
2509                 /* function type 0 */
2510                 /*******************/
2511   case 0 :
2512           {  /* This is for compatibility */
2513                   /* jroot is returned in g for old type */
2514                   if(block->nevprt<0) {
2515                           for (j =0;j<block->ng;++j) {
2516                                   block->g[j] = (double)block->jroot[j];
2517                           }
2518                   }
2519
2520                   /* concatenated entries and concatened outputs */
2521                   /* catenate inputs if necessary */
2522                   ni=0;
2523                   if (block->nin>1) {
2524                           ki=0;
2525                           for (in=0;in<block->nin;in++) {
2526                                   szi=block->insz[in]*block->insz[in+block->nin];
2527                                   for (ii=0;ii<szi;ii++) {
2528                                           intabl[ki++]= *((double *)(block->inptr[in]) + ii);
2529                                   }
2530                                   ni=ni+szi;
2531                           }
2532                           args[0]=&(intabl[0]);
2533                   }
2534                   else {
2535                           if (block->nin==0) {
2536                                   args[0]=NULL;
2537                           }
2538                           else {
2539                                   args[0]= (double *)(block->inptr[0]);
2540                                   ni=block->insz[0]*block->insz[1];
2541                           }
2542                   }
2543
2544                   /* catenate outputs if necessary */
2545                   no=0;
2546                   if (block->nout>1) {
2547                           ko=0;
2548                           for (out=0;out<block->nout;out++) {
2549                                   szi=block->outsz[out]*block->outsz[out+block->nout];
2550                                   for (ii=0;ii<szi;ii++) {
2551                                           outabl[ko++]= *((double *)(block->outptr[out]) + ii);
2552                                   }
2553                                   no=no+szi;
2554                           }
2555                           args[1]=&(outabl[0]);
2556                   }
2557                   else {
2558                           if (block->nout==0) {
2559                                   args[1]=NULL;
2560                           }
2561                           else {
2562                                   args[1]= (double *)(block->outptr[0]);
2563                                   no=block->outsz[0]*block->outsz[1];
2564                           }
2565                   }
2566
2567                   loc0 = (ScicosF0) loc;
2568
2569                   (*loc0)(flag,&block->nevprt,t,block->xd,block->x,&block->nx,
2570                           block->z,&block->nz,
2571                           block->evout,&block->nevout,block->rpar,&block->nrpar,
2572                           block->ipar,&block->nipar,(double *)args[0],&ni,
2573                           (double *)args[1],&no);
2574
2575                   /* split output vector on each port if necessary */
2576                   if (block->nout>1) {
2577                           ko=0;
2578                           for (out=0;out<block->nout;out++) {
2579                                   szi=block->outsz[out]*block->outsz[out+block->nout];
2580                                   for (ii=0;ii<szi;ii++) {
2581                                           *((double *)(block->outptr[out]) + ii) = outabl[ko++];
2582                                   }
2583                           }
2584                   }
2585
2586                   /* adjust values of output register */
2587                   for(in=0;in<block->nevout;++in) {
2588                           block->evout[in]=block->evout[in]-*t;
2589                   }
2590
2591                   break;
2592           }
2593
2594           /*******************/
2595           /* function type 1 */
2596           /*******************/
2597   case 1 :
2598           { /* This is for compatibility */
2599                   /* jroot is returned in g for old type */
2600                   if(block->nevprt<0) {
2601                           for (j =0;j<block->ng;++j) {
2602                                   block->g[j] = (double)block->jroot[j];
2603                           }
2604                   }
2605
2606                   /* one entry for each input or output */
2607                   for (in = 0 ; in < block->nin ; in++) {
2608                           args[in]=block->inptr[in];
2609                           sz[in]=block->insz[in];
2610                   }
2611                   for (out=0;out<block->nout;out++) {
2612                           args[in+out]=block->outptr[out];
2613                           sz[in+out]=block->outsz[out];
2614                   }
2615                   /* with zero crossing */
2616                   if(block->ztyp>0) {
2617                           args[block->nin+block->nout]=block->g;
2618                           sz[block->nin+block->nout]=block->ng;
2619                   }
2620
2621                   loc1 = (ScicosF) loc;
2622
2623                   (*loc1)(flag,&block->nevprt,t,block->xd,block->x,&block->nx,
2624                           block->z,&block->nz,
2625                           block->evout,&block->nevout,block->rpar,&block->nrpar,
2626                           block->ipar,&block->nipar,
2627                           (double *)args[0],&sz[0],
2628                           (double *)args[1],&sz[1],(double *)args[2],&sz[2],
2629                           (double *)args[3],&sz[3],(double *)args[4],&sz[4],
2630                           (double *)args[5],&sz[5],(double *)args[6],&sz[6],
2631                           (double *)args[7],&sz[7],(double *)args[8],&sz[8],
2632                           (double *)args[9],&sz[9],(double *)args[10],&sz[10],
2633                           (double *)args[11],&sz[11],(double *)args[12],&sz[12],
2634                           (double *)args[13],&sz[13],(double *)args[14],&sz[14],
2635                           (double *)args[15],&sz[15],(double *)args[16],&sz[16],
2636                           (double *)args[17],&sz[17]);
2637
2638                   /* adjust values of output register */
2639                   for(in=0;in<block->nevout;++in) {
2640                           block->evout[in]=block->evout[in]-*t;
2641                   }
2642
2643                   break;
2644           }
2645
2646           /*******************/
2647           /* function type 2 */
2648           /*******************/
2649   case 2 :
2650           { /* This is for compatibility */
2651                   /* jroot is returned in g for old type */
2652                   if(block->nevprt<0) {
2653                           for (j =0;j<block->ng;++j) {
2654                                   block->g[j] = (double)block->jroot[j];
2655                           }
2656                   }
2657
2658                   /* no zero crossing */
2659                   if (block->ztyp==0) {
2660                           loc2 = (ScicosF2) loc;
2661                           (*loc2)(flag,&block->nevprt,t,block->xd,block->x,&block->nx,
2662                                   block->z,&block->nz,
2663                                   block->evout,&block->nevout,block->rpar,&block->nrpar,
2664                                   block->ipar,&block->nipar,(double **)block->inptr,
2665                                   block->insz,&block->nin,
2666                                   (double **)block->outptr,block->outsz,&block->nout);
2667                   }
2668                   /* with zero crossing */
2669                   else {
2670                           loc2z = (ScicosF2z) loc;
2671                           (*loc2z)(flag,&block->nevprt,t,block->xd,block->x,&block->nx,
2672                                   block->z,&block->nz,
2673                                   block->evout,&block->nevout,block->rpar,&block->nrpar,
2674                                   block->ipar,&block->nipar,(double **)block->inptr,
2675                                   block->insz,&block->nin,
2676                                   (double **)block->outptr,block->outsz,&block->nout,
2677                                   block->g,&block->ng);
2678                   }
2679
2680                   /* adjust values of output register */
2681                   for(in=0;in<block->nevout;++in) {
2682                           block->evout[in]=block->evout[in]-*t;
2683                   }
2684
2685                   break;
2686           }
2687
2688           /*******************/
2689           /* function type 4 */
2690           /*******************/
2691   case 4 :
2692           { /* get pointer of the function type 4*/
2693                   loc4 = (ScicosF4) loc;
2694
2695                   (*loc4)(block,*flag);
2696
2697                   break;
2698           }
2699
2700           /***********************/
2701           /* function type 10001 */
2702           /***********************/
2703   case 10001 :
2704           { /* This is for compatibility */
2705                   /* jroot is returned in g for old type */
2706                   if(block->nevprt<0) {
2707                           for (j =0;j<block->ng;++j) {
2708                                   block->g[j] = (double)block->jroot[j];
2709                           }
2710                   }
2711
2712                   /* implicit block one entry for each input or output */
2713                   for (in = 0 ; in < block->nin ; in++) {
2714                           args[in]=block->inptr[in];
2715                           sz[in]=block->insz[in];
2716                   }
2717                   for (out=0;out<block->nout;out++) {
2718                           args[in+out]=block->outptr[out];
2719                           sz[in+out]=block->outsz[out];
2720                   }
2721                   /* with zero crossing */
2722                   if(block->ztyp>0) {
2723                           args[block->nin+block->nout]=block->g;
2724                           sz[block->nin+block->nout]=block->ng;
2725                   }
2726
2727                   loci1 = (ScicosFi) loc;
2728                   (*loci1)(flag,&block->nevprt,t,block->res,block->xd,block->x,
2729                           &block->nx,block->z,&block->nz,
2730                           block->evout,&block->nevout,block->rpar,&block->nrpar,
2731                           block->ipar,&block->nipar,
2732                           (double *)args[0],&sz[0],
2733                           (double *)args[1],&sz[1],(double *)args[2],&sz[2],
2734                           (double *)args[3],&sz[3],(double *)args[4],&sz[4],
2735                           (double *)args[5],&sz[5],(double *)args[6],&sz[6],
2736                           (double *)args[7],&sz[7],(double *)args[8],&sz[8],
2737                           (double *)args[9],&sz[9],(double *)args[10],&sz[10],
2738                           (double *)args[11],&sz[11],(double *)args[12],&sz[12],
2739                           (double *)args[13],&sz[13],(double *)args[14],&sz[14],
2740                           (double *)args[15],&sz[15],(double *)args[16],&sz[16],
2741                           (double *)args[17],&sz[17]);
2742
2743                   /* adjust values of output register */
2744                   for(in=0;in<block->nevout;++in) {
2745                           block->evout[in]=block->evout[in]-*t;
2746                   }
2747
2748                   break;
2749           }
2750
2751           /***********************/
2752           /* function type 10002 */
2753           /***********************/
2754   case 10002 :
2755           { /* This is for compatibility */
2756                   /* jroot is returned in g for old type */
2757                   if(block->nevprt<0) {
2758                           for (j =0;j<block->ng;++j) {
2759                                   block->g[j] = (double)block->jroot[j];
2760                           }
2761                   }
2762
2763                   /* implicit block, inputs and outputs given by a table of pointers */
2764                   /* no zero crossing */
2765                   if(block->ztyp==0) {
2766                           loci2 = (ScicosFi2) loc;
2767                           (*loci2)(flag,&block->nevprt,t,block->res,
2768                                   block->xd,block->x,&block->nx,
2769                                   block->z,&block->nz,
2770                                   block->evout,&block->nevout,block->rpar,&block->nrpar,
2771                                   block->ipar,&block->nipar,(double **)block->inptr,
2772                                   block->insz,&block->nin,
2773                                   (double **)block->outptr,block->outsz,&block->nout);
2774                   }
2775                   /* with zero crossing */
2776                   else {
2777                           loci2z = (ScicosFi2z) loc;
2778                           (*loci2z)(flag,&block->nevprt,t,block->res,
2779                                   block->xd,block->x,&block->nx,
2780                                   block->z,&block->nz,
2781                                   block->evout,&block->nevout,block->rpar,&block->nrpar,
2782                                   block->ipar,&block->nipar,
2783                                   (double **)block->inptr,block->insz,&block->nin,
2784                                   (double **)block->outptr,block->outsz,&block->nout,
2785                                   block->g,&block->ng);
2786                   }
2787
2788                   /* adjust values of output register */
2789                   for(in=0;in<block->nevout;++in) {
2790                           block->evout[in]=block->evout[in]-*t;
2791                   }
2792
2793                   break;
2794           }
2795
2796           /***********************/
2797           /* function type 10004 */
2798           /***********************/
2799   case 10004 :
2800           { /* get pointer of the function type 4*/
2801                   loc4 = (ScicosF4) loc;
2802
2803                   (*loc4)(block,*flag);
2804
2805                   break;
2806           }
2807
2808           /***********/
2809           /* default */
2810           /***********/
2811   default :
2812           {
2813                   sciprint(_("Undefined Function type\n"));
2814                   *flag=-1000;
2815                   return; /* exit */
2816           }
2817         }
2818         // sciprint("callf end  flag=%d\n",*flag);
2819         /* Implicit Solver & explicit block & flag==0 */
2820         /* adjust continuous state vector after call */
2821         if(solver==100 && block->type<10000 && *flag==0) {
2822                 block->xd  = ptr_d;
2823                 if(flagi!=7) {
2824                         for (k=0;k<block->nx;k++) {
2825                                 block->res[k]=block->res[k]-block->xd[k];
2826                         }
2827                 }
2828                 else {
2829                         for (k=0;k<block->nx;k++) {
2830                                 block->xd[k]=block->res[k];
2831                         }
2832                 }
2833         }
2834
2835         /* debug block */
2836         if (cosd > 1) {
2837                 if(debug_block>-1) {
2838                         if (*flag<0) return;  /* error in block */
2839                         if (cosd != 3) sciprint(_("Leaving block %d \n"),C2F(curblk).kfun);
2840                         call_debug_scicos(block,flag,flagi,debug_block);
2841                         /*call_debug_scicos(flag,kf,flagi,debug_block);*/
2842                 }
2843         }
2844 } /* callf */
2845 /*--------------------------------------------------------------------------*/
2846 /* call_debug_scicos */
2847 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk)
2848 {
2849         voidf loc ;
2850         int solver=C2F(cmsolver).solver, k = 0;
2851         ScicosF4 loc4;
2852         double *ptr_d=NULL;
2853
2854         C2F(cosdebugcounter).counter = C2F(cosdebugcounter).counter+1;
2855         C2F(scsptr).ptr = Blocks[deb_blk].scsptr;
2856
2857         loc  = Blocks[deb_blk].funpt; /* GLOBAL */
2858         loc4 = (ScicosF4) loc;
2859
2860         /* continuous state */
2861         if(solver==100 && block->type<10000 && *flag==0) {
2862                 ptr_d = block->xd;
2863                 block->xd  = block->res;
2864         }
2865
2866         (*loc4)(block,*flag);
2867
2868         /* Implicit Solver & explicit block & flag==0 */
2869         /* adjust continuous state vector after call */
2870         if(solver==100 && block->type<10000 && *flag==0) {
2871                 block->xd  = ptr_d;
2872                 if(flagi!=7) {
2873                         for (k=0;k<block->nx;k++) {
2874                                 block->res[k]=block->res[k]-block->xd[k];
2875                         }
2876                 }
2877                 else {
2878                         for (k=0;k<block->nx;k++) {
2879                                 block->xd[k]=block->res[k];
2880                         }
2881                 }
2882         }
2883
2884         if (*flag<0) sciprint(_("Error in the Debug block \n"));
2885 } /* call_debug_scicos */
2886 /*--------------------------------------------------------------------------*/
2887 /* simblk */
2888 static int simblk(realtype t,N_Vector yy,N_Vector yp, void *f_data)
2889 {
2890         double tx = 0., *x = NULL, *xd = NULL;
2891         int i = 0, nantest = 0;
2892
2893         tx= (double) t;
2894         x=  (double *) NV_DATA_S(yy);
2895         xd= (double *) NV_DATA_S(yp);
2896
2897         for(i=0;i<*neq;i++)   xd[i]=0;   /* à la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
2898         C2F(ierode).iero = 0;
2899         *ierr= 0;
2900         odoit(&tx, x, xd, xd);
2901         C2F(ierode).iero = *ierr;
2902
2903         if (*ierr==0) {
2904                 nantest=0;
2905                 for (i=0;i<*neq;i++) { /* NaN checking */
2906                         if ((xd[i]-xd[i]!=0)) {
2907                                 sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"),i);
2908                                 nantest=1;
2909                                 break;
2910                         }
2911                 }
2912                 if (nantest==1) return 349;/* recoverable error; */
2913         }
2914
2915         return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
2916
2917 } /* simblk */
2918 /*--------------------------------------------------------------------------*/
2919 /* grblk */
2920 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data)
2921 {
2922         double tx = 0., *x = NULL;
2923         int jj = 0, nantest = 0;
2924
2925         tx= (double) t;
2926         x=  (double *) NV_DATA_S(yy);
2927
2928         C2F(ierode).iero = 0;
2929         *ierr= 0;
2930
2931         zdoit(&tx, x, x, (double*) gout);
2932
2933         if (*ierr==0) {
2934                 nantest=0;
2935                 for (jj=0;jj<ng;jj++)
2936                         if (gout[jj]-gout[jj]!=0){
2937                                 sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"),jj);
2938                                 nantest=1;break;} /* NaN checking */
2939                         if (nantest==1) return 350;/* recoverable error; */
2940         }
2941         C2F(ierode).iero = *ierr;
2942
2943         return 0;
2944 } /* grblk */
2945 /*--------------------------------------------------------------------------*/
2946 /* simblkdaskr */
2947 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata)
2948 {
2949         double tx = 0.;
2950         double *xc = NULL, *xcdot = NULL, *residual = NULL;
2951         realtype alpha = 0.;
2952
2953         UserData data;
2954
2955         realtype hh = 0.;
2956         int qlast = 0;
2957         int jj = 0,flag = 0, nantest = 0;
2958
2959         data = (UserData) rdata; 
2960
2961         if(get_phase_simulation()==1) {
2962                 /* Just to update mode in a very special case, i.e., when initialization using modes fails.
2963                 in this case, we relax all modes and try again one more time.
2964                 */
2965                 zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), (double *)data->gwork);
2966         }
2967
2968         hh=ZERO;
2969         flag=IDAGetCurrentStep(data->ida_mem, &hh);
2970         if (flag<0) {  *ierr=200+(-flag); return (*ierr);};
2971
2972         qlast=0;
2973         flag=IDAGetCurrentOrder(data->ida_mem, &qlast);
2974         if (flag<0) {  *ierr=200+(-flag); return (*ierr);};
2975
2976         alpha=ZERO;
2977         for (jj=0;jj<qlast;jj++)
2978                 alpha=alpha -ONE/(jj+1);
2979         if (hh!=0) 
2980                 // CJ=-alpha/hh;  // For function Get_Jacobian_cj 
2981                 CJJ=-alpha/hh; 
2982         else {
2983                 *ierr= 217;return (*ierr);
2984         }
2985         xc=(double *)  NV_DATA_S(yy);
2986         xcdot= (double *) NV_DATA_S(yp);
2987         residual=(double *) NV_DATA_S(resval);
2988         tx=(double) tres;
2989
2990         C2F(dcopy)(neq, xcdot, &c__1, residual, &c__1);
2991         *ierr= 0;
2992         C2F(ierode).iero = 0;
2993         odoit(&tx, xc, xcdot, residual);
2994
2995         C2F(ierode).iero = *ierr;
2996
2997         if (*ierr==0) {
2998                 nantest=0;
2999                 for (jj=0;jj<*neq;jj++)
3000                         if (residual[jj]-residual[jj]!=0){/* NaN checking */
3001                                 //sciprint(_("\nWarning: The residual function #%d returns a NaN"),jj);
3002                                 nantest=1;
3003                                 break;
3004                         } 
3005                         if (nantest==1) return 257;/* recoverable error; */
3006         }
3007
3008         return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
3009 }/* simblkdaskr */
3010 /*--------------------------------------------------------------------------*/
3011 /* grblkdaskr */
3012 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data)
3013 {
3014         double tx = 0.;
3015         int jj = 0, nantest = 0;
3016
3017         tx=(double) t;
3018
3019         *ierr= 0;
3020         C2F(ierode).iero = 0;
3021         zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), (double *)gout);
3022         if (*ierr==0) {
3023                 nantest=0; /* NaN checking */
3024                 for (jj=0;jj<ng;jj++) {
3025                         if (gout[jj]-gout[jj]!=0) {
3026                                 sciprint(_("\nWarning: The zero-crossing function #%d returns a NaN"),jj);
3027                                 nantest=1;
3028                                 break;
3029                         }
3030                 }
3031                 if (nantest==1) {
3032                         return 258; /* recoverable error; */
3033                 }
3034         }
3035         C2F(ierode).iero = *ierr;
3036         return (*ierr);
3037 }/* grblkdaskr */
3038 /*--------------------------------------------------------------------------*/
3039 /* Subroutine addevs */
3040 static void addevs(double t, int *evtnb, int *ierr1)
3041 {
3042         static int i = 0, j = 0;
3043
3044         /* Function Body */
3045         *ierr1 = 0;
3046         if (evtspt[*evtnb] != -1) {
3047                 if ((evtspt[*evtnb] == 0) && (*pointi ==*evtnb)) {
3048                         tevts[*evtnb] = t;
3049                         return;
3050                 }else{
3051                         if (*pointi == *evtnb) {
3052                                 *pointi =evtspt[*evtnb]; /* remove from chain */
3053                         }else{
3054                                 i= *pointi;
3055                                 while (*evtnb != evtspt[i]){
3056                                         i=evtspt[i];
3057                                 }
3058                                 evtspt[i]=evtspt[*evtnb]; /* remove old evtnb from chain */
3059                                 if (TCritWarning==0){
3060                                         sciprint(_("\n Warning: an event is reprogrammed at t=%g by removing another"),t );
3061                                         sciprint(_("\n         (already programmed) event. There may be an error in"));
3062                                         sciprint(_("\n         your model. Please check your model\n"));
3063                                         TCritWarning=1;
3064                                 }
3065                                 do_cold_restart(); /* the erased event could be a critical
3066                                                                    event, so do_cold_restart is added to
3067                                                                    refresh the critical event table */
3068                         }
3069                         evtspt[*evtnb] = 0;
3070                         tevts[*evtnb] = t;
3071                 }
3072         }else {
3073                 evtspt[*evtnb] = 0;
3074                 tevts[*evtnb] = t;
3075         }
3076         if (*pointi == 0) {
3077                 *pointi = *evtnb;
3078                 return;
3079         }
3080         if (t < tevts[*pointi]) {
3081                 evtspt[*evtnb] = *pointi;
3082                 *pointi = *evtnb;
3083                 return;
3084         }
3085         i = *pointi;
3086
3087 L100:
3088         if (evtspt[i] == 0) {
3089                 evtspt[i] = *evtnb;
3090                 return;
3091         }
3092         if (t >= tevts[evtspt[i]]) {
3093                 j = evtspt[i];
3094                 if (evtspt[j] == 0) {
3095                         evtspt[j] = *evtnb;
3096                         return;
3097                 }
3098                 i = j;
3099                 goto L100;
3100         } else {
3101                 evtspt[*evtnb] = evtspt[i];
3102                 evtspt[i] = *evtnb;
3103         }
3104 } /* addevs */
3105 /*--------------------------------------------------------------------------*/
3106 /* Subroutine putevs */
3107 void putevs(double *t, int *evtnb, int *ierr1)
3108 {
3109         /* Function Body */
3110         *ierr1 = 0;
3111         if (evtspt[*evtnb] != -1) {
3112                 *ierr1 = 1;
3113                 return;
3114         } else {
3115                 evtspt[*evtnb] = 0;
3116                 tevts[*evtnb] = *t;
3117         }
3118         if (*pointi == 0) {
3119                 *pointi = *evtnb;
3120                 return;
3121         }
3122         evtspt[*evtnb] = *pointi;
3123         *pointi = *evtnb;
3124 } /* putevs */
3125 /*--------------------------------------------------------------------------*/
3126 /* Subroutine idoit */
3127 static void idoit(double *told)
3128 { /* initialisation (propagation of constant blocks outputs) */
3129         /*     Copyright INRIA */
3130
3131         int i2 = 0;
3132         scicos_flag flag = 0;
3133         int i = 0,j = 0;
3134         int ierr1 = 0;
3135
3136         ScicosImport *scs_imp = NULL;
3137         int *kf = NULL;
3138
3139         scs_imp = getscicosimportptr();
3140
3141         flag = 1;
3142         for (j=0;j<*(scs_imp->niord);j++) {
3143                 kf = &scs_imp->iord[j];
3144                 C2F(curblk).kfun = *kf; /* */
3145                 if (scs_imp->funtyp[*kf-1] > -1) {
3146                         /* continuous state */
3147                         if(scs_imp->blocks[*kf-1].nx > 0) {
3148                                 scs_imp->blocks[*kf-1].x  = &scs_imp->x[scs_imp->xptr[*kf-1]-1];
3149                                 scs_imp->blocks[*kf-1].xd = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3150                         }
3151                         scs_imp->blocks[*kf-1].nevprt = scs_imp->iord[j + *(scs_imp->niord)];
3152                         /*callf(told, xd, x, x,x,&flag);*/
3153                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3154                         if (flag < 0) {
3155                                 *ierr = 5 - flag;
3156                                 return;
3157                         }
3158                 }
3159                 if (scs_imp->blocks[*kf-1].nevout > 0) {
3160                         if (scs_imp->funtyp[*kf-1] < 0) {
3161                                 i = synchro_nev(scs_imp,*kf,ierr);
3162                                 if (*ierr != 0) {
3163                                         return;
3164                                 }
3165                                 i2 = i + scs_imp->clkptr[*kf-1] - 1;
3166                                 putevs(told, &i2, &ierr1);
3167                                 if (ierr1 != 0) {
3168                                         /* event conflict */
3169                                         *ierr = 3;
3170                                         return;
3171                                 }
3172                                 doit(told);
3173                                 if (*ierr != 0) {
3174                                         return;
3175                                 }
3176                         }
3177                 }
3178         }
3179 } /* idoit_ */
3180 /*--------------------------------------------------------------------------*/
3181 /* Subroutine doit */
3182 static void doit(double *told)
3183 { /* propagation of blocks outputs on discrete activations */
3184         /*     Copyright INRIA */
3185
3186         int i = 0,i2 = 0;
3187         scicos_flag flag = 0;
3188         int nord = 0;
3189         int ierr1 = 0;
3190         int ii = 0, kever = 0;
3191
3192         ScicosImport *scs_imp = NULL;
3193         int *kf = NULL;
3194
3195         scs_imp = getscicosimportptr();
3196
3197         kever = *pointi;
3198         *pointi = evtspt[kever];
3199         evtspt[kever] = -1;
3200
3201         nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever-1];
3202         if (nord == 0) {
3203                 return;
3204         }
3205
3206         for (ii = scs_imp->ordptr[kever-1]; ii <=scs_imp->ordptr[kever] - 1 ; ii++) {
3207                 kf = &scs_imp->ordclk[ii-1];
3208                 C2F(curblk).kfun = *kf;
3209                 if (scs_imp->funtyp[*kf-1] > -1) {
3210                         /* continuous state */
3211                         if(scs_imp->blocks[*kf-1].nx > 0) {
3212                                 scs_imp->blocks[*kf-1].x  = &scs_imp->x[scs_imp->xptr[*kf-1]-1];
3213                                 scs_imp->blocks[*kf-1].xd = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3214                         }
3215                         scs_imp->blocks[*kf-1].nevprt = abs(scs_imp->ordclk[ii + *(scs_imp->nordclk) - 1]);
3216                         flag = 1;
3217                         /*callf(told, xd, x, x,x,&flag);*/
3218                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3219                         if (flag < 0) {
3220                                 *ierr = 5 - flag;
3221                                 return;
3222                         }
3223                 }
3224
3225                 /* Initialize tvec */
3226                 if (scs_imp->blocks[*kf-1].nevout > 0) {
3227                         if (scs_imp->funtyp[*kf-1] < 0) {
3228                                 i = synchro_nev(scs_imp,*kf,ierr);
3229                                 if (*ierr != 0) {
3230                                         return;
3231                                 }
3232                                 i2 = i + scs_imp->clkptr[*kf-1] - 1;
3233                                 putevs(told, &i2, &ierr1);
3234                                 if (ierr1 != 0) {
3235                                         /* event conflict */
3236                                         *ierr = 3;
3237                                         return;
3238                                 }
3239                                 doit(told);
3240                                 if (*ierr != 0) {
3241                                         return;
3242                                 }
3243                         }
3244                 }
3245         }
3246 } /* doit_ */
3247 /*--------------------------------------------------------------------------*/
3248 /* Subroutine cdoit */
3249 static void cdoit(double *told)
3250 { /* propagation of continuous blocks outputs */
3251         /*     Copyright INRIA */
3252
3253         int i2 = 0;
3254         scicos_flag flag = 0;
3255         int ierr1 = 0;
3256         int i = 0,j = 0;
3257
3258         ScicosImport *scs_imp = NULL;
3259         int *kf = NULL;
3260
3261         scs_imp = getscicosimportptr();
3262
3263         /* Function Body */
3264         for (j=0;j<*(scs_imp->ncord);j++) {
3265                 kf = &scs_imp->cord[j];
3266                 C2F(curblk).kfun = *kf;
3267                 /* continuous state */
3268                 if(scs_imp->blocks[*kf-1].nx > 0) {
3269                         scs_imp->blocks[*kf-1].x  = &scs_imp->x[scs_imp->xptr[*kf-1]-1];
3270                         scs_imp->blocks[*kf-1].xd = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3271                 }
3272                 scs_imp->blocks[*kf-1].nevprt = scs_imp->cord[j + *(scs_imp->ncord)];
3273                 if (scs_imp->funtyp[*kf-1] > -1) {
3274                         flag = 1;
3275                         /*callf(told, xd, x, x,x,&flag);*/
3276                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3277                         if (flag < 0) {
3278                                 *ierr = 5 - flag;
3279                                 return;
3280                         }
3281                 }
3282
3283                 /* Initialize tvec */
3284                 if (scs_imp->blocks[*kf-1].nevout > 0) {
3285                         if (scs_imp->funtyp[*kf-1] < 0) {
3286                                 i = synchro_nev(scs_imp,*kf,ierr);
3287                                 if (*ierr != 0) {
3288                                         return;
3289                                 }
3290                                 i2 = i + scs_imp->clkptr[*kf-1] - 1;
3291                                 putevs(told, &i2, &ierr1);
3292                                 if (ierr1 != 0) {
3293                                         /* event conflict */
3294                                         *ierr = 3;
3295                                         return;
3296                                 }
3297                                 doit(told);
3298                                 if (*ierr != 0) {
3299                                         return;
3300                                 }
3301                         }
3302                 }
3303         }
3304 } /* cdoit_ */
3305 /*--------------------------------------------------------------------------*/
3306 /* Subroutine ddoit */
3307 static void ddoit(double *told)
3308 { /* update states & event out on discrete activations */
3309         /*     Copyright INRIA */
3310
3311         int i2 = 0,j = 0;
3312         scicos_flag flag = 0;
3313         int kiwa = 0;
3314         int i = 0,i3 = 0,ierr1 = 0;
3315         int ii = 0, keve = 0;
3316
3317         ScicosImport *scs_imp = NULL;
3318         int *kf = NULL;
3319
3320         scs_imp = getscicosimportptr();
3321
3322         /* Function Body */
3323         kiwa = 0;
3324         edoit(told,&kiwa);
3325         if (*ierr != 0) {
3326                 return;
3327         }
3328
3329         /* update continuous and discrete states on event */
3330         if (kiwa == 0) {
3331                 return;
3332         }
3333         for (i=0;i<kiwa;i++) {
3334                 keve = iwa[i];
3335                 if(critev[keve]!= 0){
3336                         hot = 0;
3337                 }
3338                 i2 = scs_imp->ordptr[keve] - 1;
3339                 for (ii = scs_imp->ordptr[keve-1]; ii <= i2; ii++) {
3340                         kf = &scs_imp->ordclk[ii-1];
3341                         C2F(curblk).kfun = *kf;
3342                         /* continuous state */
3343                         if(scs_imp->blocks[*kf-1].nx > 0) {
3344                                 scs_imp->blocks[*kf-1].x  = &scs_imp->x[scs_imp->xptr[*kf-1]-1];
3345                                 scs_imp->blocks[*kf-1].xd = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3346                         }
3347
3348                         scs_imp->blocks[*kf-1].nevprt=scs_imp->ordclk[ii + *(scs_imp->nordclk) - 1];
3349
3350                         if (scs_imp->blocks[*kf-1].nevout > 0) {
3351                                 if (scs_imp->funtyp[*kf-1] >= 0) {
3352                                         /* initialize evout */
3353                                         for (j=0;j<scs_imp->blocks[*kf-1].nevout;j++) {
3354                                                 scs_imp->blocks[*kf-1].evout[j]=-1;
3355                                         }
3356                                         flag = 3;
3357
3358                                         if(scs_imp->blocks[*kf-1].nevprt>0) { /* if event has continuous origin don't call*/
3359                                                 /*callf(told, xd, x, x ,x,&flag);*/
3360                                                 callf(told, &scs_imp->blocks[*kf-1], &flag);
3361                                                 if (flag < 0) {
3362                                                         *ierr = 5 - flag;
3363                                                         return;
3364                                                 }
3365                                         }
3366
3367                                         for (j=0;j<scs_imp->blocks[*kf-1].nevout;j++) {
3368                                                 if (scs_imp->blocks[*kf-1].evout[j] >= 0.) {
3369                                                         i3 = j + scs_imp->clkptr[*kf-1] ;
3370                                                         addevs(scs_imp->blocks[*kf-1].evout[j]+(*told), &i3, &ierr1);
3371                                                         if (ierr1 != 0) {
3372                                                                 /* event conflict */
3373                                                                 *ierr = 3;
3374                                                                 return;
3375                                                         }
3376                                                 }
3377                                         }
3378                                 }
3379                         }
3380
3381                         if(scs_imp->blocks[*kf-1].nevprt> 0) {
3382                                 if (scs_imp->blocks[*kf-1].nx+scs_imp->blocks[*kf-1].nz+scs_imp->blocks[*kf-1].noz > 0 || \
3383                                         *scs_imp->blocks[*kf-1].work !=NULL) {
3384                                                 /*  if a hidden state exists, must also call (for new scope eg)  */
3385                                                 /*  to avoid calling non-real activations */
3386                                                 flag = 2;
3387                                                 /*callf(told, xd, x, x,x,&flag);*/
3388                                                 callf(told, &scs_imp->blocks[*kf-1], &flag);
3389                                                 if (flag < 0) {
3390                                                         *ierr = 5 - flag;
3391                                                         return;
3392                                                 }
3393                                 }
3394                         }
3395                         else {
3396                                 if (*scs_imp->blocks[*kf-1].work !=NULL) {
3397                                         flag = 2;
3398                                         scs_imp->blocks[*kf-1].nevprt = 0; /* in case some hidden continuous blocks need updating */
3399                                         /*callf(told, xd, x, x,x,&flag);*/
3400                                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3401                                         if (flag < 0) {
3402                                                 *ierr = 5 - flag;
3403                                                 return;
3404                                         }
3405                                 }
3406                         }
3407                 }
3408         }
3409 } /* ddoit_ */
3410 /*--------------------------------------------------------------------------*/
3411 /* Subroutine edoit */
3412 static void edoit(double *told, int *kiwa)
3413 { /* update blocks output on discrete activations */
3414         /*     Copyright INRIA */
3415
3416         int i2 = 0;
3417         scicos_flag flag = 0;
3418         int ierr1 = 0, i = 0;
3419         int kever = 0, ii = 0;
3420
3421         ScicosImport *scs_imp = NULL;
3422         int *kf = NULL;
3423         int nord = 0;
3424
3425         scs_imp = getscicosimportptr();
3426
3427         /* Function Body */
3428         kever = *pointi;
3429
3430         *pointi = evtspt[kever];
3431         evtspt[kever] = -1;
3432
3433         nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever-1];
3434         if (nord == 0) {
3435                 return;
3436         }
3437         iwa[*kiwa] = kever;
3438         ++(*kiwa);
3439         for (ii = scs_imp->ordptr[kever-1]; ii <= scs_imp->ordptr[kever] - 1; ii++) {
3440                 kf = &scs_imp->ordclk[ii-1];
3441                 C2F(curblk).kfun = *kf;
3442
3443                 if (scs_imp->funtyp[*kf-1] > -1) {
3444                         /* continuous state */
3445                         if(scs_imp->blocks[*kf-1].nx > 0) {
3446                                 scs_imp->blocks[*kf-1].x  = &scs_imp->x[scs_imp->xptr[*kf-1]-1];
3447                                 scs_imp->blocks[*kf-1].xd = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3448                         }
3449
3450                         scs_imp->blocks[*kf-1].nevprt = abs(scs_imp->ordclk[ii + *(scs_imp->nordclk) - 1]);
3451
3452                         flag = 1;
3453                         /*callf(told, xd, x, x,x,&flag);*/
3454                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3455                         if (flag < 0) {
3456                                 *ierr = 5 - flag;
3457                                 return;
3458                         }
3459                 }
3460
3461                 /* Initialize tvec */
3462                 if (scs_imp->blocks[*kf-1].nevout > 0) {
3463                         if (scs_imp->funtyp[*kf-1] < 0) {
3464                                 i = synchro_nev(scs_imp,*kf,ierr);
3465                                 if (*ierr != 0) {
3466                                         return;
3467                                 }
3468                                 i2 = i + scs_imp->clkptr[*kf-1] - 1;
3469                                 putevs(told, &i2, &ierr1);
3470                                 if (ierr1 != 0) {
3471                                         /* event conflict */
3472                                         *ierr = 3;
3473                                         return;
3474                                 }
3475                                 edoit(told,kiwa);
3476                                 if (*ierr != 0) {
3477                                         return;
3478                                 }
3479                         }
3480                 }
3481         }
3482 } /* edoit_ */
3483 /*--------------------------------------------------------------------------*/
3484 /* Subroutine odoit */
3485 static void odoit(double *told, double *xt, double *xtd, double *residual)
3486 { /* update blocks derivative of continuous time block */
3487         /*     Copyright INRIA */
3488
3489         int i2 = 0;
3490         scicos_flag flag = 0;
3491         int keve = 0, kiwa = 0;
3492         int ierr1 = 0, i = 0;
3493         int ii = 0, jj = 0;
3494
3495         ScicosImport *scs_imp = NULL;
3496         int *kf = NULL;
3497
3498         scs_imp = getscicosimportptr();
3499
3500         /* Function Body */
3501         kiwa = 0;
3502         for (jj = 0; jj < *(scs_imp->noord); jj++) {
3503                 kf = &scs_imp->oord[jj];
3504                 C2F(curblk).kfun = *kf;
3505                 /* continuous state */
3506                 if(scs_imp->blocks[*kf-1].nx > 0) {
3507                         scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
3508                         scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
3509                         scs_imp->blocks[*kf-1].res = &residual[scs_imp->xptr[*kf-1]-1];
3510                 }
3511
3512                 scs_imp->blocks[*kf-1].nevprt = scs_imp->oord[jj + *(scs_imp->noord)];
3513                 if (scs_imp->funtyp[*kf-1] > -1) {
3514                         flag = 1;
3515                         /*callf(told, xtd, xt, residual,x,&flag);*/
3516                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3517                         if (flag < 0) {
3518                                 *ierr = 5 - flag;
3519                                 return;
3520                         }
3521                 }
3522
3523                 if (scs_imp->blocks[*kf-1].nevout > 0) {
3524                         if (scs_imp->funtyp[*kf-1] < 0) {
3525                                 if(scs_imp->blocks[*kf-1].nmode > 0) {
3526                                         i2 = scs_imp->blocks[*kf-1].mode[0] + scs_imp->clkptr[*kf-1] - 1;
3527                                 }
3528                                 else {
3529                                         i = synchro_nev(scs_imp,*kf,ierr);
3530                                         if (*ierr != 0) {
3531                                                 return;
3532                                         }
3533                                         i2 = i + scs_imp->clkptr[*kf-1] - 1;
3534                                 }
3535                                 putevs(told, &i2, &ierr1);
3536                                 if (ierr1 != 0) {
3537                                         /* event conflict */
3538                                         *ierr = 3;
3539                                         return;
3540                                 }
3541                                 ozdoit(told, xt, xtd, &kiwa);
3542                                 if (*ierr != 0) {
3543                                         return;
3544                                 }
3545                         }
3546                 }
3547         }
3548
3549         /*  update states derivatives */
3550         for (ii = 0; ii < *(scs_imp->noord); ii++) {
3551                 kf = &scs_imp->oord[ii];
3552                 C2F(curblk).kfun = *kf;
3553                 if (scs_imp->blocks[*kf-1].nx > 0 || \
3554                         *scs_imp->blocks[*kf-1].work !=NULL) {
3555                                 /* work tests if a hidden state exists, used for delay block */
3556                                 flag = 0;
3557                                 /* continuous state */
3558                                 if(scs_imp->blocks[*kf-1].nx > 0) {
3559                                         scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
3560                                         scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
3561                                         scs_imp->blocks[*kf-1].res = &residual[scs_imp->xptr[*kf-1]-1];
3562                                 }
3563                                 scs_imp->blocks[*kf-1].nevprt = scs_imp->oord[ii + *(scs_imp->noord)];
3564                                 /*callf(told, xtd, xt, residual,xt,&flag);*/
3565                                 callf(told, &scs_imp->blocks[*kf-1], &flag);
3566
3567                                 if (flag < 0) {
3568                                         *ierr = 5 - flag;
3569                                         return;
3570                                 }
3571                 }
3572         }
3573
3574         for (i = 0; i < kiwa; i++) {
3575                 keve = iwa[i];
3576                 for (ii = scs_imp->ordptr[keve-1]; ii <= scs_imp->ordptr[keve] - 1; ii++) {
3577                         kf = &scs_imp->ordclk[ii-1];
3578                         C2F(curblk).kfun = *kf;
3579                         if (scs_imp->blocks[*kf-1].nx > 0 || \
3580                                 *scs_imp->blocks[*kf-1].work !=NULL) {
3581                                         /* work tests if a hidden state exists */
3582                                         flag = 0;
3583                                         /* continuous state */
3584                                         if(scs_imp->blocks[*kf-1].nx > 0) {
3585                                                 scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
3586                                                 scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
3587                                                 scs_imp->blocks[*kf-1].res = &residual[scs_imp->xptr[*kf-1]-1];
3588                                         }
3589                                         scs_imp->blocks[*kf-1].nevprt = abs(scs_imp->ordclk[ii + *(scs_imp->nordclk) - 1]);
3590                                         /*callf(told, xtd, xt, residual,xt,&flag);*/
3591                                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3592
3593                                         if (flag < 0) {
3594                                                 *ierr = 5 - flag;
3595                                                 return;
3596                                         }
3597                         }
3598                 }
3599         }
3600 } /* odoit_ */
3601 /*--------------------------------------------------------------------------*/
3602 /* Subroutine reinitdoit */
3603 static void reinitdoit(double *told)
3604 { /* update blocks xproperties of continuous time block */
3605         /*     Copyright INRIA */
3606
3607         int i2 = 0;
3608         scicos_flag flag = 0;
3609         int keve = 0, kiwa = 0;
3610         int ierr1 = 0, i = 0;
3611         int ii = 0, jj = 0;
3612
3613         ScicosImport *scs_imp = NULL;
3614         int *kf = NULL;
3615
3616         scs_imp = getscicosimportptr();
3617
3618         /* Function Body */
3619         kiwa = 0;
3620         for (jj = 0; jj < *(scs_imp->noord); jj++) {
3621                 kf = &scs_imp->oord[jj];
3622                 C2F(curblk).kfun = *kf;
3623                 /* continuous state */
3624                 if(scs_imp->blocks[*kf-1].nx > 0) {
3625                         scs_imp->blocks[*kf-1].x  = &scs_imp->x[scs_imp->xptr[*kf-1]-1];
3626                         scs_imp->blocks[*kf-1].xd = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3627                 }
3628                 scs_imp->blocks[*kf-1].nevprt = scs_imp->oord[jj + *(scs_imp->noord)];
3629                 if (scs_imp->funtyp[*kf-1] > -1) {
3630                         flag = 1;
3631                         /*callf(told, xd, x, x,x,&flag);*/
3632                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3633                         if (flag < 0) {
3634                                 *ierr = 5 - flag;
3635                                 return;
3636                         }
3637                 }
3638
3639                 if (scs_imp->blocks[*kf-1].nevout > 0 && scs_imp->funtyp[*kf-1] < 0) {
3640                         i = synchro_nev(scs_imp,*kf,ierr);
3641                         if (*ierr != 0) {
3642                                 return;
3643                         }
3644                         if(scs_imp->blocks[*kf-1].nmode>0) {
3645                                 scs_imp->blocks[*kf-1].mode[0] = i;
3646                         }
3647                         i2 = i + scs_imp->clkptr[*kf-1] - 1;
3648                         putevs(told, &i2, &ierr1);
3649                         if (ierr1 != 0) {
3650                                 /* event conflict */
3651                                 *ierr = 3;
3652                                 return;
3653                         }
3654                         doit(told);
3655                         if (*ierr != 0) {
3656                                 return;
3657                         }
3658                 }
3659         }
3660
3661         /* re-initialize */
3662         for (ii=0;ii<*(scs_imp->noord);ii++) {
3663                 kf = &scs_imp->oord[ii];
3664                 C2F(curblk).kfun = *kf;
3665                 if (scs_imp->blocks[*kf-1].nx > 0) {
3666                         flag = 7;
3667                         scs_imp->blocks[*kf-1].x  = &scs_imp->x[scs_imp->xptr[*kf-1]-1];
3668                         scs_imp->blocks[*kf-1].xd = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3669                         scs_imp->blocks[*kf-1].res = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3670                         scs_imp->blocks[*kf-1].nevprt = scs_imp->oord[ii + *(scs_imp->noord)];
3671                         scs_imp->blocks[*kf-1].xprop = &scs_imp->xprop[-1+scs_imp->xptr[*kf-1]];
3672                         /*callf(told, xd, x, xd,x,&flag);*/
3673                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3674                         if (flag < 0) {
3675                                 *ierr = 5 - flag;
3676                                 return;
3677                         }
3678                 }
3679         }
3680
3681         for (i=0;i<kiwa;i++) {
3682                 keve = iwa[i];
3683                 for (ii = scs_imp->ordptr[keve-1]; ii <= scs_imp->ordptr[keve] - 1; ii++) {
3684                         kf = &scs_imp->ordclk[ii-1];
3685                         C2F(curblk).kfun = *kf;
3686                         if (scs_imp->blocks[*kf-1].nx > 0) {
3687                                 flag = 7;
3688                                 scs_imp->blocks[*kf-1].x  = &scs_imp->x[scs_imp->xptr[*kf-1]-1];
3689                                 scs_imp->blocks[*kf-1].xd = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3690                                 scs_imp->blocks[*kf-1].res = &scs_imp->xd[scs_imp->xptr[*kf-1]-1];
3691                                 scs_imp->blocks[*kf-1].nevprt = abs(scs_imp->ordclk[ii + *(scs_imp->nordclk) - 1]);
3692                                 scs_imp->blocks[*kf-1].xprop = &scs_imp->xprop[-1+scs_imp->xptr[*kf-1]];
3693                                 /*callf(told, xd, x, xd,x,&flag);*/
3694                                 callf(told, &scs_imp->blocks[*kf-1], &flag);
3695                                 if (flag < 0) {
3696                                         *ierr = 5 - flag;
3697                                         return;
3698                                 }
3699                         }
3700                 }
3701         }
3702 } /* reinitdoit_ */
3703 /*--------------------------------------------------------------------------*/
3704 /* Subroutine ozdoit */
3705 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa)
3706 { /* update blocks output of continuous time block on discrete activations */
3707         /*     Copyright INRIA */
3708
3709         int i2 = 0;
3710         scicos_flag flag = 0;
3711         int nord = 0;
3712         int ierr1 = 0, i = 0;
3713         int ii = 0, kever = 0;
3714
3715         ScicosImport *scs_imp = NULL;
3716         int *kf = NULL;
3717
3718         scs_imp = getscicosimportptr();
3719
3720         /* Function Body */
3721         kever = *pointi;
3722         *pointi = evtspt[kever];
3723         evtspt[kever] = -1;
3724
3725         nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever-1];
3726         if (nord == 0) {
3727                 return;
3728         }
3729         iwa[*kiwa] = kever;
3730         ++(*kiwa);
3731
3732         for (ii = scs_imp->ordptr[kever-1]; ii <= scs_imp->ordptr[kever] - 1; ii++) {
3733                 kf = &scs_imp->ordclk[ii-1];
3734                 C2F(curblk).kfun=*kf;
3735                 if (scs_imp->funtyp[*kf-1] > -1) {
3736                         /* continuous state */
3737                         if(scs_imp->blocks[*kf-1].nx > 0) {
3738                                 scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
3739                                 scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
3740                         }
3741                         scs_imp->blocks[*kf-1].nevprt = abs(scs_imp->ordclk[ii + *(scs_imp->nordclk) - 1]);
3742                         flag = 1;
3743                         /*callf(told, xtd, xt, xt,x,&flag);*/
3744                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3745                         if (flag < 0) {
3746                                 *ierr = 5 - flag;
3747                                 return;
3748                         }
3749                 }
3750
3751                 /* Initialize tvec */
3752                 if (scs_imp->blocks[*kf-1].nevout > 0) {
3753                         if (scs_imp->funtyp[*kf-1] < 0) {
3754                                 if (phase==1 || scs_imp->blocks[*kf-1].nmode==0) {
3755                                         i = synchro_nev(scs_imp,*kf,ierr);
3756                                         if (*ierr != 0) {
3757                                                 return;
3758                                         }
3759                                 }
3760                                 else {
3761                                         i = scs_imp->blocks[*kf-1].mode[0];
3762                                 }
3763                                 i2 = i + scs_imp->clkptr[*kf-1] - 1;
3764                                 putevs(told, &i2, &ierr1);
3765                                 if (ierr1 != 0) {
3766                                         /* event conflict */
3767                                         *ierr = 3;
3768                                         return;
3769                                 }
3770                                 ozdoit(told, xt, xtd, kiwa);
3771                         }
3772                 }
3773         }
3774 } /* ozdoit_ */
3775 /*--------------------------------------------------------------------------*/
3776 /* Subroutine zdoit */
3777 static void zdoit(double *told, double *xt, double *xtd, double *g)
3778 { /* update blocks zcross of continuous time block  */
3779         /*     Copyright INRIA */
3780         int i2 = 0;
3781         scicos_flag flag = 0;
3782         int keve = 0, kiwa = 0;
3783         int ierr1 = 0, i = 0,j = 0;
3784         int ii = 0, jj = 0;
3785
3786         ScicosImport *scs_imp = NULL;
3787         int *kf = NULL;
3788
3789         scs_imp = getscicosimportptr();
3790
3791         /* Function Body */
3792         for(i=0;i<*(scs_imp->ng);i++) {
3793                 g[i]=0.;
3794         }
3795
3796         kiwa = 0;
3797         for (jj = 0; jj < *(scs_imp->nzord); jj++) {
3798                 kf = &scs_imp->zord[jj];
3799                 C2F(curblk).kfun = *kf;
3800                 /* continuous state */
3801                 if(scs_imp->blocks[*kf-1].nx > 0) {
3802                         scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
3803                         scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
3804                 }
3805                 scs_imp->blocks[*kf-1].nevprt = scs_imp->zord[jj + *(scs_imp->nzord)];
3806
3807                 if (scs_imp->funtyp[*kf-1] > -1) {
3808                         flag = 1;
3809                         /*callf(told, xtd, xt, xt,xt,&flag);*/
3810                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3811                         if (flag < 0) {
3812                                 *ierr = 5 - flag;
3813                                 return;
3814                         }
3815                 }
3816
3817                 /* Initialize tvec */
3818                 if (scs_imp->blocks[*kf-1].nevout > 0) {
3819                         if (scs_imp->funtyp[*kf-1] < 0) {
3820                                 if (phase==1 || scs_imp->blocks[*kf-1].nmode==0) {
3821                                         i = synchro_nev(scs_imp,*kf,ierr);
3822                                         if (*ierr != 0) {
3823                                                 return;
3824                                         }
3825                                 }
3826                                 else {
3827                                         i = scs_imp->blocks[*kf-1].mode[0];
3828                                 }
3829                                 i2 = i + scs_imp->clkptr[*kf-1] - 1;
3830                                 putevs(told, &i2, &ierr1);
3831                                 if (ierr1 != 0) {
3832                                         /* event conflict */
3833                                         *ierr = 3;
3834                                         return;
3835                                 }
3836                                 ozdoit(told, xt, xtd, &kiwa);
3837                                 if (*ierr != 0) {
3838                                         return;
3839                                 }
3840                         }
3841                 }
3842         }
3843
3844         /* update zero crossing surfaces */
3845         for (ii = 0; ii < *(scs_imp->nzord); ii++) {
3846                 kf = &scs_imp->zord[ii];
3847                 C2F(curblk).kfun = *kf;
3848                 if (scs_imp->blocks[*kf-1].ng > 0) {
3849                         /* update g array ptr */
3850                         scs_imp->blocks[*kf-1].g = &g[scs_imp->zcptr[*kf-1]-1];
3851                         if (scs_imp->funtyp[*kf-1] > 0) {
3852                                 flag = 9;
3853                                 /* continuous state */
3854                                 if(scs_imp->blocks[*kf-1].nx > 0) {
3855                                         scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
3856                                         scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
3857                                 }
3858                                 scs_imp->blocks[*kf-1].nevprt = scs_imp->zord[ii + *(scs_imp->nzord)];
3859                                 /*callf(told, xtd, xt, xtd,g,&flag);*/
3860                                 callf(told, &scs_imp->blocks[*kf-1], &flag);
3861                                 if (flag < 0) {
3862                                         *ierr = 5 - flag;
3863                                         return;
3864                                 }
3865                         }
3866                         else {
3867                                 j = synchro_g_nev(scs_imp,g,*kf,ierr);
3868                                 if (*ierr != 0) {
3869                                         return;
3870                                 }
3871                                 if( (phase==1) && (scs_imp->blocks[*kf-1].nmode>0) ) {
3872                                         scs_imp->blocks[*kf-1].mode[0]= j;
3873                                 }
3874                         }
3875
3876                         // scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
3877
3878                 }
3879         }
3880
3881         for (i=0;i<kiwa;i++) {
3882                 keve = iwa[i];
3883                 for (ii = scs_imp->ordptr[keve-1]; ii <= scs_imp->ordptr[keve] - 1; ii++) {
3884                         kf = &scs_imp->ordclk[ii-1];
3885                         C2F(curblk).kfun = *kf;
3886                         if (scs_imp->blocks[*kf-1].ng > 0) {
3887                                 /* update g array ptr */
3888                                 scs_imp->blocks[*kf-1].g = &g[scs_imp->zcptr[*kf-1]-1];
3889                                 if (scs_imp->funtyp[*kf-1] > 0) {
3890                                         flag = 9;
3891                                         /* continuous state */
3892                                         if(scs_imp->blocks[*kf-1].nx > 0) {
3893                                                 scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
3894                                                 scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
3895                                         }
3896                                         scs_imp->blocks[*kf-1].nevprt = abs(scs_imp->ordclk[ii + *(scs_imp->nordclk) -1]);
3897                                         /*callf(told, xtd, xt, xtd,g,&flag);*/
3898                                         callf(told, &scs_imp->blocks[*kf-1], &flag);
3899                                         if (flag < 0) {
3900                                                 *ierr = 5 - flag;
3901                                                 return;
3902                                         }
3903                                 }
3904                                 else {
3905                                         j = synchro_g_nev(scs_imp,g,*kf,ierr);
3906                                         if (*ierr != 0) {
3907                                                 return;
3908                                         }
3909                                         if((phase==1) && (scs_imp->blocks[*kf-1].nmode>0)) {
3910                                                 scs_imp->blocks[*kf-1].mode[0]= j;
3911                                         }
3912                                 }
3913
3914                                 //scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
3915                         }
3916                 }
3917         }
3918 } /* zdoit_ */
3919 /*--------------------------------------------------------------------------*/
3920 /* Subroutine Jdoit */
3921 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job)
3922 { /* update blocks jacobian of continuous time block  */
3923         /*     Copyright INRIA */
3924
3925         int i2 = 0;
3926         scicos_flag flag = 0;
3927         int keve = 0, kiwa = 0;
3928         int ierr1 = 0, i = 0;
3929         int ii = 0, jj = 0;
3930
3931         ScicosImport *scs_imp = NULL;
3932         int *kf = NULL;
3933
3934         scs_imp = getscicosimportptr();
3935
3936         /* Function Body */
3937         kiwa = 0;
3938         for (jj = 0; jj < *(scs_imp->noord); jj++) {
3939                 kf = &scs_imp->oord[jj];
3940                 C2F(curblk).kfun = *kf;
3941                 scs_imp->blocks[*kf-1].nevprt = scs_imp->oord[jj + *(scs_imp->noord)];
3942                 if (scs_imp->funtyp[*kf-1] > -1) {
3943                         flag = 1;
3944                         /* applying desired output */
3945                         if ((*job==2)&&(scs_imp->oord[jj]==AJacobian_block)) {
3946                         }
3947                         else
3948                                 /* continuous state */
3949                                 if(scs_imp->blocks[*kf-1].nx > 0) {
3950                                         scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
3951                                         scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
3952                                         scs_imp->blocks[*kf-1].res = &residual[scs_imp->xptr[*kf-1]-1];
3953                                 }
3954
3955                                 /*callf(told, xtd, xt, residual,x,&flag);*/
3956                                 callf(told, &scs_imp->blocks[*kf-1], &flag);
3957                                 if (flag < 0) {
3958                                         *ierr = 5 - flag;
3959                                         return;
3960                                 }
3961                 }
3962
3963                 if (scs_imp->blocks[*kf-1].nevout > 0) {
3964                         if (scs_imp->funtyp[*kf-1] < 0) {
3965                                 if(scs_imp->blocks[*kf-1].nmode > 0){
3966                                         i2 = scs_imp->blocks[*kf-1].mode[0] + scs_imp->clkptr[*kf-1] - 1;
3967                                 }
3968                                 else {
3969                                         i = synchro_nev(scs_imp,*kf,ierr);
3970                                         if (*ierr != 0) {
3971                                                 return;
3972                                         }
3973                                         i2 = i + scs_imp->clkptr[*kf-1] - 1;
3974                                 }
3975                                 putevs(told, &i2, &ierr1);
3976                                 if (ierr1 != 0) {
3977                                         /* event conflict */
3978                                         *ierr = 3;
3979                                         return;
3980                                 }
3981                                 ozdoit(told, xt, xtd, &kiwa);
3982                         }
3983                 }
3984         }
3985
3986         /* update states derivatives */
3987         for (ii = 0; ii < *(scs_imp->noord); ii++) {
3988                 kf = &scs_imp->oord[ii];
3989                 C2F(curblk).kfun = *kf;
3990                 if (scs_imp->blocks[*kf-1].nx > 0 || \
3991                         *scs_imp->blocks[*kf-1].work !=NULL) {
3992                                 /* work tests if a hidden state exists, used for delay block */
3993                                 flag = 0;
3994                                 if (((*job==1)&&(scs_imp->oord[ii]==AJacobian_block)) || (*job!=1)) {
3995                                         if (*job==1)  flag = 10;
3996                                         /* continuous state */
3997                                         if(scs_imp->blocks[*kf-1].nx > 0) {
3998                                                 scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
3999                                                 scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
4000                                                 scs_imp->blocks[*kf-1].res = &residual[scs_imp->xptr[*kf-1]-1];
4001                                         }
4002                                         scs_imp->blocks[*kf-1].nevprt = scs_imp->oord[ii + *(scs_imp->noord)];
4003                                         /*callf(told, xtd, xt, residual,xt,&flag);*/
4004                                         callf(told, &scs_imp->blocks[*kf-1], &flag);
4005                                 }
4006                                 if (flag < 0) {
4007                                         *ierr = 5 - flag;
4008                                         return;
4009                                 }
4010                 }
4011         }
4012
4013         for (i = 0; i < kiwa; i++) {
4014                 keve = iwa[i];
4015                 for (ii = scs_imp->ordptr[keve-1]; ii <= scs_imp->ordptr[keve] - 1;  ii++) {
4016                         kf = &scs_imp->ordclk[ii-1];
4017                         C2F(curblk).kfun = *kf;
4018                         if (scs_imp->blocks[*kf-1].nx > 0 || \
4019                                 *scs_imp->blocks[*kf-1].work !=NULL) {
4020                                         /* work tests if a hidden state exists */
4021                                         flag = 0;
4022                                         if (((*job==1)&&(scs_imp->oord[ii-1]==AJacobian_block))||(*job!=1)) {
4023                                                 if (*job==1)  flag = 10;
4024                                                 /* continuous state */
4025                                                 if(scs_imp->blocks[*kf-1].nx > 0) {
4026                                                         scs_imp->blocks[*kf-1].x  = &xt[scs_imp->xptr[*kf-1]-1];
4027                                                         scs_imp->blocks[*kf-1].xd = &xtd[scs_imp->xptr[*kf-1]-1];
4028                                                         scs_imp->blocks[*kf-1].res = &residual[scs_imp->xptr[*kf-1]-1];
4029                                                 }
4030                                                 scs_imp->blocks[*kf-1].nevprt = abs(scs_imp->ordclk[ii + *(scs_imp->nordclk) - 1]);
4031                                                 /*callf(told, xtd, xt, residual,xt,&flag);*/
4032                                                 callf(told, &scs_imp->blocks[*kf-1], &flag);
4033                                         }
4034                                         if (flag < 0) {
4035                                                 *ierr = 5 - flag;
4036                                                 return;
4037                                         }
4038                         }
4039                 }
4040         }
4041 } /* Jdoit_ */
4042 /*--------------------------------------------------------------------------*/
4043 /* Subroutine synchro_nev */
4044 static int synchro_nev(ScicosImport *scs_imp,int kf,int *ierr)
4045 { /* synchro blocks computation  */
4046         /*     Copyright INRIA */
4047         SCSREAL_COP *outtbdptr = NULL;     /*to store double of outtb*/
4048         SCSINT8_COP *outtbcptr = NULL;     /*to store int8 of outtb*/
4049         SCSINT16_COP *outtbsptr = NULL;    /*to store int16 of outtb*/
4050         SCSINT32_COP *outtblptr = NULL;    /*to store int32 of outtb*/
4051         SCSUINT8_COP *outtbucptr = NULL;   /*to store unsigned int8 of outtb */
4052         SCSUINT16_COP *outtbusptr = NULL;  /*to store unsigned int16 of outtb */
4053         SCSUINT32_COP *outtbulptr = NULL;  /*to store unsigned int32 of outtb */
4054
4055         int cond = 0;
4056         int i=0; /* return 0 by default */
4057
4058         /* variable for param */
4059         int *outtbtyp = 0;
4060         void **outtbptr = NULL;
4061         int *funtyp = 0;
4062         int *inplnk = 0;
4063         int *inpptr = 0;
4064
4065         /* get param ptr */
4066         outtbtyp = scs_imp->outtbtyp;
4067         outtbptr = scs_imp->outtbptr;
4068         funtyp = scs_imp->funtyp;
4069         inplnk = scs_imp->inplnk;
4070         inpptr = scs_imp->inpptr;
4071
4072         /* if-then-else blk */
4073         if (funtyp[kf-1] == -1) {
4074                 switch(outtbtyp[-1+inplnk[inpptr[kf-1]-1]])
4075                 {
4076                 case SCSREAL_N    : outtbdptr=(SCSREAL_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4077                         cond = (*outtbdptr <= 0.);
4078                         break;
4079
4080                 case SCSCOMPLEX_N : outtbdptr=(SCSCOMPLEX_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4081                         cond = (*outtbdptr <= 0.);
4082                         break;
4083
4084                 case SCSINT8_N    : outtbcptr=(SCSINT8_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4085                         cond = (*outtbcptr <= 0);
4086                         break;
4087
4088                 case SCSINT16_N   : outtbsptr=(SCSINT16_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4089                         cond = (*outtbsptr <= 0);
4090                         break;
4091
4092                 case SCSINT32_N   : outtblptr=(SCSINT32_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4093                         cond = (*outtblptr <= 0);
4094                         break;
4095
4096                 case SCSUINT8_N   : outtbucptr=(SCSUINT8_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4097                         cond = (*outtbucptr <= 0);
4098                         break;
4099
4100                 case SCSUINT16_N  : outtbusptr=(SCSUINT16_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4101                         cond = (*outtbusptr <= 0);
4102                         break;
4103
4104                 case SCSUINT32_N  : outtbulptr=(SCSUINT32_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4105                         cond = (*outtbulptr <= 0);
4106                         break;
4107
4108                 default  : /* Add a message here */
4109                         *ierr = 25;
4110                         return 0;
4111                         break;
4112                 }
4113                 if (cond) {
4114                         i=2;
4115                 }
4116                 else {
4117                         i=1;
4118                 }
4119         }
4120         /* eselect blk */
4121         else if (funtyp[kf-1] == -2) {
4122                 switch(outtbtyp[-1+inplnk[inpptr[kf-1]-1]])
4123                 {
4124                 case SCSREAL_N    : outtbdptr=(SCSREAL_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4125                         i=max(min((int) *outtbdptr,scs_imp->blocks[kf-1].nevout),1);
4126                         break;
4127
4128                 case SCSCOMPLEX_N : outtbdptr=(SCSCOMPLEX_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4129                         i=max(min((int) *outtbdptr,scs_imp->blocks[kf-1].nevout),1);
4130                         break;
4131
4132                 case SCSINT8_N    : outtbcptr=(SCSINT8_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4133                         i=max(min((int) *outtbcptr,
4134                                 scs_imp->blocks[kf-1].nevout),1);
4135                         break;
4136
4137                 case SCSINT16_N   : outtbsptr=(SCSINT16_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4138                         i=max(min((int) *outtbsptr,scs_imp->blocks[kf-1].nevout),1);
4139                         break;
4140
4141                 case SCSINT32_N   : outtblptr=(SCSINT32_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4142                         i=max(min((int) *outtblptr,scs_imp->blocks[kf-1].nevout),1);
4143                         break;
4144
4145                 case SCSUINT8_N   : outtbucptr=(SCSUINT8_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4146                         i=max(min((int) *outtbucptr,scs_imp->blocks[kf-1].nevout),1);
4147                         break;
4148
4149                 case SCSUINT16_N  : outtbusptr=(SCSUINT16_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4150                         i=max(min((int) *outtbusptr,scs_imp->blocks[kf-1].nevout),1);
4151                         break;
4152
4153                 case SCSUINT32_N  : outtbulptr=(SCSUINT32_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4154                         i=max(min((int) *outtbulptr,scs_imp->blocks[kf-1].nevout),1);
4155                         break;
4156
4157                 default  : /* Add a message here */
4158                         *ierr = 25;
4159                         return 0;
4160                         break;
4161                 }
4162         }
4163         return i;
4164 } /* synchro_nev */
4165 /*--------------------------------------------------------------------------*/
4166 /* Subroutine synchro_g_nev */
4167 static int synchro_g_nev(ScicosImport *scs_imp,double *g,int kf,int *ierr)
4168 { /* synchro blocks with zcross computation  */
4169         /*     Copyright INRIA */
4170         SCSREAL_COP *outtbdptr = NULL;     /*to store double of outtb*/
4171         SCSINT8_COP *outtbcptr = NULL;     /*to store int8 of outtb*/
4172         SCSINT16_COP *outtbsptr = NULL;    /*to store int16 of outtb*/
4173         SCSINT32_COP *outtblptr = NULL;    /*to store int32 of outtb*/
4174         SCSUINT8_COP *outtbucptr = NULL;   /*to store unsigned int8 of outtb */
4175         SCSUINT16_COP *outtbusptr = NULL;  /*to store unsigned int16 of outtb */
4176         SCSUINT32_COP *outtbulptr = NULL;  /*to store unsigned int32 of outtb */
4177
4178         int cond = 0;
4179         int i=0; /* return 0 by default */
4180         int jj=0;
4181
4182         /* variable for param */
4183         int *outtbtyp = NULL;
4184         void **outtbptr = NULL;
4185         int *funtyp = NULL;
4186         int *inplnk = NULL;
4187         int *inpptr = NULL;
4188         int *zcptr = NULL;
4189
4190         /* get param ptr */
4191         outtbtyp = scs_imp->outtbtyp;
4192         outtbptr = scs_imp->outtbptr;
4193         funtyp = scs_imp->funtyp;
4194         inplnk = scs_imp->inplnk;
4195         inpptr = scs_imp->inpptr;
4196         zcptr = scs_imp->zcptr;
4197
4198         /* if-then-else blk */
4199         if (funtyp[kf-1] == -1) {
4200                 switch(outtbtyp[-1+inplnk[inpptr[kf-1]-1]])
4201                 {
4202                 case SCSREAL_N    : outtbdptr=(SCSREAL_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4203                         g[zcptr[kf-1]-1]=*outtbdptr;
4204                         cond = (*outtbdptr <= 0.);
4205                         break;
4206
4207                 case SCSCOMPLEX_N : outtbdptr=(SCSCOMPLEX_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4208                         g[zcptr[kf-1]-1]=*outtbdptr;
4209                         cond = (*outtbdptr <= 0.);
4210                         break;
4211
4212                 case SCSINT8_N    : outtbcptr=(SCSINT8_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4213                         g[zcptr[kf-1]-1]=(double) *outtbcptr;
4214                         cond = (*outtbcptr <= 0);
4215                         break;
4216
4217                 case SCSINT16_N   : outtbsptr=(SCSINT16_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4218                         g[zcptr[kf-1]-1]=(double) *outtbsptr;
4219                         cond = (*outtbsptr <= 0);
4220                         break;
4221
4222                 case SCSINT32_N   : outtblptr=(SCSINT32_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4223                         g[zcptr[kf-1]-1]=(double) *outtblptr;
4224                         cond = (*outtblptr <= 0);
4225                         break;
4226
4227                 case SCSUINT8_N   : outtbucptr=(SCSUINT8_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4228                         g[zcptr[kf-1]-1]=(double) *outtbucptr;
4229                         cond = (*outtbucptr <= 0);
4230                         break;
4231
4232                 case SCSUINT16_N  : outtbusptr=(SCSUINT16_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4233                         g[zcptr[kf-1]-1]=(double) *outtbusptr;
4234                         cond = (*outtbusptr <= 0);
4235                         break;
4236
4237                 case SCSUINT32_N  : outtbulptr=(SCSUINT32_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4238                         g[zcptr[kf-1]-1]=(double) *outtbulptr;
4239                         cond = (*outtbulptr <= 0);
4240                         break;
4241
4242                 default  : /* Add a message here */
4243                         *ierr = 25;
4244                         return 0;
4245                         break;
4246                 }
4247                 if (cond) {
4248                         i=2;
4249                 }
4250                 else {
4251                         i=1;
4252                 }
4253         }
4254         /* eselect blk */
4255         else if (funtyp[kf-1] == -2) {
4256                 switch(outtbtyp[-1+inplnk[inpptr[kf-1]-1]])
4257                 {
4258                 case SCSREAL_N    : outtbdptr=(SCSREAL_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4259                         for (jj=0;jj<scs_imp->blocks[kf-1].nevout-1;jj++) {
4260                                 g[zcptr[kf-1]-1+jj]=*outtbdptr-(double)(jj+2);
4261                         }
4262                         i=max(min((int) *outtbdptr,scs_imp->blocks[kf-1].nevout),1);
4263                         break;
4264
4265                 case SCSCOMPLEX_N : outtbdptr=(SCSCOMPLEX_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4266                         for (jj=0;jj<scs_imp->blocks[kf-1].nevout-1;jj++) {
4267                                 g[zcptr[kf-1]-1+jj]=*outtbdptr-(double)(jj+2);
4268                         }
4269                         i=max(min((int) *outtbdptr,scs_imp->blocks[kf-1].nevout),1);
4270                         break;
4271
4272                 case SCSINT8_N    : outtbcptr=(SCSINT8_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4273                         for (jj=0;jj<scs_imp->blocks[kf-1].nevout-1;jj++) {
4274                                 g[zcptr[kf-1]-1+jj]=(double) *outtbcptr-(double)(jj+2);
4275                         }
4276                         i=max(min((int) *outtbcptr,scs_imp->blocks[kf-1].nevout),1);
4277                         break;
4278
4279                 case SCSINT16_N   : outtbsptr=(SCSINT16_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4280                         for (jj=0;jj<scs_imp->blocks[kf-1].nevout-1;jj++) {
4281                                 g[zcptr[kf-1]-1+jj]=(double) *outtbsptr-(double)(jj+2);
4282                         }
4283                         i=max(min((int) *outtbsptr,scs_imp->blocks[kf-1].nevout),1);
4284                         break;
4285
4286                 case SCSINT32_N   : outtblptr=(SCSINT32_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4287                         for (jj=0;jj<scs_imp->blocks[kf-1].nevout-1;jj++) {
4288                                 g[zcptr[kf-1]-1+jj]=(double) *outtblptr-(double)(jj+2);
4289                         }
4290                         i=max(min((int) *outtblptr,scs_imp->blocks[kf-1].nevout),1);
4291                         break;
4292
4293                 case SCSUINT8_N   : outtbucptr=(SCSUINT8_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4294                         for (jj=0;jj<scs_imp->blocks[kf-1].nevout-1;jj++) {
4295                                 g[zcptr[kf-1]-1+jj]=(double) *outtbucptr-(double)(jj+2);
4296                         }
4297                         i=max(min((int) *outtbucptr,scs_imp->blocks[kf-1].nevout),1);
4298                         break;
4299
4300                 case SCSUINT16_N  : outtbusptr=(SCSUINT16_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4301                         for (jj=0;jj<scs_imp->blocks[kf-1].nevout-1;jj++) {
4302                                 g[zcptr[kf-1]-1+jj]=(double) *outtbusptr-(double)(jj+2);
4303                         }
4304                         i=max(min((int) *outtbusptr,scs_imp->blocks[kf-1].nevout),1);
4305                         break;
4306
4307                 case SCSUINT32_N  : outtbulptr=(SCSUINT32_COP *)outtbptr[-1+inplnk[inpptr[kf-1]-1]];
4308                         for (jj=0;jj<scs_imp->blocks[kf-1].nevout-1;jj++) {
4309                                 g[zcptr[kf-1]-1+jj]=(double) *outtbulptr-(double)(jj+2);
4310                         }
4311                         i=max(min((int) *outtbulptr,scs_imp->blocks[kf-1].nevout),1);
4312                         break;
4313
4314                 default  : /* Add a message here */
4315                         *ierr = 25;
4316                         return 0;
4317                         break;
4318                 }
4319         }
4320         return i;
4321 } /* synchro_g_nev */
4322 /*--------------------------------------------------------------------------*/
4323 /* FREE_blocks */
4324 static void FREE_blocks()
4325 {
4326         int kf = 0;
4327         for (kf = 0; kf < nblk; ++kf) {
4328                 if (Blocks[kf].insz!=NULL) {
4329                         FREE(Blocks[kf].insz);
4330                 }else {
4331                         break;
4332                 }
4333                 if (Blocks[kf].inptr!=NULL){
4334                         FREE(Blocks[kf].inptr);
4