clear all
clf
QNOX = 2.9; % source pollutant emission rate, in g/s
QCO = 13.7; % source pollutant emission rate, in g/
u = 2; % horizontal wind velocity along the plume centerline, m/s 
H = 300;
%run gaussian Plume
C=gaussianPlume(QCO,u,H,'reflection',false);
dimc=size(C);
cxz=squeeze(C(1,:,:));
%Building downwash parameter setup
Lb=22;
Hb=55;
Wb=40;
BS=min(Hb,Wb);
BL=max(Hb,Wb);
R=BS^(2/3)*BL^(1/3);
xb=40;
LR=round(1.8*Wb/((Lb/Hb)^0.3*(1+0.24*Wb/Hb)));
HR=Hb+0.22*R;
Hc=zeros(1,dimc(2));
%Near Wake boundry caculation
for i = xb:(xb+round(R/2+1))
    Hc(i)=HR+4*((i-xb)-R/2)^2*(Hb-HR)/R^2;%The rooftop cavity doesnt attach to the building
end
for i = (xb+round(R/2+1)):(LR+xb)
    Hc(i)=HR*(1-((i-xb)-R/2)^2/(Lb+LR-R)^2)^0.5;
    if Hc(i)<=0
        break
    end
end
Hc=real(Hc);
%Far wake boundry caculation
for i = xb:300
    Hw(i)=1.2*R*((i-xb)/R+(Hb/1.2/R)^3)^(1/3);
end

%Near/far wake concentration
    %setting parameters
B=4;
for i = 1:dimc(2)
    Q(i)=sum(cxz(i,:));
end
yc=Wb/(2*pi)^0.5;
UH=u;
    %caculation of Hp
for i = 1:dimc(2)
    for j = 1:dimc(3)
    if cxz(i,j)>1.0e-07
        Hp(i)=j;
        break
    end
    end 
    for j = 1:dimc(3)
    if cxz(i,(dimc(3)+1-j))>1.0e-07
        Hpt(i)=dimc(3)-j+1;
        break
    end
    end
    
end
    %caculation of vertical dispersion coefficient
for i = 2:dimc(2)
    sz(i)=std(cxz(i,Hp(i):Hpt(i)));
    dif=abs(cxz(i,Hp(i):Hpt(i))-sz(i));
    [cmax index]=max(cxz(i,Hp(i):Hpt(i)));
    [cmax index2]=min(dif);
    sz(i)=abs(index-index2);
end
    %caculation of fraction z
fz= erf((Hc-Hp)./(2^0.5)./sz);
for i = 1:dimc(2)
    if fz(i) <0
        fz(i)=0;
    end
end
    %caculation of Cn, finally!
CN=zeros(dimc(2),dimc(3));CF=zeros(dimc(2),dimc(3));
for i = xb:round(xb+1.15*LR)
    for j =1:Hc(i)
        CN(i,j)=B*fz(i)*Q(i)*exp(-0.5*(20/yc)^2)/(UH*Hc(i)*Wb);
        
    end
end
for i = xb:300
    for j =1:Hw(i)
    CF(i,j)=0.3*Q(i)*exp(-0.5*(20/yc)^2)/(pi*UH*Wb);
    end
end
cxz2=cxz;
for i =xb:300
    shift=round(Hw(i)-55);
    cxz2(i,:)=0;
    cxz2(i,round((Hp(i)+shift)):141)=cxz(i,round(Hp(i)):(141-shift));
end
lamda=1/24;
for i = xb:round(xb+Lb+0.85*LR)
    for j = 1:dimc(3)
        cxz(i,j)=CN(i,j)*100+cxz2(i,j);
    end
end
for i = round(xb+Lb+0.85*LR):round(xb+Lb+LR)
    for j = 1:dimc(3)
        cxz(i,j)=cxz2(i,j)+CN(i,j)*100*(1-(i-xb-Lb-0.85*LR)*lamda)+CF(i,j)*80*(i-xb-Lb-0.85*LR)*lamda;
    end
end
for i = round(xb+Lb+LR):round(xb+Lb+1.15*LR)
    for j = 1:dimc(3)
        cxz(i,j)=cxz2(i,j)+CN(i,j)*100*(1-(i-xb-Lb-0.85*LR)*lamda)+CF(i,j)*80*(i-xb-Lb-0.85*LR)*lamda;
    end
end
for i = round(xb+Lb+1.15*LR+1):dimc(2)
    for j = 1:dimc(3)
        cxz(i,j)=cxz2(i,j)+CF(i,j)*80;
    end
end

% end of caculation for near wake concentration
cxz3=cxz';cxz4=zeros(141,302);
f=1;
for i =2:302
pcolor(cxz4)
hold on
pcolor(cxz3(:,1:i));
axis([0 302 1 141]);
shading interp;
rectangle('Position',[xb 0 Lb Hb],'FaceColor', 'k');
%end
%plot wake field
sizeHc=size(Hc);
plot([xb:sizeHc(2)],Hc(xb:sizeHc(2)));
plot([xb:300],Hw(xb:300));
if mod(i,10)==0
    saveas(gcf,['downwash_',num2str(i),'.png'])
end
hold off
end


