Partiel fix of bug 4658
[scilab.git] / scilab / modules / graphics / macros / contourf.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA
3 // This file must be used under the terms of the CeCILL.
4 // This source file is licensed as described in the file COPYING, which
5 // you should have received as part of this distribution.  The terms
6 // are also available at    
7 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
8
9 function contourf(x,y,z,nv,style,strf,leg,rect,nax)
10         
11         [nout,nin]=argn(0);
12         
13         if nin == 0 then   // demo
14                 t = -%pi:0.1:%pi;
15                 m = sin(t)'*cos(t);
16                 contourf(t,t,m);
17                 return;
18   end
19
20 if nin <= 0 then x=1:10;end
21 if nin <= 1 then y=1:10;end
22 if nin <= 2 then z=rand(size(x,'*'),size(y,'*'));end
23 if nin <= 3 then zmin=mini(z);zmax=maxi(z);nv = zmin + (1:10)*(zmax-zmin)/(11);end
24 if nin <= 5 then strf="121";end
25 if nin <= 6 then leg=" ";end
26 if nin <= 7 then rect=[0,0,1,1];end
27 if nin <= 8 then nax=[1,10,1,10];end
28 if x==[] then x=1:size(z,'r');end 
29 if y==[] then y=1:size(z,'c');end 
30
31 nvs=size(nv,'*') ;
32 if nvs==1 then nvs=nv;zmin=mini(z);zmax=maxi(z);nv = zmin + (1:nvs)*(zmax-zmin)/(nvs+1);end;
33 if nin <= 4 then style = -1*ones(1,nvs);end
34 if nin <= 7 then rect=[mini(x),mini(y),maxi(x),maxi(y)]; end 
35 nv1=nv
36 [mz,nz] = size(z);
37 minz = min(z);
38 maxz = max(z);
39 // Surround the matrix by a very low region to get closed contours, and
40 // replace any NaN with low numbers as well.
41 zz=[ %nan*ones(1,nz+2); %nan*ones(mz,1),z,%nan*ones(mz,1);%nan*ones(1,nz+2)];
42 kk=find(isnan(zz(:)));
43 zz(kk)=minz-1e4*(maxz-minz)+zeros(kk);
44
45 xx = [2*x(1)-x(2); x(:); 2*x(mz)-x(mz-1)];
46 yy = [2*y(1)-y(2); y(:); 2*y(nz)-y(nz-1)];
47
48 // Internal call to get the contours 
49 [x1,y1]=contour2di(xx,yy,zz,nv);
50 CS=[x1;y1];
51 // Find the indices of the curves in the c matrix, and get the
52 // area of closed curves in order to draw patches correctly. 
53 ii = 1;
54 ncurves = 0;
55 I = [];
56 Area=[];
57
58 while (ii < size(CS,2)),
59   nl=CS(2,ii);
60   ncurves = ncurves + 1;
61   I(ncurves) = ii;
62   xp=CS(1,ii+(1:nl));  // First patch
63   yp=CS(2,ii+(1:nl));
64   Area(ncurves)=sum( mtlb_diff(xp).*(yp(1:nl-1)+yp(2:nl))/2 );
65   ii = ii + nl + 1;
66 end
67
68 lp=xget('lastpattern');
69
70 if size(nv,'*') > 1 // case where nv is a vector defining the level curve values
71   if  size(nv,'*') > lp ; write(%io(2),'Colormap too small');return ;end 
72 else
73   if nv > lp ; write(%io(2),'Colormap too small');return ;end 
74 end
75
76 min_nv=mini(nv);
77 max_nv=maxi(nv);
78
79 plot2d([mini(xx);maxi(xx)],[mini(yy);maxi(yy)],0,strf,leg,rect,nax);
80
81 // Plot patches in order of decreasing size. This makes sure that
82 // all the lev1es get drawn, not matter if we are going up a hill or
83 // down into a hole. When going down we shift levels though, you can
84 // tell whether we are going up or down by checking the sign of the
85 // area (since curves are oriented so that the high side is always
86 // the same side). Lowest curve is largest and encloses higher data
87 // always.
88
89 draw_min=1;
90 H=[];
91 [FA,IA]=sort(abs(Area));
92
93   drawlater(); // postpon the drawing here
94   a=gca();
95   old_foreground = a.foreground;
96   pat=xget('pattern');
97   for jj=IA',
98     nl=CS(2,I(jj));
99     lev1=CS(1,I(jj));
100     if (lev1 ~= minz | draw_min),
101       xp=CS(1,I(jj)+(1:nl));  
102       yp=CS(2,I(jj)+(1:nl)); 
103       pat=size(find( nv <= lev1),'*');
104       xset("pattern",pat);
105       xfpoly(xp,yp)
106     end;
107   end
108   
109   if style(1)<>-1 then 
110     contour2d(xx,yy,zz,nv,style,"000",leg,rect,nax);
111   end
112   a.foreground = old_foreground;
113   drawnow(); // draw all now!
114
115 endfunction