11692ee69e8fb495544b4f6023cfcb0904831cf3
[scilab.git] / scilab / modules / graphics / src / c / Contour.c
1 /*------------------------------------------------------------------------
2  *    Graphic library
3  *    Copyright (C) 1998-2001 Enpc/Jean-Philippe Chancelier
4  *    jpc@cermics.enpc.fr 
5  --------------------------------------------------------------------------*/
6
7 #include <string.h> /* in case of dbmalloc use */
8 #include <stdio.h>
9 #include <math.h>
10 #include "math_graphics.h"
11 #include "Graphics.h" 
12 #include "PloEch.h"
13 #include "Xcall1.h"
14 #include "sciprint.h"
15 #include "MALLOC.h"
16
17 #include "GetProperty.h"
18
19 /* #include "Entities.h" /\* F.Leray 21.04.04 : for update_2dbounds call*\/ */
20 /*extern void compute_data_bounds(int cflag,char dataflag,double *x,double *y,int n1,int n2,double *drect);*/
21 extern void compute_data_bounds2(int cflag,char dataflag,char *logflags,double *x,double *y,int n1,int n2,double *drect);
22 extern BOOL update_specification_bounds(sciPointObj *psubwin, double *rect,int flag);
23 extern int re_index_brect(double * brect, double * drect);
24 extern BOOL strflag2axes_properties(sciPointObj * psubwin, char * strflag);
25
26 #ifdef _MSC_VER
27 extern void Scistring (char *str);
28 #endif
29
30 typedef void (level_f) __PARAMS((integer ival, double Cont, double xncont,
31                                double yncont));
32 typedef void (*ptr_level_f) __PARAMS((integer ival, double Cont, double xncont,
33                                double yncont));
34
35 static int 
36 Contour2D __PARAMS((ptr_level_f,char *,double *x,double *y,double *z,integer *n1,
37                    integer *n2,integer *flagnz,integer *nz,double *zz,
38                    integer *style,char *strflag,char *legend,double *brect,
39                    integer *aaint,integer lstr1,integer lstr2));
40 static void 
41 contourI __PARAMS((ptr_level_f,double *, double *, double *,
42                   double *, integer *, integer *, integer *));
43 static void
44 look __PARAMS((ptr_level_f, integer i, integer j, integer ib,
45               integer jb, integer qq,double Cont, integer style));
46
47 static integer ffnd __PARAMS((ptr_level_f,integer,integer,integer,integer,integer,
48                              integer,integer,integer,integer,integer,
49                              double,integer *));
50
51 static int Gcont_size = 0;
52
53 static void ContourTrace __PARAMS((double Cont, integer style));
54 static level_f ContStore, ContStore1, ContStore2,GContStore2;
55 static void GContStore2Last __PARAMS((void));
56 static double x_cont __PARAMS((integer i));
57 static double y_cont __PARAMS((integer i));
58 static double phi_cont __PARAMS((integer, integer)); 
59 static double f_intercept  __PARAMS((double, double, double, double, double ));
60 static integer not_same_sign  __PARAMS((double val1, double val2)); 
61 static integer get_itg_cont  __PARAMS((integer i, integer j)); 
62 static void inc_itg_cont  __PARAMS((integer i, integer j, integer val)); 
63 static integer oddp  __PARAMS((integer i));
64
65 /*-----------------------------------------------------------------------
66  *  Level curves 
67  *  The computer journal vol 15 nul 4 p 382 (1972)
68  *  from the Lisp Macsyma source (M.I.T)
69  * -------------------------------------------------------------------------*/
70
71 #define HIDDENFRAMECOLOR 2L /* default color for hidden frame */
72
73 /*---------------------------------------------------------------------------
74  * General functions (could be changed in #define or
75  *   inline functions to increase speed)
76  *---------------------------------------------------------------------------*/
77
78 static double *GX,*GY,*GZ;
79 static integer Gn1,Gn2;
80
81 static double * Gxcont = NULL ;
82 static double * Gycont = NULL ;
83
84 static integer * xcont = NULL ;
85 static integer * ycont = NULL ;
86
87 static void InitValues(double *x, double *y, double *z, integer n1, integer n2)
88 {
89   Gn1=n1;  Gn2=n2;  GX = x;  GY = y;  GZ = z;
90 }
91
92 /*--------return the  value of f for a pointeger on the grid-----*/
93
94 static double phi_cont(integer i, integer j)
95 {
96   return(GZ[i+Gn1*j]);
97 }
98
99 /*---------return the coordinates between  [xi,xj] along one axis 
100  *  for which the value of f is zCont */ 
101
102 static double f_intercept(double zCont, double fi, double xi, double fj, double xj)
103 {
104   return( xi+ (zCont-fi)*(xj-xi)/ (fj-fi));
105 }
106
107 /* check for boundary points */
108
109 static integer  bdyp(integer i, integer j)
110 {
111   return (  j == 0 || i == 0 || j == Gn2-1 || i == Gn1-1 );
112 }
113
114 /* store or get flag values */
115
116 static  integer *itg_cont, *xbd_cont,*ybd_cont;
117
118 static integer get_itg_cont(integer i, integer j)
119 {
120   return( itg_cont[i+Gn1*j]);
121 }
122
123 static void inc_itg_cont(integer i, integer j, integer val)
124 {
125   itg_cont[i+Gn1*j] += val;
126 }
127
128 static integer not_same_sign(double val1, double val2)
129 {
130   if ( ISNAN(val1) ==1 || ISNAN(val2) == 1) return(0);
131   /** 0.0 est consid\'er\'e comme positif **/
132   if ( val1 >= 0.0) 
133     {
134       if (val2 < 0.0) return(1) ; else return(0);}
135   else 
136     {
137       if ( val2 >= 0.0) return(1) ; else return(0);}
138 }
139
140 static integer oddp(integer i) { return( i == 1 || i ==3 );}
141
142 /*---------return the x-value of a grid point--------*/
143
144 static double x_cont(integer i) { return GX[i] ;}
145
146 /*---------return the y-value of a grid point --------*/
147
148 static double y_cont(integer i) { return GY[i] ;}
149
150 /*------------------------------------------------------------
151  * Draw level curves for a function f(x,y) which values 
152  * at points x(i),y(j) are given by z(i,j)
153  * - z is a (n1,n2) matrix 
154  * - x is a (1,n1) matrix 
155  * - y is a (1,n2) matrix 
156  * - x,y,z are stored as one dimensionnal array in C 
157  * - if *flagnz =0 
158  * -   then  nz is an integer pointer to the number of level curves. 
159  *     else  zz is an array which gives th requested level values.
160  *            (and nz is the size of thos array)
161  * Computed from min and max of z
162  * Exemple Contour(1:5,1:10,rand(5,10),5);
163  *---------------------------------------------------------------*/
164
165 static double ZC=0.0;
166 static char   ContNumFormat[100];
167
168 /** interface for contour2d **/
169
170 int C2F(contour2)(double *x, double *y, double *z, integer *n1, integer *n2, integer *flagnz, integer *nz, double *zz, integer *style, char *strflag, char *legend, double *brect, integer *aaint, integer lstr1, integer lstr2)
171 {
172   Contour2D(ContStore2,"contour2",x,y,z,n1,n2,flagnz,nz,zz,style,strflag,
173             legend,brect,aaint,lstr1,lstr2);
174   return 0;
175 }
176
177 /* interface for contour2di used in macro contourf 
178  * used when we want to get the values which constitues the contour inside Scilab 
179  * contour2di + c2dex 
180  * THIS PROCEDURE IS NO LONGUER USED 
181  */
182
183 static int Contour2D(ptr_level_f func, char *name, double *x, double *y, double *z, integer *n1, integer *n2, integer *flagnz, integer *nz, double *zz, integer *style, char *strflag, char *legend, double *brect, integer *aaint, integer lstr1, integer lstr2)
184 {
185   integer err=0;
186   static double *zconst;
187   double zmin,zmax;
188   integer N[3],i;
189   double drect[6];
190   sciPointObj * psubwin = NULL;  /* Adding F.Leray 22.04.04 */
191   BOOL bounds_changed = FALSE;
192   BOOL axes_properties_changed = FALSE;
193   
194   /** Boundaries of the frame **/
195   psubwin = sciGetSelectedSubWin (sciGetCurrentFigure ());
196
197   /* Force psubwin->is3d to FALSE: we are in 2D mode */
198   if (sciGetSurface(psubwin) == (sciPointObj *) NULL)
199   {
200     pSUBWIN_FEATURE (psubwin)->is3d = FALSE;
201     pSUBWIN_FEATURE (psubwin)->project[2]= 0;
202   }
203   else
204   {
205     pSUBWIN_FEATURE (psubwin)->theta_kp=pSUBWIN_FEATURE (psubwin)->theta;
206     pSUBWIN_FEATURE (psubwin)->alpha_kp=pSUBWIN_FEATURE (psubwin)->alpha;  
207   }
208
209   pSUBWIN_FEATURE (psubwin)->alpha  = 0.0;
210   pSUBWIN_FEATURE (psubwin)->theta  = 270.0;
211
212   /*****TO CHANGE F.Leray 10.09.04    for (i=0;i<4;i++) pSUBWIN_FEATURE(psubwin)->axes.aaint[i] = aaint[i]; */
213   if (sciGetGraphicMode (psubwin)->autoscaling) {
214     /* compute and merge new specified bounds with psubwin->Srect */
215     switch (strflag[1])  {
216         case '0': 
217           /* do not change psubwin->Srect */
218           break;
219         case '1' : case '3' : case '5' : case '7':
220           /* Force psubwin->Srect=brect */
221           re_index_brect(brect,drect);
222           break;
223         case '2' : case '4' : case '6' : case '8':
224           /* Force psubwin->Srect to the x and y bounds */
225           /*  compute_data_bounds(1,'g',x,y,*n1,*n2,drect); */
226           compute_data_bounds2(1,'g',pSUBWIN_FEATURE(psubwin)->logflags,x,y,*n1,*n2,drect);
227           break;
228     }
229     if (!pSUBWIN_FEATURE(psubwin)->FirstPlot &&(strflag[1] == '7' || strflag[1] == '8')) { /* merge psubwin->Srect and drect */
230       drect[0] = Min(pSUBWIN_FEATURE(psubwin)->SRect[0],drect[0]); /*xmin*/
231       drect[2] = Min(pSUBWIN_FEATURE(psubwin)->SRect[2],drect[2]); /*ymin*/
232       drect[1] = Max(pSUBWIN_FEATURE(psubwin)->SRect[1],drect[1]); /*xmax*/
233       drect[3] = Max(pSUBWIN_FEATURE(psubwin)->SRect[3],drect[3]); /*ymax*/
234     }
235     if (strflag[1] != '0')
236       bounds_changed = update_specification_bounds(psubwin, drect,2);
237   } 
238
239   if(pSUBWIN_FEATURE (psubwin)->FirstPlot == TRUE) bounds_changed = TRUE;
240
241   axes_properties_changed = strflag2axes_properties(psubwin, strflag);
242
243   pSUBWIN_FEATURE (psubwin)->FirstPlot = FALSE; /* just after strflag2axes_properties */
244
245   zmin=(double) Mini(z,*n1*(*n2)); 
246   zmax=(double) Maxi(z,*n1*(*n2));
247
248   /** Scales **/
249   if (strcmp(name,"contour2")==0 )
250     {
251       axis_draw(strflag);
252       frame_clip_on();
253     }
254   if (*flagnz==0)
255     {
256       if (( zconst = MALLOC( (*nz) * sizeof(double) ) )== 0) 
257         {
258           sciprint("Running out of memory\r\n");
259           return 0;
260         }
261       for ( i =0 ; i < *nz ; i++) 
262         zconst[i]=zmin + (i+1)*(zmax-zmin)/(*nz+1);
263       N[0]= *n1;N[1]= *n2;N[2]= *nz;
264       contourI(func,x,y,z,zconst,N,style,&err);
265       FREE(zconst) ;
266       zconst = NULL ;
267     }
268   else
269     {
270       N[0]= *n1;N[1]= *n2;N[2]= *nz;
271       contourI(func,x,y,z,zz,N,style,&err);
272     }
273
274   if (strcmp(name,"contour2")==0 )frame_clip_off();
275
276   return(0);
277 }
278
279 int C2F(contourif)(double *x, double *y, double *z, integer *n1, integer *n2, integer *flagnz, integer *nz, double *zz, integer *style)
280 {
281   integer err=0;
282   static double *zconst;
283   double zmin,zmax;
284   integer N[3],i;
285
286   zmin=(double) Mini(z,*n1*(*n2)); 
287   zmax=(double) Maxi(z,*n1*(*n2));
288
289   if (*flagnz==0)
290     {
291       if ( ( zconst = MALLOC( (*nz) * sizeof(double) ) ) == 0 ) 
292         {
293           sciprint("Running out of memory\r\n");
294           return 0;
295         }
296       for ( i =0 ; i < *nz ; i++) 
297         zconst[i]=zmin + (i+1)*(zmax-zmin)/(*nz+1);
298       N[0]= *n1;N[1]= *n2;N[2]= *nz;
299       contourI(GContStore2,x,y,z,zconst,N,style,&err);
300       FREE(zconst) ;
301       zconst = NULL ;
302     }
303   else
304     {
305       N[0]= *n1;N[1]= *n2;N[2]= *nz;
306       contourI(GContStore2,x,y,z,zz,N,style,&err);
307     }
308
309   return(0);
310 }
311
312 /*-------------------------------------------------------
313  *  The function f is given on a grid and we want the level curves 
314  *  for the zCont[N[2]] values 
315  *  x : of size N[0] gives the x-values of the grid 
316  *  y : of size N[1] gives the y-values of the grid 
317  *  z : of size N[0]*N[1]  gives the f-values on the grid 
318  *  style: size ncont (=N[2]) or empty integer pointer 
319  *  gives the dash style for contour i
320  *-------------------------------------------------------*/
321
322 static void contourI(ptr_level_f func, double *x, double *y, double *z, double *zCont, integer *N, integer *style, integer *err)
323 {
324   int check = 1;
325   char *F;
326   integer n1,n2,ncont,i,c,j,k,n5;
327   integer stylec;
328   n1=N[0];n2=N[1];ncont=N[2];
329   F=getFPF();
330   if ( F[0] == '\0') 
331     ChoixFormatE1(ContNumFormat,zCont,N[2]);
332   InitValues(x,y,z,n1,n2);
333   n5 =  2*(n1)+2*(n2)-3;
334   /* Allocation */
335   Gcont_size = 0; /** initialize the array indices for storing contours **/
336   xbd_cont = MALLOC( n5 * sizeof(int) ) ;
337   ybd_cont = MALLOC( n5 * sizeof(int) ) ;
338   itg_cont = MALLOC( n1*n2 * sizeof(int) ) ;
339   if ( (xbd_cont == NULL) && n5 != 0) check= 0;
340   if ( (ybd_cont == NULL) && n5 != 0) check= 0;
341   if ( (itg_cont == NULL) && n1*n2 != 0) check= 0;
342   if ( check == 0) 
343     {
344       FREE( xbd_cont ) ;
345       FREE( ybd_cont ) ;
346       FREE( itg_cont ) ;
347       sciprint("contourI_: Running out of memory\n");
348       return;
349     }
350   /* just a parametrization of the boundary points */
351   for ( i = 0 ; i <  n2 ; i++)
352         {
353           ybd_cont[i] = i ;
354           xbd_cont[i] = 0 ;
355         }
356   for ( i = 1 ; i <  n1 ; i++)
357         {
358           ybd_cont[n2+i-1] = n2-1 ;
359           xbd_cont[n2+i-1] = i  ;
360         }
361   for ( i = n2-2;  i >= 0  ; i--)
362         {
363           ybd_cont[2*n2 +n1-3-i] = i ;
364           xbd_cont[2*n2 +n1-3-i] = n1-1  ;
365         }
366   for ( i = n1-2 ; i >= 0 ; i--)
367         {
368           ybd_cont[2*n2 +2*n1-4-i] = 0 ;
369           xbd_cont[2*n2 +2*n1-4-i] = i  ;
370         }
371   for ( c= 0 ; c < ncont ; c++)
372     {
373       stylec = ( style != (integer *) 0) ? style[c] : c;
374       /** itg-cont is a flag array to memorize checked parts of the grid **/
375       for ( i = 0 ; i < n1; i++)
376         for ( j =0 ; j < n2 ; j++)
377           itg_cont[i+n1*j]=0 ;
378       /** all the boundary segments **/
379       for ( k = 1 ; k < n5 ; k++)
380         { integer ib,jb;
381         i = xbd_cont[k] ; j = ybd_cont[k];
382         ib = xbd_cont[k-1] ; jb= ybd_cont[k-1];
383         if  (not_same_sign (phi_cont(i,j)-zCont[c] , 
384                             phi_cont(ib,jb)-zCont[c]))
385           look(func,i,j,ib,jb,1L,zCont[c],stylec);
386         }
387       /** inside segments **/
388       for ( i = 1 ; i < n1-1; i++)
389         for ( j = 1 ; j < n2-1 ; j++)
390           if  (not_same_sign ( phi_cont(i,j)-zCont[c] , 
391                                phi_cont(i, j-1)-zCont[c]))
392             look(func,i,j,i,j-1,2L,zCont[c],stylec);
393     }
394   FREE( xbd_cont ) ;
395   FREE( ybd_cont ) ;
396   FREE( itg_cont ) ;
397
398 }
399
400 /*--------------------------------------------------------------------
401  *  the level curve is crossing the segment (i,j) (ib,jb)
402  *  look store the level curve point and try to find the next segment to look at
403  *  Cont: value of f along the contour 
404  *  ncont: number of contour 
405  *  c: indice of the contour Cont 
406  *---------------------------------------------------------------------*/
407
408 static void look(ptr_level_f func, integer i, integer j, integer ib, integer jb, integer qq, double Cont, integer style)
409 {
410   integer ip,jp,im,jm,zds,ent=0,flag=0,wflag;
411   jp= j+1; ip= i+1; jm=j-1;im=i-1;
412   /*  on regarde comment est le segment de depart */
413   if  ( jb == jm)  flag = 1; 
414   else  { 
415     if ( ib == im ) flag = 2 ;
416     else  {
417       if ( jb == jp ) flag = 3 ;
418       else  if ( ib == ip ) flag = 4;}}
419   switch  (  flag)
420     {
421     case  1 :
422       if  (get_itg_cont(i,jm) > 1) return;
423       ent=1 ; /* le segment est vertical vers le bas */
424       /* Storing intersection point */
425       (*func)(0,Cont, x_cont(i), 
426                 f_intercept(Cont,phi_cont(i,jm),
427                             y_cont(jm),phi_cont(i,j),y_cont(j)));
428       break;
429     case 2 : 
430       if  (get_itg_cont(im,j) == 1 || get_itg_cont(im,j)==3 ) return;
431       ent=2 ; /* le segment est horizontal gauche */
432       /* Storing intersection point */
433       (*func)( 0,Cont,
434                 f_intercept(Cont,phi_cont(im,j),
435                             x_cont(im),phi_cont(i,j),x_cont(i)), y_cont(j));
436       break ; 
437     case 3 :
438       if  (get_itg_cont(i,j) > 1 ) return;
439       ent=3 ; /* le segment est vertical haut */
440       /* Storing intersection point */
441       (*func)(0,Cont,x_cont(i), f_intercept(Cont,phi_cont(i,j),
442                                             y_cont(j),phi_cont(i,jp),y_cont(jp)));
443       break ;
444     case 4 :
445       if  (get_itg_cont(i,j) == 1 || get_itg_cont(i,j)==3 ) return;
446       ent=4 ; /* le segment est horizontal droit */
447       /* Storing intersection point */
448       (*func)(0,Cont,f_intercept(Cont,phi_cont(i,j),
449                                  x_cont(i),phi_cont(ip,j),x_cont(ip)),
450               y_cont(j));
451       break;
452       default :
453         sciprint(" Error in case wrong value ");
454       break;
455     }
456   wflag=1;
457   while ( wflag) 
458     { 
459       jp= j+1; ip= i+1; jm=j-1;im=i-1;
460       switch  ( ent) 
461         {case 1 :
462           inc_itg_cont(i,jm,2L);
463           ent = ffnd(func,i,ip,ip,i,j,j,jm,jm,ent,qq,Cont,&zds);
464           /* on calcule le nouveau point, ent donne la 
465              direction du segment a explorer */
466           switch ( ent)
467             {
468             case -1: wflag=0; break;
469             case 1 : i=ip ; break ;
470             case 2 : i=ip;j=jm; break ;
471             }
472           break ;
473         case 2  :
474           inc_itg_cont(im,j,1L);
475           ent = ffnd(func,i,i,im,im,j,jm,jm,j,ent,qq,Cont,&zds);
476           switch ( ent)
477             { 
478             case -1: wflag=0; break;
479             case 2 : j = jm ;break ;
480             case  3  : i=im;j=jm; break ;
481             }
482           break ;
483         case 3 :
484           inc_itg_cont(i,j,2L);
485           ent = ffnd(func,i,im,im,i,j,j,jp,jp,ent,qq,Cont,&zds);
486           switch ( ent)
487             { 
488             case -1: wflag=0; break;
489             case 3 : i=im; break ;
490             case 4 : i=im;j=jp; break ;
491             }
492           break ;
493         case 4 :
494           inc_itg_cont(i,j,1L);
495           ent = ffnd(func,i,i,ip,ip,j,jp,jp,j,ent,qq,Cont,&zds);
496           switch ( ent)
497             {
498             case -1: wflag=0; break;
499             case 4 :j=jp;break ;
500             case 1 :i=ip;j=jp;break ;
501             }
502           break ;
503         }
504      
505       /** new segment is on the boundary **/
506       if ( zds == 1) 
507         {
508           switch ( ent) 
509             {
510             case 1 : inc_itg_cont(i,(j-1),2L); break ; 
511             case 2 : inc_itg_cont(i-1,j,1L);  break ; 
512             case 3 : inc_itg_cont(i,j,2L); break ; 
513             case 4 : inc_itg_cont(i,j,1L); break ; 
514             }
515           /** we must quit the while loop **/
516           wflag = 0 ;
517           }
518       /**  init point was inside the domain */
519       if ( qq == 2) 
520         {
521           switch ( ent) 
522             {
523             case 1 : if  ( get_itg_cont (i,j-1)  > 1) wflag = 0 ; break ; 
524             case 2 : if  ( oddp(get_itg_cont(i-1,j))) wflag = 0 ; break ; 
525             case 3 : if  ( get_itg_cont(i,j) > 1)     wflag = 0 ; break ; 
526             case 4 : if  ( oddp(get_itg_cont(i,j)))   wflag = 0 ; break ; 
527             }
528         }
529     }
530   if ( func == GContStore2 ) 
531     GContStore2Last();
532   else 
533     ContourTrace(Cont,style);
534 }
535
536
537 /*-----------------------------------------------------------------------
538  *   ffnd : cette fonction  recoit en entree quatre points 
539  *       on sait que la courbe de niveau passe entre le point 1 et le quatre 
540  *       on cherche a savoir ou elle resort, 
541  *       et on fixe une nouvelle valeur de ent qui indiquera le segment 
542  *       suivant a explorer 
543  *-----------------------------------------------------------------------*/
544
545 static integer ffnd (ptr_level_f func, integer i1, integer i2, integer i3, integer i4, integer jj1, integer jj2, integer jj3, integer jj4, integer ent, integer qq, double Cont, integer *zds)
546 {
547   double phi1,phi2,phi3,phi4,xav,yav,phiav;
548   integer revflag,i;
549   phi1=phi_cont(i1,jj1)-Cont;
550   phi2=phi_cont(i2,jj2)-Cont;
551   phi3=phi_cont(i3,jj3)-Cont;
552   phi4=phi_cont(i4,jj4)-Cont;
553   revflag = 0;
554   *zds = 0;
555   /* le point au centre du rectangle */
556   xav = ( x_cont(i1)+ x_cont(i3))/2.0 ; 
557   yav = ( y_cont(jj1)+ y_cont(jj3))/2.0 ; 
558   phiav = ( phi1+phi2+phi3+phi4) / 4.0;
559   if (ISNAN(phiav)==1) 
560     {
561       return -1;
562     }
563   if (  not_same_sign( phiav,phi4)) 
564     {
565       integer l1, k1; 
566       double phi;
567       revflag = 1 ; 
568       l1= i4; k1= jj4;
569       i4=i1; jj4 = jj1; i1= l1; jj1= k1;
570       l1= i3; k1= jj3;
571       i3=i2; jj3= jj2; i2=l1; jj2= k1;
572       phi = phi1; phi1 = phi4; phi4= phi;
573       phi = phi3; phi3 = phi2; phi2= phi;
574     }
575   /* on stocke un nouveau point  */
576   (*func)(1,Cont,f_intercept(0.0,phi1,x_cont(i1),phiav,xav),
577             f_intercept(0.0,phi1,y_cont(jj1),phiav,yav));
578   /*
579    * on parcourt les segments du rectangle pour voir sur quelle face
580    * on sort 
581    */
582   for  ( i = 0 ;  ; i++)
583     { integer l1,k1;
584       double phi;
585       if ( not_same_sign ( phi1,phi2))   /** sortir du for **/ break ; 
586       if  ( phiav != 0.0 ) 
587         {
588           (*func)(1,Cont,f_intercept(0.0,phi2,x_cont(i2),phiav,xav),
589                     f_intercept(0.0,phi2,y_cont(jj2),phiav,yav));
590         } 
591       /** on permutte les points du rectangle **/
592       l1=i1; k1= jj1;
593       i1=i2;jj1=jj2;i2=i3;jj2=jj3;i3=i4;jj3=jj4;i4=l1;jj4=k1;
594       phi=phi1; phi1=phi2;phi2=phi3;phi3=phi4;phi4=phi;
595     }
596   (*func)(1,Cont,f_intercept(0.0,phi1,x_cont(i1),phi2,x_cont(i2)),
597             f_intercept(0.0,phi1,y_cont(jj1),phi2,y_cont(jj2)));
598   if ( qq==1 && bdyp(i1,jj1) && bdyp(i2,jj2)) *zds = 1 ;
599   if ( revflag == 1  &&  ! oddp (i) )  i = i+2;
600   return ( 1 + (  ( i + ent + 2) % 4));
601 }
602
603 /*--------------------------------------------------------------
604  * Storing and tracing level curves 
605  *----------------------------------------------------------------*/
606
607 static integer cont_size ;
608
609 /*
610  * store a point in the current level curve if ival == 0 the level 
611  * curve is reinitialized 
612  * used for a contour in a 3D drawing 
613  */
614
615 static void
616 G_ContStore(integer ival, int xncont, int yncont)
617 {
618   int n;
619   /* nouveau contour */
620   if ( ival == 0) cont_size =0 ;
621   n= cont_size + 1;
622   xcont = REALLOC( xcont, n * sizeof(int) ) ;
623   ycont = REALLOC( ycont, n * sizeof(int) ) ;
624   if ( ( xcont == NULL || ycont == NULL ) && n != 0 )
625   {
626     FREE( xcont ) ;
627     FREE( ycont ) ;
628
629     return ;
630   }
631   xcont[cont_size]= xncont;
632   ycont[cont_size++]= yncont;
633 }
634
635 /*
636  * store a point in the current level curve if ival == 0 the level 
637  * curve is reinitialized 
638  * used for a contour in a 3D drawing 
639  */
640
641 static void
642 ContStore(integer ival, double Cont, double xncont, double yncont)
643 {
644   G_ContStore(ival,GEOX(xncont,yncont,Cont),
645               GEOY(xncont,yncont,Cont));
646 }
647
648 /*
649  * store a point in the current level curve if ival == 0 the level 
650  * curve is reinitialized 
651  * used for a contour in a 3D drawing with projection at level ZC 
652  */
653
654 static void
655 ContStore1(integer ival, double Cont, double xncont, double yncont)
656 {
657   G_ContStore(ival,GEOX(xncont,yncont,ZC),
658               GEOY(xncont,yncont,ZC));
659 }
660
661 /*
662  * store a point in the current level curve if ival == 0 the level 
663  * curve is reinitialized 
664  * used for a contour in a 2D drawing 
665  */
666
667 static void
668 ContStore2(integer ival, double Cont, double xncont, double yncont)
669 {
670   G_ContStore(ival,XScale(xncont),YScale(yncont));
671 }
672
673 /* 
674  * Explicit drawing of the current level curve with a dash style 
675  * The curve level is also drawn as a string according to current 
676  * floating point format 
677  */
678
679 static void ContourTrace(double Cont, integer style)
680
681   char *F;
682   integer verbose=0 ,Dnarg,Dvalue[10];
683   integer close=0, flag=0, uc;
684   double angle=0.0;
685   char str[100];
686
687   C2F(dr)("xget","use color",&verbose,&uc,&Dnarg,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
688   if (uc) {
689     C2F(dr)("xget","color",&verbose,Dvalue,&Dnarg,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
690     C2F(dr)("xset","color",&style,PI0,PI0,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
691     C2F(dr)("xlines","void",&cont_size,xcont,ycont,&close,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
692     C2F(dr)("xset","color",Dvalue,PI0,PI0,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
693   }
694   else {
695     C2F(dr)("xget","line style",&verbose,Dvalue,&Dnarg,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
696     C2F(dr)("xset","line style",&style,PI0,PI0,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
697     C2F(dr)("xlines","void",&cont_size,xcont,ycont,&close,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
698     C2F(dr)("xset","line style",Dvalue,PI0,PI0,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
699   }
700
701   F=getFPF();
702   if ( F[0] == '\0') 
703     sprintf(str,ContNumFormat,Cont);
704   else 
705     sprintf(str,F,Cont);
706   C2F(dr)("xstring",str, &xcont[cont_size / 2],&ycont[cont_size /2],
707           PI0,&flag,PI0,PI0, &angle,PD0,PD0,PD0,0L,0L);
708 }
709
710 /*--------------------------------------------------------------
711  * Following code is used to store the current level curves as 
712  * double in order to access to the stored data at Scilab level 
713  *----------------------------------------------------------------*/
714
715 static int last=-1;
716 static int count=0; 
717  
718 /** used to bring back data to Scilab Stack **/
719
720 int C2F(getconts)(double **x, double **y, integer *mm, integer *n)
721 {
722   *x = Gxcont;
723   *y = Gycont;
724   *mm= 1;
725   *n= Gcont_size;
726   return 0;
727 }
728
729 static void GContStore2(integer ival, double Cont, double xncont, double yncont)
730 {
731   int n;
732   if ( ival == 0) 
733     {
734       /* Here : ival == 0 means stop the current level curve and 
735        * store data at the end but do reset Gcont_size to zero 
736        */
737       n= Gcont_size + 2;
738       if ( Gxcont == NULL )
739       {
740         Gxcont = MALLOC( n * sizeof(double) ) ;
741         Gycont = MALLOC( n * sizeof(double) ) ;
742       }
743       else
744       {
745         Gxcont = REALLOC( Gxcont, n * sizeof(double) ) ;
746         Gycont = REALLOC( Gycont, n * sizeof(double) ) ;
747       }
748       if ( (Gxcont == NULL) && n != 0) return ; 
749       if ( (Gycont == NULL) && n != 0) return ;
750       Gxcont[Gcont_size] = Cont;
751       if ( last != -1 && last < n ) Gycont[last]= count;
752       last = Gcont_size;
753       Gcont_size++;
754       count = 0;
755     }
756   else
757   {
758     n = Gcont_size + 1 ;
759     if ( Gxcont == NULL )
760     {
761       Gxcont = MALLOC( n * sizeof(double) ) ;
762       Gycont = MALLOC( n * sizeof(double) ) ;
763     }
764     else
765     {
766       Gxcont = REALLOC( Gxcont, n * sizeof(double) ) ;
767       Gycont = REALLOC( Gycont, n * sizeof(double) ) ;
768     }
769     if ( (Gxcont == NULL) && n != 0) return ; 
770     if ( (Gycont == NULL) && n != 0) return ;
771   }
772   Gxcont[Gcont_size]=xncont;
773   Gycont[Gcont_size++]=yncont;
774   count++;
775 }
776
777 static void GContStore2Last(void)
778 {
779   if ( last != -1 ) Gycont[last]= count;
780 }
781
782
783
784
785
786
787
788
789
790
791