quicwind/plot_wind_above_terrain.m

   1 function plot_wind_above_terrain(U,X,n_d)
   2 % plot_wind(U,X,d)
   3 %   U  cell array size 3 with the wind components
   4 %   X  cell array size 3 with the mesh coordinates
   5 %   n_d  row vector size 3, mesh display size
   6
   7 if ~exist('n_d','var')
   8     n_d = size(X{1});
   9 end
  10
  11 check_mesh(X);
  12
  13 figure,
  14 mesh(X{1}(:,:,1),X{2}(:,:,1),X{3}(:,:,1)); % terrain
  15
  16 [nx1,ny1,nz1]=size(X{1}); % mesh size
  17 nx=nx1-1; ny=ny1-1; nz=nz1-1;
  18 d = (size(X{1})-1)./n_d;
  19
  20 [ix,jx,kx]=meshgrid(1:nx1,1:ny1,1:nz1);           % grid indices of corners xy z
  21 [iu,ju,ku]=meshgrid(1:nx1,0.5+[1:ny],0.5+[1:nz]); % grid indices of u
  22 [iv,jv,kv]=meshgrid(0.5+[1:nx],1:ny1,0.5+[1:nz]); % grid indices of v
  23 [iw,jw,kw]=meshgrid(0.5+[1:nx],0.5+[1:ny],1:nz1);  % grid indices of w
  24
  25 px=1+d(1)/2:d(1):nx;
  26 py=1+d(2)/2:d(2):ny;
  27 pz=1+d(3)/2:d(3):nz;
  28
  29 [iq,jq,kq]=ndgrid(px,py,pz);  % query in index space
  30 nq=size(iq);
  31 iq=iq(:);
  32 jq=jq(:);
  33 kq=kq(:);
  34
  35 % interpolate coordinates and vectors to query
  36 p=[2,1,3];
  37 xd = interp3(ix,jx,kx,permute(X{1},p),iq,jq,kq);
  38 yd = interp3(ix,jx,kx,permute(X{2},p),iq,jq,kq);
  39 zd = interp3(ix,jx,kx,permute(X{3},p),iq,jq,kq);
  40 ud = interp3(iu,ju,ku,permute(U{1},p),iq,jq,kq);
  41 vd = interp3(iv,jv,kv,permute(U{2},p),iq,jq,kq);
  42 wd = interp3(iw,jw,kw,permute(U{3},p),iq,jq,kq);
  43
  44 hold on
  45 quiver3(xd,yd,zd,ud,vd,wd,0)
  46
  47 axis square
  48 xlabel('x'),ylabel('y'),zlabel('z')
  49 title('Wind field')
  50 a=gcf;fprintf('Figure %i: wind field 3D\n',a.Number)
  51 hold off
  52
  53 ud2 = reshape(ud,nq);
  54 wd2 = reshape(wd,nq);
  55 xd2 = reshape(xd,nq);
  56 zd2 = reshape(zd,nq);
  57
  58 mid = ceil(ny1/2);
  59 midq = ceil(nq(2)/2);
  60
  61 figure,
  62 plot(X{1}(:,mid,1),X{3}(:,mid,1))
  63 hold on
  64 quiver(xd2(:,midq,:),zd2(:,midq,:),ud2(:,midq,:),wd2(:,midq,:),0)
  65 xlabel('x'),ylabel('z')
  66 title('Wind field')
  67 a=gcf;fprintf('Figure %i: wind field 2D\n',a.Number)
  68 hold off
  69
  70 end
  71