scilab/modules/graphics/macros/polarplot.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) INRIA
   3 // Copyright (C) 2010 - DIGITEO - Manuel Juliachs
   4 // This file must be used under the terms of the CeCILL.
   5 // This source file is licensed as described in the file COPYING, which
   6 // you should have received as part of this distribution.  The terms
   7 // are also available at
   8 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
   9
  10 function polarplot(theta,rho,style,strf,leg,rect)
  11     [lhs,rhs]=argn(0)
  12     if rhs<=0 then
  13         theta=0:.01:2*%pi;
  14         rho=sin(2*theta).*cos(2*theta)
  15         scf();
  16         polarplot(theta,rho)
  17         return
  18     end
  19     if size(theta,1)==1 then
  20         theta=theta(:),
  21     end
  22     if size(rho,1)==1 then
  23         rho=rho(:),
  24     end
  25     rm=max(abs(rho))
  26     x=rho.*cos(theta)
  27     y=rho.*sin(theta)
  28
  29     opts=[]
  30     isstrf=%f;
  31     isframeflag=%f;
  32     isrect=%f;
  33     if exists("style","local")==1 then
  34         opts=[opts,"style=style"]
  35     end
  36     if exists("strf","local")==1 then
  37         opts=[opts,"strf=strf"]
  38         isstrf=%t
  39     end
  40     if exists("leg","local")==1 then
  41         opts=[opts,"leg=leg"]
  42     end
  43     if exists("rect","local")==1 then
  44         opts=[opts,"rect=rect"]
  45         isrect=%t
  46     end
  47     if exists("frameflag","local")==1 then
  48         opts=[opts,"frameflag=frameflag"]
  49         isframeflag=%t
  50     end
  51
  52     if size(opts,2)<rhs-2 then
  53         error(msprintf(gettext("%s: Wrong value for input argument: ''%s'', ''%s'', ''%s'', ''%s'' or ''%s'' expected.\n"),"polarplot","style","strf","leg","rect","frameflag"));
  54     end
  55
  56     // Some default values:
  57     Amin=0 // starting angle for the frame
  58     dA=360 // span of the angular frame
  59     nn=4    // number of quadrants to be drawn
  60
  61     xmin=min(x);
  62     xmax=max(x);
  63     L=(xmax-xmin)*1.07;
  64     ymin=min(y);
  65     ymax=max(y);
  66     H=(ymax-ymin)*1.07;
  67     // Angle at which Radial labels will be displayed
  68     A=round(atan((ymin+ymax)/2,(xmin+xmax)/2)/%pi*180/45)*45;
  69     dx=-0.5, dy=-0.5  // H & V shifts in string-width and string-height units
  70
  71     // Case without rect=
  72     if ~isrect then
  73         // Determines quadrant(s) to be drawn
  74         Q=[%T %T %T %T];
  75         e=rm/500;
  76
  77         if min(x)<-e then
  78             xmin=-rm;
  79         else
  80             xmin=0; Q([2 3])=%F;
  81         end
  82
  83         if max(x)>e then
  84             xmax= rm;
  85         else
  86             xmax=0; Q([1 4])=%F;
  87         end
  88
  89         if min(y)<-e then
  90             ymin=-rm;
  91         else
  92             ymin=0; Q([3 4])=%F;
  93         end
  94
  95         if max(y)>e then
  96             ymax= rm;
  97         else
  98             ymax=0; Q([1 2])=%F;
  99         end
 100
 101         L=(xmax-xmin)*1.1; if L==0, L=2*rm*1.1; end
 102         H=(ymax-ymin)*1.1; if H==0, H=2*rm*1.1; end
 103         x0=(xmin+xmax)/2; y0=(ymin+ymax)/2;
 104         rect=[x0-L/2 y0-H/2 x0+L/2 y0+H/2]
 105
 106         // Special case: data aligned on the x or y axis
 107         if Q==[%F %F %F %F],
 108             if (ymax-ymin)<2*e, // on x axis
 109                 if xmin<-e then
 110                     Q([2 3])=%T
 111                 end
 112                 if xmax> e  then
 113                     Q([1 4])=%T
 114                 end
 115             else  // on y axis
 116                 if ymin<-e  then
 117                     Q([3 4])=%T
 118                 end
 119                 if ymax> e then
 120                     Q([1 2])=%T
 121                 end
 122             end
 123         end
 124
 125         n=find(Q);   // id numbers of quadrants to be drawn
 126         nn=length(n) // number of quadrants to be drawn
 127         Amin=(n(1)-1)*90
 128
 129         select nn
 130         case 1,
 131             dA=90;
 132             if n==1, A=90, dx=-1.1, dy=-0.5
 133             elseif n==2, A=90, dx=0.2, dy=-0.5
 134             elseif n==3, A=270, dx=0.2, dy=-0.5
 135             else A=270, dx=-1.1, dy=-0.5
 136             end
 137         case 2
 138             dA=180;
 139             if n(1)==1
 140                 if n(2)==2, //A=90, dx=0.1, dy=-0.5
 141                 else Amin=-90, A=90, dx=-1.2, dy=-0.5, end
 142             elseif n(1)==2, A=90, dx=0.2, dy=-0.5
 143             else A=0, dx=-0.5, dy=0.2
 144             end
 145         else
 146             Amin=0, dA=360
 147         end
 148         opts=[opts,"rect=rect"]
 149     end // if ~isrect
 150
 151     if isstrf& isframeflag then
 152         error(msprintf(gettext("%s: ''%s'' cannot be used with ''%s''.\n"),"polarplot","frameflag","strf"));
 153     end
 154     if ~(isstrf) then
 155         axesflag=0
 156         opts=[opts,"axesflag=axesflag"],
 157     end
 158     if ~(isstrf|isframeflag) then
 159         frameflag=4
 160         opts=[opts,"frameflag=frameflag"],
 161     end
 162     drawlater()
 163     execstr("plot2d(x,y,"+strcat(opts,",")+")")
 164     a=gca();
 165     a.data_bounds=[rect(1:2);rect(3:4)]
 166     a.margins=[0.07 0.07 0.12 0.07]
 167
 168     fcolor=color("grey70");
 169     xset("dashes",1)
 170
 171     // CIRCULAR FRAME:
 172     // Radial values for the frame:
 173     fmt_in=format(), format("v",9)
 174     // Tunning for smart values:
 175     p=floor(log10(abs(rm)));
 176     m=rm/10^p;
 177     if m<1.3, dm=0.2
 178     elseif m<=2, dm=0.3
 179     elseif m<4, dm=0.5
 180     else dm=1,
 181     end
 182     k=fix(m/dm)
 183     if m-k*dm < dm/5, k=k-1, end
 184     R=[(1:k)*dm*10^p ]
 185     // Tuning for smart 10^ display using LaTeX instead of D+## exponential display
 186     if abs(p)<4,
 187         Rtxt=string(R)
 188         [v,k]=max(length(Rtxt)), tmp=xstringl(0,0,Rtxt(k))
 189     else
 190         if dm<1, dm=dm*10, p=p-1, end
 191         tmp=string(R/10^p)+"108"
 192         [v,k]=max(length(tmp)), tmp=xstringl(0,0,tmp(k))
 193         Rtxt="$\scriptstyle "+string(R/10^p)+"\:.10^{"+string(p)+"}$";
 194     end
 195     w=tmp(3); h=tmp(4);
 196     format(fmt_in(2),fmt_in(1))  // Restoring entrance format
 197     R = [ R  rm ]
 198
 199     // Drawing & labelling the radial frame
 200     kM=size(R,"*");
 201     for k=1:kM
 202         r=R(k)
 203         xarc(-r,r,2*r,2*r,Amin*64,dA*64)
 204         e = gce();,e.line_style=3
 205         e.foreground=fcolor;
 206         if k==kM
 207             e.line_style=1;  // solid outer arc
 208         else
 209             xstring(r*cosd(A)+w*dx, r*sind(A)+h*dy, Rtxt(k))
 210             e = gce();
 211             e.clip_state = "off";
 212         end
 213     end
 214
 215     // ANGULAR FRAME:
 216     if nn<3, eA=10, else eA=30; end // adaptative angular sampling
 217     an=linspace(Amin,Amin+dA,round(dA/eA)+1);
 218     // avoiding 360 == 0
 219     if nn>2, tmp=find(abs(an-360)<eA/10); an(tmp)=[]; end
 220     // Adjusting H-shifts of angular labels
 221     tmp=xstringl(0,0,"360");
 222     w=tmp(3); h=tmp(4);
 223     rL=rm*1.03;  // Radius of angular labels
 224     for k=an  // draws and labels angular rays
 225         xsegs([0;rm*cosd(k)],[0;rm*sind(k)])
 226         e = gce(); e.segs_color=fcolor; e.line_style=3;
 227         xstring((rL+w/2)*cosd(k)-w/2, (rL+h/2)*sind(k)-h/2, string(k))
 228         e = gce();
 229         e.clip_state = "off";
 230     end
 231     drawnow()
 232 endfunction