quicwind/plot_wind.m

   1 function plot_wind(u3,h3,d)
   2 % plot_wind(u3,h3,s)
   3 %   u3  cell array size 3 with the wind components
   4 %   h3  vector of mesh spacing (m)
   5 %   d   (optional) distance between points to display, in all 3 directions
   6
   7 if ~exist('d','var')
   8     d = h3(3)*(size(u3{1},3)-1)/3;
   9 end
  10
  11 % interpolate to center points
  12 u = 0.5*(u3{1}(1:end-1,:,:)+u3{1}(2:end,:,:));
  13 v = 0.5*(u3{2}(:,1:end-1,:)+u3{2}(:,2:end,:));
  14 w = 0.5*(u3{3}(:,:,1:end-1)+u3{3}(:,:,2:end));
  15
  16 % dimension of the center point array
  17 n=size(u);
  18
  19 for i=1:3
  20     c{i}=h3(i)*(1/2 + [1:n(i)]);     % coordinates of s
  21     l(i) = h3(i)*(n(i)-1);           % length of the dimension
  22     p{i} = [d/2 : d : l(i)];         % display points in physical space [0, h3(i)*(n(i)-1)]
  23     q{i} = 1+p{i}/h3(i);             % display points in index space [1,n(i)]
  24 end
  25
  26
  27 [x, y ,z ] = ndgrid(c{1},c{2},c{3});
  28 [xd,yd,zd] = ndgrid(p{1},p{2},p{3});
  29 [iq,jq,kq] = ndgrid(q{1},q{2},q{3});
  30
  31 ud = interp3(u,iq,jq,kq);
  32 vd = interp3(v,iq,jq,kq);
  33 wd = interp3(w,iq,jq,kq);
  34
  35 quiver3(xd,yd,zd,ud,vd,wd,0.5)
  36
  37 disp('for some reason, streamlines do not work')
  38 %hold on
  39 %[startx,starty,startz]=ndgrid(p{1}(1),p{2}(1),p{3}(1));
  40 %streamline(x,y,z,u,v,w,startx,starty,startz)
  41 %hold off
  42
  43 axis([0,l(1),0,l(2),0,l(3)])
  44 axis equal
  45 xlabel('m'),ylabel('m'),zlabel('m')
  46 title('Wind field')
  47 a=gcf;fprintf('Figure %i: wind field\n',a.Number)
  48 mid=round(size(xd,2)/2);
  49 figure
  50 quiver(xd(:,mid,:),zd(:,mid,:),ud(:,mid,:),vd(:,mid,:),0.5)
  51 axis([0,l(1),0,l(3)])
  52 axis equal
  53 xlabel('m'),ylabel('m'),zlabel('m')
  54 title('Wind field at middle crossection')
  55 a=gcf;fprintf('Figure %i: wind field middle crossection\n',a.Number)
  56
  57
  58 end