9443d5cfc8370080af9163cc88937b5ebeeb0dbd
[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         newstyle = get('figure_style')=='new'
14         
15         if nin == 0 then   // demo
16                 t = -%pi:0.1:%pi;
17                 m = sin(t)'*cos(t);
18                 contourf(t,t,m);
19                 return;
20         end
21
22 if nin <= 8 then nax=[1,10,1,10];end 
23 if nin <= 7 then rect=[0,0,1,1];end
24 if nin <= 6 then leg=" ";end 
25 if nin <= 5 then strf="121";end 
26 if nin <= 3 then zmin=mini(z);zmax=maxi(z);nv = zmin + (1:10)*(zmax-zmin)/(11);end 
27 if nin <= 2 then z=rand(size(x,'*'),size(y,'*'));end 
28 if nin <= 1 then y=1:10;end 
29 if nin <= 0 then x=1:10;end 
30 if x==[] then x=1:size(z,'r');end 
31 if y==[] then y=1:size(z,'c');end 
32
33 nvs=size(nv,'*') ;
34 if nvs==1 then nvs=nv;zmin=mini(z);zmax=maxi(z);nv = zmin + (1:nvs)*(zmax-zmin)/(nvs+1);end;
35 if nin <= 4 then style = -1*ones(1,nvs);end
36 if nin <= 7 then rect=[mini(x),mini(y),maxi(x),maxi(y)]; end 
37 nv1=nv
38 [mz,nz] = size(z);
39 minz = min(z);
40 maxz = max(z);
41 // Surround the matrix by a very low region to get closed contours, and
42 // replace any NaN with low numbers as well.
43 zz=[ %nan*ones(1,nz+2); %nan*ones(mz,1),z,%nan*ones(mz,1);%nan*ones(1,nz+2)];
44 kk=find(isnan(zz(:)));
45 zz(kk)=minz-1e4*(maxz-minz)+zeros(kk);
46
47 xx = [2*x(1)-x(2); x(:); 2*x(mz)-x(mz-1)];
48 yy = [2*y(1)-y(2); y(:); 2*y(nz)-y(nz-1)];
49
50 // Internal call to get the contours 
51 [x1,y1]=contour2di(xx,yy,zz,nv);
52 CS=[x1;y1];
53 // Find the indices of the curves in the c matrix, and get the
54 // area of closed curves in order to draw patches correctly. 
55 ii = 1;
56 ncurves = 0;
57 I = [];
58 Area=[];
59
60 while (ii < size(CS,2)),
61   nl=CS(2,ii);
62   ncurves = ncurves + 1;
63   I(ncurves) = ii;
64   xp=CS(1,ii+(1:nl));  // First patch
65   yp=CS(2,ii+(1:nl));
66   Area(ncurves)=sum( mtlb_diff(xp).*(yp(1:nl-1)+yp(2:nl))/2 );
67   ii = ii + nl + 1;
68 end
69
70 lp=xget('lastpattern');
71
72 if size(nv,'*') > 1 // case where nv is a vector defining the level curve values
73   if  size(nv,'*') > lp ; write(%io(2),'Colormap too small');return ;end 
74 else
75   if nv > lp ; write(%io(2),'Colormap too small');return ;end 
76 end
77
78 min_nv=mini(nv);
79 max_nv=maxi(nv);
80
81 plot2d([mini(xx);maxi(xx)],[mini(yy);maxi(yy)],0,strf,leg,rect,nax);
82
83 // Plot patches in order of decreasing size. This makes sure that
84 // all the lev1es get drawn, not matter if we are going up a hill or
85 // down into a hole. When going down we shift levels though, you can
86 // tell whether we are going up or down by checking the sign of the
87 // area (since curves are oriented so that the high side is always
88 // the same side). Lowest curve is largest and encloses higher data
89 // always.
90
91 draw_min=1;
92 H=[];
93 [FA,IA]=sort(abs(Area));
94
95 if newstyle then
96   drawlater(); // postpon the drawing here
97   a=gca();
98   old_foreground = a.foreground;
99   pat=xget('pattern');
100   for jj=IA',
101     nl=CS(2,I(jj));
102     lev1=CS(1,I(jj));
103     if (lev1 ~= minz | draw_min),
104       xp=CS(1,I(jj)+(1:nl));  
105       yp=CS(2,I(jj)+(1:nl)); 
106       pat=size(find( nv <= lev1),'*');
107       xset("pattern",pat);
108       xfpoly(xp,yp)
109     end;
110   end
111   
112   if style(1)<>-1 then 
113     contour2d(xx,yy,zz,nv,style,"000",leg,rect,nax);
114   end
115   a.foreground = old_foreground;
116   drawnow(); // draw all now!
117 else
118   pat=xget('pattern');
119   for jj=IA',
120     nl=CS(2,I(jj));
121     lev1=CS(1,I(jj));
122     if (lev1 ~= minz | draw_min),
123       xp=CS(1,I(jj)+(1:nl));  
124       yp=CS(2,I(jj)+(1:nl)); 
125       pat=size(find( nv <= lev1),'*');
126       xset("pattern",pat);
127       xfpoly(xp,yp)
128     end;
129   end
130   
131   xset('pattern',pat);
132   if style(1)<>-1 then 
133     contour2d(xx,yy,zz,nv,style,"000",leg,rect,nax);
134   end
135 end
136
137 endfunction