vis3d/plot_wind_at_h.m

   1 function plot_wind_at_h(p,iheights,levels,alpha,timestep,windscale,heightscale)
   2 x=p.xlat;
   3 y=p.xlong;
   4 clf
   5 if ~exist('timestep','var'),
   6     timestep=1;
   7 end
   8 if ~isscalar(timestep),
   9     error('timestep must be scalar')
  10 end
  11 hgt=p.hgt(:,:,timestep); % terrain height
  12 % wind at given height
  13 for i=iheights(:)'
  14     z=p.heights(i)*ones(size(x));
  15     ws=sqrt(p.uch(:,:,i,timestep).^2 + p.vch(:,:,i,timestep).^2);
  16     print_layer_msg(i,z,ws)
  17     surf(x,y,z+hgt,ws,'EdgeAlpha',alpha*0,'FaceAlpha',alpha)
  18     hold on
  19 end
  20 for i=levels(:)'
  21     z=p.altitude(:,:,i,timestep)-hgt;
  22     ws=sqrt(p.uc(:,:,i,timestep).^2 + p.vc(:,:,i,timestep).^2);
  23     print_layer_msg(i,z,ws)
  24     surf(x,y,z+hgt,ws,'EdgeAlpha',alpha*0,'FaceAlpha',alpha)
  25     hold on
  26 end
  27 caxis([0,windscale])
  28 title({'Horizontal wind speed',p.times{1}},'Interpreter','none')
  29 xlabel('m'), ylabel('m')  % for now, ideal only
  30 h=colorbar
  31 xlabel(h,'m/s')
  32 axis('image')
  33 a=[min(x(:)),max(x(:)),min(y(:)),max(y(:)),0,heightscale]; % var limits
  34 axis(a);     % set var limits
  35 xysize=max([a(4)-a(3),a(2)-a(1)]);
  36 zsize=a(6)-a(5);
  37 daspect([1,1,3*zsize/xysize]);
  38 hold off, drawnow
  39 end
  40
  41 function print_layer_msg(i,z,ws)
  42     fprintf('level %i height %g to %g median %g horizontal wind speed %g to %g median %g\n',...
  43        i,min(z(:)),max(z(:)),median(z(:)),min(ws(:)),max(ws(:)),median(ws(:)) )
  44 end