disp('zeng');
z=10;
kar=.4;
g=9.8;
zx(1)=1e-4;
j=0;
al1=0.011;
al2=0.016;
umax=20;


for l=.4:.1:20;
  	j=j+1;
   u(j)=l;
   zx(j)=1e-4;
   f=((l-10)/5)^2;
   fg=1+f;
   if l<10,
      al=al1;
   else
      al=al1/fg+al2*f/fg;
   end;
   alj(j)=al;
   cdy(j)=0.29+3.1/u(j)+7.7/u(j)^2;
   if u(j)>6,
      cdy(j)=.6+.07*u(j);
   end;
   cdy(j)=0.001*cdy(j);
	for i=1:10;
   	c=kar/log(z/zx(j));
      us=u(j)*c;
      zx(j)=al*us*us/g+.11*1.5e-5/us;
   end;   
   ust(j)=us;
   rr(j)=us*zx(j)/1.5e-5;
   zq1(j)=5e-5/rr(j)^.6;
   zq2(j)=zx(j)/exp(2.67*rr(j)^.25-2.57);
   rq1(j)=zq1(j)*ust(j)/1.5e-5;
   rq2(j)=zq2(j)*ust(j)/1.5e-5;
   cq1(j)=0.16/log(z/zx(j))/log(z/zq1(j));
   cq2(j)=0.16/log(z/zx(j))/log(z/zq2(j));
   cdn(j)=0.16/log(z/zx(j))/log(z/zx(j));
end;
figure;
loglog(rr(1:j),zq1(1:j),rr(1:j), zq2(1:j),'o');
title('line-cwf    o-zeng');
xlabel('R_r');
ylabel('z_q (m)');

figure;
loglog(rr(1:j),rq1(1:j),rr(1:j), rq2(1:j),'o');
title('line-cwf    o-zeng');
xlabel('R_r');
ylabel('R_q');

figure;
plot(u(1:j),cq1(1:j),u(1:j),cq2(1:j),'o');
title('line-cwf    o-zeng');
xlabel('U_{10} (m/s)');
ylabel('C_{En}');
axis([0 20 0.001 0.0015]);


figure;
plot(u(1:j),cdn(1:j),u(1:j),cq2(1:j),'o', u(1:j),cdy(1:j),'.');
title('line-cwf    ');
xlabel('U_{10} (m/s)');
ylabel('C_{dn}');
axis([0 20 0.0008 0.0020]);
