File: polarplot.sci

package info (click to toggle)
scilab 5.3.3-10
  • links: PTS, VCS
  • area: main
  • in suites: wheezy
  • size: 330,656 kB
file content (226 lines) | stat: -rw-r--r-- 5,868 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA
// Copyright (C) 2010 - DIGITEO - Manuel Juliachs
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at    
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

function polarplot(theta,rho,style,strf,leg,rect)
  [lhs,rhs]=argn(0)
  if rhs<=0 then
    theta=0:.01:2*%pi;
    rho=sin(2*theta).*cos(2*theta)
    clf();
    polarplot(theta,rho)
    return
  end
  if size(theta,1)==1 then
    theta=theta(:),
  end
  if size(rho,1)==1 then
    rho=rho(:),
  end
  rm=max(abs(rho))
  x=rho.*cos(theta)
  y=rho.*sin(theta)

  opts=[]
  isstrf=%f;
  isframeflag=%f;
  isrect=%f;
  if exists('style','local')==1 then
    opts=[opts,'style=style']
  end
  if exists('strf','local')==1 then
    opts=[opts,'strf=strf']
    isstrf=%t
  end
  if exists('leg','local')==1 then 
    opts=[opts,'leg=leg']
  end
  if exists('rect','local')==1 then
    opts=[opts,'rect=rect']
    isrect=%t
  end
  if exists('frameflag','local')==1 then
    opts=[opts,'frameflag=frameflag']
    isframeflag=%t
  end

  if size(opts,2)<rhs-2 then
    error(msprintf(gettext("%s: Wrong value for input argument: ''%s'', ''%s'', ''%s'', ''%s'' or ''%s'' expected.\n"),"polarplot","style","strf","leg","rect","frameflag"));
  end
  
  // Some default values:
  Amin=0 // starting angle for the frame
  dA=360 // span of the angular frame
  nn=4    // number of quadrants to be drawn
  
  xmin=min(x);
  xmax=max(x);
  L=(xmax-xmin)*1.07;
  ymin=min(y);
  ymax=max(y);
  H=(ymax-ymin)*1.07;
  // Angle at which Radial labels will be displayed
  A=round(atan((ymin+ymax)/2,(xmin+xmax)/2)/%pi*180/45)*45;
  dx=-0.5, dy=-0.5  // H & V shifts in string-width and string-height units

  // Case without rect=
  if ~isrect then
    // Determines quadrant(s) to be drawn
    Q=[%T %T %T %T];
    e=rm/500;
    
    if min(x)<-e then
      xmin=-rm; 
    else
      xmin=0; Q([2 3])=%F; 
    end
    
    if max(x)>e then
      xmax= rm; 
    else 
      xmax=0; Q([1 4])=%F; 
    end
    
    if min(y)<-e then
      ymin=-rm; 
    else 
      ymin=0; Q([3 4])=%F; 
    end
    
    if max(y)>e then
      ymax= rm; 
    else 
      ymax=0; Q([1 2])=%F; 
    end 
    
    L=(xmax-xmin)*1.1; if L==0, L=2*rm*1.1; end
    H=(ymax-ymin)*1.1; if H==0, H=2*rm*1.1; end
    x0=(xmin+xmax)/2; y0=(ymin+ymax)/2;
    rect=[x0-L/2 y0-H/2 x0+L/2 y0+H/2]
    
    // Special case: data aligned on the x or y axis
    if Q==[%F %F %F %F],
        if (ymax-ymin)<2*e, // on x axis
          if xmin<-e then
            Q([2 3])=%T
          end
          if xmax> e  then
            Q([1 4])=%T
          end
        else  // on y axis
          if ymin<-e  then
            Q([3 4])=%T
          end
          if ymax> e then 
            Q([1 2])=%T 
          end
        end
    end
      
    n=find(Q);   // id numbers of quadrants to be drawn
    nn=length(n) // number of quadrants to be drawn
    Amin=(n(1)-1)*90

    select nn
      case 1, 
        dA=90;
        if n==1, A=90, dx=-1.1, dy=-0.5
        elseif n==2, A=90, dx=0.2, dy=-0.5
        elseif n==3, A=270, dx=0.2, dy=-0.5
        else A=270, dx=-1.1, dy=-0.5 
        end
      case 2
        dA=180;
        if n(1)==1
          if n(2)==2, //A=90, dx=0.1, dy=-0.5 
          else Amin=-90, A=90, dx=-1.2, dy=-0.5, end
        elseif n(1)==2, A=90, dx=0.2, dy=-0.5 
        else A=0, dx=-0.5, dy=0.2
        end
      else 
         Amin=0, dA=360
      end
    opts=[opts,'rect=rect']
  end // if ~isrect

  if isstrf& isframeflag then
    error(msprintf(gettext("%s: ''%s'' cannot be used with ''%s''.\n"),"polarplot","frameflag","strf"));
  end
  if ~(isstrf) then
    axesflag=0
    opts=[opts,'axesflag=axesflag'],
  end
  if ~(isstrf|isframeflag) then
    frameflag=4
    opts=[opts,'frameflag=frameflag'],
  end
  drawlater()
  execstr('plot2d(x,y,'+strcat(opts,',')+')')
  a=gca(); 
  a.data_bounds=[rect(1:2);rect(3:4)]
  a.margins=[0.07 0.07 0.12 0.07]
 
  fcolor=color("grey70");
  xset("dashes",1)
  
  // CIRCULAR FRAME:
  // Radial values for the frame:
  fmt_in=format(), format("v",9)
  // Tunning for smart values:
  p=floor(log10(abs(rm)));
  m=rm/10^p;
  if m<1.3, dm=0.2
    elseif m<=2, dm=0.3
    elseif m<4, dm=0.5
    else dm=1, 
  end
  k=fix(m/dm)
  if m-k*dm < dm/5, k=k-1, end
  R=[(1:k)*dm*10^p ]
  // Tuning for smart 10^ display using LaTeX instead of D+## exponential display
  if abs(p)<4, 
    Rtxt=string(R)
    [v,k]=max(length(Rtxt)), tmp=xstringl(0,0,Rtxt(k))
  else 
    if dm<1, dm=dm*10, p=p-1, end
    tmp=string(R/10^p)+"108"
    [v,k]=max(length(tmp)), tmp=xstringl(0,0,tmp(k))
    Rtxt="$\scriptstyle "+string(R/10^p)+"\:.10^{"+string(p)+"}$";
  end
  w=tmp(3); h=tmp(4);
  format(fmt_in(2),fmt_in(1))  // Restoring entrance format
  R = [ R  rm ]
  
  // Drawing & labelling the radial frame
  kM=size(R,"*");
  for k=1:kM
    r=R(k)      
    xarc(-r,r,2*r,2*r,Amin*64,dA*64)
    e = gce();,e.line_style=3
    e.foreground=fcolor;
    if k==kM, e.line_style=1;  // solid outer arc
    else xstring(r*cosd(A)+w*dx, r*sind(A)+h*dy, Rtxt(k))
    end
  end
  
  // ANGULAR FRAME:
  if nn<3, eA=10, else eA=30; end // adaptative angular sampling
  an=linspace(Amin,Amin+dA,round(dA/eA)+1);
  // avoiding 360 == 0
  if nn>2, tmp=find(abs(an-360)<eA/10); an(tmp)=[]; end
  // Adjusting H-shifts of angular labels
  tmp=xstringl(0,0,'360');
  w=tmp(3); h=tmp(4); 
  rL=rm*1.03;  // Radius of angular labels
  for k=an  // draws and labels angular rays
    xsegs([0;rm*cosd(k)],[0;rm*sind(k)])
    e = gce(); e.segs_color=fcolor; e.line_style=3;
    xstring((rL+w/2)*cosd(k)-w/2, (rL+h/2)*sind(k)-h/2, string(k))
  end
  drawnow()
endfunction