clear
delt=10.0/16.;
delt2=1./32.;
% mean building height

h1=10.0;
d1=h1;
H1=100.0;
hm=h1;
%hm=17.2;

h2=1.0;
d2=h2;
H2=8.0;

utau=0.653;
ustar=sqrt(1-d1/H1)*utau

un=utau;
%un=ustar;
un2=un*un;

utau2=1.0;
ustar2=sqrt(1-d2/H2)*utau2

u2n=utau2;
%u2n=ustar2
u2n2=u2n*u2n;

[k1,u1,u2r1,v2r1,w2r1,tke1,uwd1,u2d1,v2d1,w2d1,dke1,sigu1,sigtke1,siguwd1,sigdke1]=textread('1dspav.data','%d %f %f %f %f %f %f %f %f %f %f %f %f %f %f');

%z1=(k1-1)*delt;
z1=k1*delt;

mke1=0.5*u1.*u1;
totke1=mke1+dke1+tke1;

[k2,u2,uwr2,u2r2,v2r2,w2r2,tke2,uwd2,u2d2,v2d2,w2d2,dke2,sigu2,siguwr2,sigtke2,siguwd2,sigdke2]=textread('1dspav_16128.data','%d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f');

z2=(k2-0.5)*delt2;

mke2=0.5*u2.*u2;
totke2=mke2+dke2+tke2;

subplot(1,2,1)
    hold on
     plot(mke1/un2,z1/hm,'k--','LineWidth',2)
     plot(dke1/un2,z1/hm,'k-.','LineWidth',2)
     plot(tke1/un2,z1/hm,'k:','LineWidth',2)
     plot(totke1/un2,z1/hm,'k','LineWidth',2)

     xlabel('MKE/u_\tau^2, DKE/u_\tau^2, TKE/u_\tau^2, TOTKE/u_\tau^2')
     ylabel('z/h_m')
     %     ylabel('(z-d)/H')
     legend('MKE','DKE','TKE','TOTKE');
     title('(a) Random array');
     axis([0,4,0,10])
     plot([0 4],[1 1],'k--')
     plot([0 4],[1.72 1.72],'k--')

subplot(1,2,2)
    hold on
     plot(mke2/u2n2,z2,'k--','LineWidth',2)
     plot(dke2/u2n2,z2,'k-.','LineWidth',2)
     plot(tke2/u2n2,z2,'k:','LineWidth',2)
     plot(totke2/u2n2,z2,'k','LineWidth',2)

     xlabel('MKE/u_\tau^2, DKE/u_\tau^2, TKE/u_\tau^2, TOTKE/u_\tau^2')
     ylabel('z/h')
     %     ylabel('(z-d)/H')
   %  legend('MKE','DKE','TKE','TOTKE');
     title('(b) Uniform array');
     axis([0,4,0,10])
     plot([0 4],[1 1],'k--')
%     plot([0 4],[1.72 1.72],'k--')

     print('-depsc','-r400','spav_comp_engy.eps')