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