1## Copyright (C) 1998, 1999, 2000 Joao Cardoso.
2##
3## This program is free software; you can redistribute it and/or modify it
4## under the terms of the GNU General Public License as published by the
5## Free Software Foundation; either version 2 of the License, or (at your
6## option) any later version.
7##
8## This program is distributed in the hope that it will be useful, but
9## WITHOUT ANY WARRANTY; without even the implied warranty of
10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11## General Public License for more details.
12##
13## This file is part of plplot_octave.
14## It is based on the corresponding demo function of PLplot.
15
16## Does several contour plots using different coordinate mappings.
17
181;
19
20global	XPTS=35;		## Data points in x
21global	YPTS=46;		## Datat points in y
22
23global	tr = [2/(XPTS-1); 0.0; -1.0; 0.0; 2/(YPTS-1); -1.0];
24
25function [tx ty] = mypltr(x, y)
26
27  global	XPTS
28  global	YPTS
29  global	tr
30
31  tx = tr(1) * x + tr(2) * y + tr(3);
32  ty = tr(4) * x + tr(5) * y + tr(6);
33
34endfunction
35
36
37function _polar
38
39  plenv(-1., 1., -1., 1., 0, -2);
40  plcol0(1);
41
42  ## Perimeter
43  PERIMETERPTS = 100;
44  RPTS = 40;
45  THETAPTS = 40;
46
47  i = 0:PERIMETERPTS-1;
48  t = (2.*pi/(PERIMETERPTS-1))*i;
49  px = cos(t);
50  py = sin(t);
51  plline(px', py');
52
53  ## create data to be contoured
54  i = 0:RPTS-1;
55  r = i/(RPTS-1);
56  z = r'*ones(1,RPTS);
57
58  j = (0:THETAPTS-1)';
59  theta = (2.*pi/(THETAPTS-1))*j;
60  xg = cos(theta)*r;
61  yg = sin(theta)*r;
62
63  i = 0:9;
64  lev = 0.05 + 0.10* i';
65
66  plcol0(2);
67  plcont2(z, 1, RPTS, 1, THETAPTS, lev, xg', yg');
68  plcol0(1);
69  pllab("", "", "Polar Contour Plot");
70endfunction
71
72## shielded potential contour plot example.
73
74function potential
75
76  PPERIMETERPTS = 100;
77  PRPTS = 40;
78  PTHETAPTS = 64;
79  PNLEVEL = 20;
80
81  ## create data to be contoured
82
83  i = 0:PRPTS-1;
84  r = 0.5 + i;
85  j = (0:PTHETAPTS-1)';
86  theta = (2.*pi/(PTHETAPTS-1))*(0.5 + j);
87  xg = (cos(theta)*r)';
88  yg = (sin(theta)*r)';
89
90  rmax = max(r);
91  xmin = min(min(xg));
92  xmax = max(max(xg));
93  ymin = min(min(yg));
94  ymax = max(max(yg));
95
96  x0 = (xmin + xmax)/2;
97  y0 = (ymin + ymax)/2;
98
99  ## Expanded limits
100  peps = 0.05;
101  xpmin = xmin - abs(xmin)*peps;
102  xpmax = xmax + abs(xmax)*peps;
103  ypmin = ymin - abs(ymin)*peps;
104  ypmax = ymax + abs(ymax)*peps;
105
106  ## Potential inside a conducting cylinder (or sphere) by method of images.
107  ## Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
108  ## Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
109  ## Also put in smoothing term at small distances.
110
111  eeps = 2.;
112
113  q1 = 1.;
114  d1 = rmax/4.;
115
116  q1i = - q1*rmax/d1;
117  d1i = rmax^2/d1;
118
119  q2 = -1.;
120  d2 = rmax/4.;
121
122  q2i = - q2*rmax/d2;
123  d2i = rmax^2/d2;
124
125  div1 = sqrt((xg-d1).^2 + (yg-d1).^2 + eeps^2);
126  div1i = sqrt((xg-d1i).^2 + (yg-d1i).^2 + eeps^2);
127  div2 = sqrt((xg-d2).^2 + (yg+d2).^2 + eeps^2);
128  div2i = sqrt((xg-d2i).^2 + (yg+d2i).^2 + eeps^2);
129  z = q1./div1 + q1i./div1i + q2./div2 + q2i./div2i;
130
131  zmin = min(min(z));
132  zmax = max(max(z));
133
134  ## Positive and negative contour levels.
135  dz = (zmax-zmin)/ PNLEVEL;
136
137  i = 0:PNLEVEL-1;
138  clevel = zmin + ( i + 0.5)*dz;
139  clevelneg = clevel(clevel <= 0);
140  clevelpos = clevel(clevel > 0);
141  nlevelneg = columns( clevelneg);
142  nlevelpos = columns( clevelpos);
143
144  ## Colours!
145  ncollin = 11;
146  ncolbox = 1;
147  ncollab = 2;
148
149  ## Finally start plotting this page!
150  pladv(0);
151  plcol0(ncolbox);
152
153  plvpas(0.1, 0.9, 0.1, 0.9, 1.0);
154  plwind(xpmin, xpmax, ypmin, ypmax);
155  plbox("", 0., 0, "", 0., 0);
156
157  plcol0(ncollin);
158  if(nlevelneg >0)
159    ## Negative contours
160    pllsty(2);
161    plcont2(z, 1, PRPTS, 1, PTHETAPTS, clevelneg', xg, yg);
162  endif
163
164  if(nlevelpos >0)
165    ## Positive contours
166    pllsty(1);
167    plcont2(z, 1, PRPTS, 1, PTHETAPTS, clevelpos', xg, yg);
168  endif
169
170  ## Draw outer boundary
171  i = 0:PPERIMETERPTS-1;
172  t = (2.*pi/(PPERIMETERPTS-1))*i;
173  px = x0 + rmax*cos(t);
174  py = y0 + rmax*sin(t);
175
176  plcol0(ncolbox);
177  plline(px', py');
178
179  plcol0(ncollab);
180  pllab("", "", "Shielded potential of charges in a conducting sphere");
181
182endfunction
183
184function ix09c
185
186	  global	XPTS
187	  global	YPTS
188	  global	tr
189
190	  clevel = linspace(-1,1,11)';
191	  ## Hack to ensure that zero contour really is zero
192	  ## and not a very small number.
193	  ## For me problem only occurs on i386 32 bit Debian 3.1
194	  ## with octave 2.1.69.
195	  ## 64-bit Ubuntu Dapper with octave 2.1.72 seems ok.
196	  clevel(6) = 0.0;
197
198          mark0 = []; space0 = [];
199	  mark1 = [1500]; space1 = [1500];
200
201	  ## Parse and process command line arguments
202
203	  ##    (void) plparseopts(&argc, argv, PL_PARSE_FULL);
204
205	  ## Initialize plplot
206
207	  plinit();
208
209	  ## Set up function arrays
210
211	  for i=0:XPTS-1
212	    xx = (i - fix(XPTS / 2)) /  fix(XPTS / 2);
213	    yy = ((0:YPTS-1) - fix(YPTS / 2)) / fix(YPTS / 2) - 1.0;
214	    z(i+1,:) = xx * xx - yy .* yy;
215	    w(i+1,:) = 2 * xx .* yy;
216	  endfor
217
218	  ## Plot using identity transform
219if (0)
220	  pl_setcontlabelparam(0.006, 0.3, 0.1, 0);
221	  plenv(-1.0, 1.0, -1.0, 1.0, 0, 0);
222	  plcol0(2);
223	  plcont(z, 1, XPTS, 1, YPTS, clevel, tr);
224	  plstyl(mark1, space1);
225	  plcol0(3);
226	  plcont(w, 1, XPTS, 1, YPTS, clevel, tr);
227	  plstyl(mark0, space0);
228	  plcol0(1);
229	  pllab("X Coordinate", "Y Coordinate", "Streamlines of flow");
230endif
231	  pl_setcontlabelformat(4,3);
232	  pl_setcontlabelparam(0.006, 0.3, 0.1, 1);
233	  plenv(-1.0, 1.0, -1.0, 1.0, 0, 0);
234	  plcol0(2);
235	  plcont(z, 1, XPTS, 1, YPTS, clevel, tr);
236	  plstyl(mark1, space1);
237	  plcol0(3);
238	  plcont(w, 1, XPTS, 1, YPTS, clevel, tr);
239	  plstyl(mark0, space0);
240	  plcol0(1);
241	  pllab("X Coordinate", "Y Coordinate", "Streamlines of flow");
242
243	  ## Set up grids
244
245	  for i=0:XPTS-1
246	    [xx yy] = mypltr(i, (0:YPTS-1));
247
248	    argx = xx * pi/2;
249	    argy = yy * pi/2;
250	    distort = 0.4;
251
252	    xg3(i+1,:) = xx .+ distort .* cos(argx);
253	    yg3(i+1,:) = yy .- distort .* cos(argy);
254
255	    xg2(i+1,:) = xx .+ distort .* cos(argx) .* cos(argy);
256	    yg2(i+1,:) = yy .- distort .* cos(argx) .* cos(argy);
257	  endfor
258
259	  xg1 = xg3(:,1);
260	  yg1 = yg3(XPTS,:)';
261
262
263	  ## Plot using 1d coordinate transform
264
265	  pl_setcontlabelparam(0.006, 0.3, 0.1, 0);
266	  plenv(-1.0, 1.0, -1.0, 1.0, 0, 0);
267	  plcol0(2);
268	  plcont1(z, 1, XPTS, 1, YPTS, clevel, xg1, yg1);
269	  plstyl(mark1, space1);
270	  plcol0(3);
271	  plcont1(w, 1, XPTS, 1, YPTS, clevel, xg1, yg1);
272	  plstyl(mark0, space0);
273	  plcol0(1);
274	  pllab("X Coordinate", "Y Coordinate", "Streamlines of flow");
275if(0)
276	  pl_setcontlabelparam(0.006, 0.3, 0.1, 1);
277	  plenv(-1.0, 1.0, -1.0, 1.0, 0, 0);
278	  plcol0(2);
279	  plcont1(z, 1, XPTS, 1, YPTS, clevel, xg1, yg1);
280	  plstyl(mark1, space1);
281	  plcol0(3);
282	  plcont1(w, 1, XPTS, 1, YPTS, clevel, xg1, yg1);
283	  plstyl(mark0, space0);
284	  plcol0(1);
285	  pllab("X Coordinate", "Y Coordinate", "Streamlines of flow");
286endif
287	  ## Plot using 2d coordinate transform
288
289	  pl_setcontlabelparam(0.006, 0.3, 0.1, 0);
290	  plenv(-1.0, 1.0, -1.0, 1.0, 0, 0);
291	  plcol0(2);
292	  plcont2(z, 1, XPTS, 1, YPTS, clevel, xg2, yg2);
293	  plstyl(mark1, space1);
294	  plcol0(3);
295	  plcont2(w, 1, XPTS, 1, YPTS, clevel, xg2, yg2);
296	  plstyl(mark0, space0);
297	  plcol0(1);
298	  pllab("X Coordinate", "Y Coordinate", "Streamlines of flow");
299if(0)
300	  pl_setcontlabelparam(0.006, 0.3, 0.1, 1);
301	  plenv(-1.0, 1.0, -1.0, 1.0, 0, 0);
302	  plcol0(2);
303	  plcont1(z, 1, XPTS, 1, YPTS, clevel, xg1, yg1);
304	  plstyl(mark1, space1);
305	  plcol0(3);
306	  plcont1(w, 1, XPTS, 1, YPTS, clevel, xg1, yg1);
307	  plstyl(mark0, space0);
308	  plcol0(1);
309	  pllab("X Coordinate", "Y Coordinate", "Streamlines of flow");
310endif
311	  pl_setcontlabelparam(0.006, 0.3, 0.1, 0);
312	  _polar();
313	  ## pl_setcontlabelparam(0.006, 0.3, 0.1, 1);
314	  ## _polar();
315
316
317    pl_setcontlabelparam(0.006, 0.3, 0.1, 0);
318    potential();
319    ## pl_setcontlabelparam(0.006, 0.3, 0.1, 1);
320    ## potential();
321
322
323	  plend1();
324
325	endfunction
326
327ix09c
328