1## Copyright (C) 2004-2006 Andrew Ross
2## Copyright (C) 2004 Rafael Laboissiere
3##
4## This program is free software; you can redistribute it and/or modify it
5## under the terms of the GNU General Public License as published by the
6## Free Software Foundation; either version 2 of the License, or (at your
7## option) any later version.
8##
9## This program is distributed in the hope that it will be useful, but
10## WITHOUT ANY WARRANTY; without even the implied warranty of
11## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12## General Public License for more details.
13##
14## This file is part of plplot_octave.
15## It is based on the corresponding demo function of PLplot.
16
17# Simple line plot and multiple windows demo.
18
191;
20
21global xmax;
22
23function ix22c
24
25  ## Parse and process command line arguments */
26
27  ##    plparseopts(&argc, argv, PL_PARSE_FULL);
28
29  ## Initialize plplot */
30  plinit;
31
32  ## Set up the data */
33  ## Original case */
34
35  circulation;
36
37  narr = 6;
38  fill = 0;
39
40  arrow_x = [-0.5 0.5 0.3 0.5 0.3 0.5];
41  arrow_y = [0.0 0.0 0.2 0.0 -0.2 0.0];
42  arrow2_x = [-0.5 0.3 0.3 0.5 0.3 0.3];
43  arrow2_y = [0.0 0.0   0.2 0.0 -0.2 0.0];
44
45  ## Set arrow style using arrow_x and arrow_y then
46  ## plot using these arrows.
47  plsvect(arrow_x', arrow_y', fill);
48  constriction(1);
49
50  ## Set arrow style using arrow2_x and arrow2_y then
51  ## plot using these filled arrows. */
52  fill = 1;
53  plsvect(arrow2_x', arrow2_y', fill);
54  constriction(2);
55
56  constriction2;
57
58  plsvect([],[],0);
59
60  potential;
61
62  ## Don't forget to call plend1  to finish off! */
63
64  plend1();
65
66endfunction
67
68function circulation
69  nx = 20;
70  ny = 20;
71
72  dx = 1.0;
73  dy = 1.0;
74
75  xmin = -nx/2*dx;
76  xmax = nx/2*dx;
77  ymin = -ny/2*dy;
78  ymax = ny/2*dy;
79
80  xg = [xmin+dx/2:dx:xmax-dx/2]'*ones(1,ny);
81  yg = ones(nx,1)*[ymin+dy/2:dy:ymax-dy/2];
82  u = yg;
83  v = -xg;
84
85  ## Plot vectors with default arrows
86  plenv(xmin, xmax, ymin, ymax, 0, 0);
87  pllab("(x)", "(y)", "#frPLplot Example 22 - circulation");
88  plcol0(2);
89  plvect2(u,v,0.0,xg,yg);
90  plcol0(1);
91
92end
93
94## Vector plot of flow through a constricted pipe
95function constriction( astyle )
96  nx = 20;
97  ny = 20;
98
99  dx = 1.0;
100  dy = 1.0;
101
102  xmin = -nx/2*dx;
103  xmax = nx/2*dx;
104  ymin = -ny/2*dy;
105  ymax = ny/2*dy;
106
107  Q = 2.0;
108  xg = [xmin+dx/2:dx:xmax-dx/2]'*ones(1,ny);
109  yg = ones(nx,1)*[ymin+dy/2:dy:ymax-dy/2];
110
111  b = ymax/4.0.*(3-cos(pi*xg/xmax));
112  dbdx = ymax/4.0.*sin(pi*xg/xmax)*pi/xmax.*yg./b;
113  u = Q*ymax./b.*(abs(yg)<b);
114  v = dbdx.*u.*(abs(yg)<b);
115
116  plenv(xmin, xmax, ymin, ymax, 0, 0);
117  title = sprintf( "#frPLplot Example 22 - constriction (arrow style %d)", astyle );
118
119  pllab("(x)", "(y)", title );
120  plcol0(2);
121  plvect2(u,v,-1.0,xg,yg);
122  plcol0(1);
123
124end
125
126##
127## Global transform function for a constriction using data passed in
128## This is the same transformation used in constriction.
129##
130function [xt, yt] = transform( x, y, data )
131    global xmax;
132    xt = x;
133    yt = y / 4.0 * ( 3 - cos( pi * x / xmax ) );
134end
135
136## Vector plot of flow through a constricted pipe with
137## a coordinate transformation
138function constriction2()
139
140  global xmax;
141
142  nx = 20;
143  ny = 20;
144  nc = 11;
145  nseg = 20;
146
147  dx = 1.0;
148  dy = 1.0;
149
150  xmin = -nx/2*dx;
151  xmax = nx/2*dx;
152  ymin = -ny/2*dy;
153  ymax = ny/2*dy;
154
155  plstransform( @transform, [] );
156
157  Q = 2.0;
158  xg = [xmin+dx/2:dx:xmax-dx/2]'*ones(1,ny);
159  yg = ones(nx,1)*[ymin+dy/2:dy:ymax-dy/2];
160
161  b = ymax/4.0.*(3.0-cos(pi*xg/xmax));
162  u = Q*ymax./b;
163  v = zeros(nx,ny);
164
165  clev = Q + (0:(nc-1))*Q/(nc-1);
166
167  plenv(xmin, xmax, ymin, ymax, 0, 0);
168  pllab("(x)", "(y)", "#frPLplot Example 22 - constriction with plstransform" );
169  plcol0(2);
170  plshades(u, xmin+dx/2,xmax-dx/2,ymin+dy/2,ymax-dy/2,clev',0.0,1,1.0,0);
171  plvect2(u,v,-1.0,xg,yg);
172  ## Plot edges using plpath (which accounts for coordinate transformation) rather than plline
173  plpath( nseg, xmin, ymax, xmax, ymax );
174  plpath( nseg, xmin, ymin, xmax, ymin );
175  plcol0(1);
176
177  plstransform( [], [] );
178
179end
180
181## Vector plot of the gradient of a shielded potential (see example 9)
182
183function potential
184  nper = 100;
185  nlevel = 10;
186  nr = 20;
187  ntheta = 20;
188
189  ## Potential inside a conducting cylinder (or sphere) by method of images.
190  ## Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
191  ## Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
192  ## Also put in smoothing term at small distances.
193
194  rmax = nr;
195
196  eps = 2.;
197
198  q1 = 1.;
199  d1 = rmax/4.;
200
201  q1i = - q1*rmax/d1;
202  d1i = rmax^2/d1;
203
204  q2 = -1.;
205  d2 = rmax/4.;
206
207  q2i = - q2*rmax/d2;
208  d2i = rmax^2/d2;
209
210  r = [0.5:1:nr-0.5];
211  dtheta = 2*pi/(ntheta-1);
212  theta = [dtheta/2:dtheta:2*pi+dtheta/2];
213  x = r'*cos(theta);
214  y = r'*sin(theta);
215  div1 = sqrt((x-d1).^2 + (y-d1).^2 + eps^2);
216  div1i = sqrt((x-d1i).^2 + (y-d1i).^2 + eps^2);
217  div2 = sqrt((x-d2).^2 + (y+d2).^2 + eps^2);
218  div2i = sqrt((x-d2i).^2 + (y+d2i).^2 + eps^2);
219  z = q1*div1.^-1 + q1i*div1i.^-1 + q2*div2.^-1 + q2i*div2i.^-1;
220  u = -q1*(x-d1)./div1.^3 - q1i*(x-d1i)./div1i.^3 - q2*(x-d2)./div2.^3 - q2i*(x-d2i)./div2i.^3;
221  v = -q1*(y-d1)./div1.^3 - q1i*(y-d1i)./div1i.^3 - q2*(y+d2)./div2.^3 - q2i*(y+d2i)./div2i.^3;
222
223  xmin = min(min(x));
224  xmax = max(max(x));
225  ymin = min(min(y));
226  ymax = max(max(y));
227  zmin = min(min(z));
228  zmax = max(max(z));
229
230  plenv(xmin, xmax, ymin, ymax, 0, 0);
231  pllab("(x)", "(y)", "#frPLplot Example 22 - potential gradient vector plot");
232  ## Plot contours of the potential
233  dz = (zmax-zmin)/nlevel;
234  clevel = [zmin+dz/2:dz:zmax-dz/2];
235  plcol0(3);
236  pllsty(2);
237  plcont2(z,1,nr,1,ntheta,clevel',x,y);
238  pllsty(1);
239  plcol0(1);
240
241  ## Plot the vectors of the gradient of the potential
242  plcol0(2);
243  plvect2(u,v,25.0,x,y);
244  plcol0(1);
245
246  ## Plot the perimeter of the cylinder
247  theta = [0:2*pi/(nper-1):2*pi];
248  px = rmax*cos(theta);
249  py = rmax*sin(theta);
250  plline(px',py');
251end
252
253ix22c
254
255
256