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;

[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;

subplot(2,2,1)
    hold on
    plot(u2d1/un2,z1/hm,'k','LineWidth',2)
%     plot(u2d1/un2,(z1-d1)/H1,'k','LineWidth',2)
%
     plot(u2d2/u2n2,z2,'k-.','LineWidth',2)
%    plot(u2d2/u2n2,(z2-d2)/H2,'LineWidth',2)

%     xlabel('<u2d>/u_*^2')
     xlabel('\sigma^2_{ud}/u_\tau^2')
     ylabel('z/h_m, z/h')
%     ylabel('(z-d)/H')
   legend('Random','Uniform')
     title('(a)');
     axis([-0.1,2.0,0,2.0])
     plot([-0.1 2.0],[1 1],'k--')
     plot([-0.1 2.0],[1.72 1.72],'k--')
     plot([0 0],[0 2],'k:')


     subplot(2,2,2)
    hold on
        plot(v2d1/un2,z1/hm,'k','LineWidth',2)
        %     plot(v2d1/un2,(z1-d1)/H1,'k--','LineWidth',2)
   %
     plot(v2d2/u2n2,z2,'k-.','LineWidth',2)
     %    plot(v2d2/u2n2,(z2-d2)/H2,'--','LineWidth',2)

     %     xlabel('<v2d>/u_*^2')
     xlabel('\sigma^2_{vd}/u_\tau^2')
     ylabel('z/h_m, z/h')
     %     ylabel('(z-d)/H')
      title('(b)');
      axis([-0.1,0.5,0,2.0])
     plot([-0.1 0.5],[1 1],'k--')
     plot([-0.1 0.5],[1.72 1.72],'k--')
     plot([0 0],[0 2],'k:')


     subplot(2,2,3)
    hold on
     plot(w2d1/un2,z1/hm,'k','LineWidth',2)
     %     plot(w2d1/un2,(z1-d1)/H1,'k-.','LineWidth',2)
%
     plot(w2d2/u2n2,z2,'k-.','LineWidth',2)
     %    plot(w2d2/u2n2,(z2-d2)/H2,'-.','LineWidth',2)

     %     xlabel('<w2d>/u_*^2')
     xlabel('\sigma^2_{wd}/u_\tau^2')
     ylabel('z/h_m, z/h')
     %     ylabel('(z-d)/H')
     title('(c)');
     axis([-0.1,0.5,0,2.0])
     plot([-0.1 0.5],[1 1],'k--')
     plot([-0.1 0.5],[1.72 1.72],'k--')
     plot([0 0],[0 2],'k:')


     subplot(2,2,4)
    hold on
     plot(-uwd1/un2,z1/hm,'k','LineWidth',2)
     %     plot(uwd1/un2,(z1-d1)/H1,'k','LineWidth',2)
%
     plot(-uwd2/u2n2,z2,'k-.','LineWidth',2)
     %    plot(uwd2/u2n2,(z2-d2)/H2,'LineWidth',2)

     %     xlabel('-<uwd>/u_*^2')
     xlabel('')

     ylabel('z/h_m, z/h')
     %     ylabel('(z-d)/H')
     title('(d)');
     axis([-0.1,0.2,0,2.0])
     plot([-0.1 0.2],[1 1],'k--')
     plot([-0.1 0.2],[1.72 1.72],'k--')
     plot([0 0],[0 2],'k:')

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