vis3d/frame3d.m

   1 function frame3d(swind,amin,amax,astep,qstep,qs,...
   2     fxlong,fxlat,xlong,xlat,zsf,fgrnhfx,fuel_frac,uf,vf,u,v,w,ph,phb,hgt)
   3
   4 clf,hold off
   5
   6 r = size(fxlong)./(size(xlong));  % refinement ratio
   7 % see if xlong and xlat are bogus
   8 ideal = all(xlong(:) ==0)|any(abs(xlong(:))>180)|any(abs(xlat(:))>180);
   9
  10 if ideal,
  11     % do not have xlong and xlat in ideal case, but we made the coordinates up
  12     % for fxlong and fxlat, so interpolate
  13     n=size(xlong);
  14     [xa,ya]=ndgrid(r(1)*[1/2:1:n(1)],r(2)*[1/2:1:n(2)]); % centers of atm cells in fire grid
  15     axlong = interp2(fxlong,xa,ya); % fake coordinates for display
  16     axlat = interp2(fxlat,xa,ya); % fake coordinates for display
  17     xscale=1;
  18     yscale=1;
  19 else
  20     axlong=xlong; axlat=xlat; % in real case, xlong xlat populated
  21         er=6378e3; % earth radius in m
  22     deg2rad=2*pi/360; % 1 degree in radians
  23     % degrees to meters, long and lat
  24     lat_deg2m=er*deg2rad;
  25     long_deg2m=er*deg2rad*cos(deg2rad*mean(xlat(find(xlat))));
  26     xscale=1/long_deg2m;
  27     yscale=1/lat_deg2m;
  28 end
  29 aspect_ratio = [xscale yscale 1];
  30
  31 % compute the corresponding part of the fire grid
  32 dmin=1+r.*(amin(1:2)-1);
  33 dmax=r.*amax(1:2);
  34
  35 % surface colored by heat flux and fuel fraction
  36 i=dmin(1):dmax(1);j=dmin(2):dmax(2);
  37 hs=surf(fxlong(i,j),fxlat(i,j),zsf(i,j),fgrnhfx(i,j),'EdgeColor','none','FaceAlpha',0.7);
  38 % caxis([0,1e6]);   % for heat flux color
  39 axis tight, colorbar
  40 hold on
  41 [c,hc]=contour3(fxlong(i,j),fxlat(i,j),zsf(i,j));
  42 for h=hc(:)', set(h,'EdgeColor','black');end  % color the contours black
  43
  44
  45 % surface wind arrows
  46 if swind,
  47     ii=dmin(1):qstep(1):dmax(1);jj=dmin(2):qstep(2):dmax(2);
  48     hq=quiver3(fxlong(ii,jj),fxlat(ii,jj),zsf(ii,jj),...
  49         xscale*uf(ii,jj),yscale*vf(ii,jj),zeros(size(vf(ii,jj))),qs);
  50     set(hq,'Color','black')
  51 end
  52
  53 % geopotential altitude of the centers of lower cell face
  54 % ph is perturbation from the background phb
  55 % note ph and phb are staggered just as w
  56 a = (ph+phb)/9.81;
  57 % check if the lowest background geopotential height is terrain height
  58 err_hgt = big(hgt-phb(:,:,1)/9.81)
  59 % interpolate the geopotential and w velocity from staggered grid to cell centers
  60 ac = (a(:,:,1:end-1)+a(:,:,2:end))/2;
  61 wc = (w(:,:,1:end-1)+w(:,:,2:end))/2;
  62 % interpolate u and v components of wind velocity from staggered grid to cell centers
  63 uc = (u(1:end-1,:,:)+u(2:end,:,:))/2;  % from front and rear face to the center
  64 vc = (v(:,1:end-1,:)+v(:,2:end,:))/2;  % from front and rear face to the center
  65
  66 % atmosphere wind arrows
  67 ia=amin(1):astep(1):amax(1);
  68 ja=amin(2):astep(2):amax(2);
  69 colors={'red','green','blue','black'};
  70 for k=amin(3):astep(3):amax(3),
  71     q1=u(ia,ja,k)*xscale;
  72     q2=v(ia,ja,k)*yscale;
  73     q3=w(ia,ja,k);
  74     hqa=quiver3(axlong(ia,ja),axlat(ia,ja),ac(ia,ja,k),q1,q2,q3,qs);
  75     set(gca,'DataAspectRatio',[xscale yscale 1]);
  76     i=1+mod(k-1,length(colors));
  77     set(hqa,'Color',colors{i});
  78 end
  79
  80 if ideal
  81     xlabel('x (m)')
  82     ylabel('y (m)')
  83 else
  84     % set(gca,'PlotBoxAspectRatioMode','auto');
  85     xlabel('longitude (deg)')
  86     ylabel('latitude (deg)')
  87 end
  88 zlabel('z (m)')
  89 set(gca,'DataAspectRatio',[xscale yscale 1]);
  90 hold off
  91 drawnow
  92
  93 end