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