%clf;
clear;
clc;

%% Setting up parameters for Gaussian air pollutant dispersion equation

QNOX = 2.9; % source pollutant emission rate, in g/s
QCO = 13.7; % source pollutant emission rate, in g/
u = 3; % horizontal wind velocity along the plume centerline, m/s 
H = 300; % height of emission plume centerline above ground level, in m
format long

%% Run gaussian dispersion model
C=gaussianPlume(QCO,u,H);

%% Get the color graph
cubedim=size(C);
len_x=size(C,2);
len_y=size(C,1);
len_z=size(C,3);
cubexyz=len_x*len_y*len_z;
cubexy=len_x*len_y;
cubexz=len_x*len_z;

x=0:100:100*(len_x-1);
y=0:100:100*(len_y-1);
z=0:100:100*(len_z-1);

yy=repmat(y,[1,cubexz]);
zz=reshape(repmat(z,[cubexy,1]),[1,cubexyz]);
xx=repmat(reshape(repmat(x,[len_y,1]),[1,cubexy]),[1,len_z]);
%%
Cvector=reshape(C,[1,cubexyz]);
Bins = quantile(Cvector,100);
[num, ind]= histc(Cvector,Bins);
cc = jet(100);
% Ccolor=zeros(length(ind),3);
% for i=1:length(ind)
%     Ccolor(i,:) = cc(ind(i)+1,:);
% end
% scatter3(xx,yy,zz,5,Ccolor,'filled', 'MarkerEdgeColor','none');
% alpha(.001)

%% initial condition
x=50; % 100 meter
y=40;
z=10;
noparticle=10000;
time=70;
dt=1;
% initial position
p_x=zeros(noparticle,1);
p_y=5*rand(noparticle,1)+2000;
p_z=ones(noparticle,1)+300;
% intial velocity
u=zeros(noparticle,1);
v=zeros(noparticle,1);
w=10*ones(noparticle,1);
% velocity & position update
pos_x=p_x;
pos_y=p_y;
pos_z=p_z;
ut=zeros(noparticle,1);
vt=zeros(noparticle,1);
wt=w;

pout=100;
f=0;
% renewing postion, velocity and plotting
for iter = 1:time
    [ut(1:pout), vt(1:pout), wt(1:pout)] = changev(p_x(1:pout),p_y(1:pout),p_z(1:pout),u(1:pout),v(1:pout),w(1:pout));
    [pos_x(1:pout), pos_y(1:pout), pos_z(1:pout)] = changep(p_x(1:pout),p_y(1:pout),p_z(1:pout),ut(1:pout),vt(1:pout),wt(1:pout),dt);
    p_x(1:pout)=pos_x(1:pout);
    p_y(1:pout)=pos_y(1:pout);
    p_z(1:pout)=pos_z(1:pout);
    u(1:pout)=ut(1:pout);
    v(1:pout)=vt(1:pout);
    w(1:pout)=wt(1:pout);
    % color
    zone=zeros(pout,3);
    color=zeros(pout,3);
    for i=1:pout
        zone(i,:)=[round(p_x(i)/100),round(p_y(i)/100),round(p_z(i)/100)];
        color(i,:)= cc(ind(zone(i,3)*x*y+zone(i,3)*40+zone(i,3)),:);
    end
    % plot
    if mod(iter,1)==0
        f=f+1;
        clf;
        scatter3(pos_x(1:pout),pos_y(1:pout),pos_z(1:pout),10,color,'filled', 'MarkerEdgeColor','none');
        view([0 90])
        axis([0 5000 0 4000 0 1000])
        grid on;
        pause(0.001)
        drawnow     
        saveas(gcf,['plum2_',sprintf('%03d',f),'.png']);
        hold off
    end
    if pout < noparticle
        pout = pout+100;
    end
end
