Data;
while(~StopFlag)
    Cur_Iter=Cur_Iter+1;
    %宏觀量計(jì)算
    rho=sum(F,3);                         %我的理解：得到每個(gè)格點(diǎn)的局部密度
    
    if Cur_Iter >1                        %有什么必要么？
        ux=zeros(LN,WN);
        uy=zeros(LN,WN);
        for ic=1:N_c-1;                   %求出了宏觀動(dòng)量
            ux = ux + C_x(ic).*F(:,:,ic);
            uy = uy + C_y(ic).*F(:,:,ic);
        end
    end
    ux(ijWhole)=ux(ijWhole)./rho(ijWhole);
    uy(ijWhole)=uy(ijWhole)./rho(ijWhole);%求出了宏觀速度
    
    uxsq(ijWhole)=ux(ijWhole).^2;
    uysq(ijWhole)=uy(ijWhole).^2; 
    usq(ijWhole)=uxsq(ijWhole)+uysq(ijWhole);%求出了個(gè)格點(diǎn)宏觀速度的平方，為計(jì)算平衡分布函數(shù)做準(zhǔn)備
    %***************************************
    %          粒子碰撞（整個(gè)區(qū)域）
    %***************************************
    %方向信息
    %  6 2 5        NW  N  NE
    %  3 9 1 ---->  W   RP  E
    %  7 4 8        SW  S  SE
    East=1;North=2;West=3;South=4;NE=5;NW=6;SW=7;SE=8;RP=9;
    rt0=w0.*rho;rt1=w1.*rho;rt2=w2.*rho;
    F_EQ(ijWhole+NxM*(1-1))=rt1(ijWhole).*(1+f1*ux(ijWhole)+f2*uxsq(ijWhole)-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(9-1))=rt0(ijWhole).*(1-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(2-1))=rt1(ijWhole).*(1+f1*uy(ijWhole)+f2*uysq(ijWhole)-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(3-1))=rt1(ijWhole).*(1-f1*ux(ijWhole)+f2*uxsq(ijWhole)-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(4-1))=rt1(ijWhole).*(1-f1*uy(ijWhole)+f2*uysq(ijWhole)-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(5-1))=rt2(ijWhole).*(1+f1*(+ux(ijWhole)+uy(ijWhole))+f2*(+ux(ijWhole)+uy(ijWhole)).^2-f3.*usq(ijWhole));
    F_EQ(ijWhole+NxM*(6-1))=rt2(ijWhole).*(1+f1*(-ux(ijWhole)+uy(ijWhole))+f2*(-ux(ijWhole)+uy(ijWhole)).^2-f3.*usq(ijWhole));
    F_EQ(ijWhole+NxM*(7-1))=rt2(ijWhole).*(1+f1*(-ux(ijWhole)-uy(ijWhole))+f2*(-ux(ijWhole)-uy(ijWhole)).^2-f3.*usq(ijWhole));
    F_EQ(ijWhole+NxM*(8-1))=rt2(ijWhole).*(1+f1*(+ux(ijWhole)-uy(ijWhole))+f2*(+ux(ijWhole)-uy(ijWhole)).^2-f3.*usq(ijWhole)); 
    F=(1.-omega).*F+omega.*F_EQ;%注意omega的具體意義
    %在流體上施加壓力
    for ic=1:N_c                              %外層為ic方向循環(huán)，內(nèi)層為格子點(diǎn)循環(huán)
        for ia=1:lenWhole
            i=iWhole(ia);  
            j=jWhole(ia);
            F(i,j,ic)= F(i,j,ic) + force(ic); %force已經(jīng)在data.m中定義好了，注意這是常力矩陣，我的需要在這個(gè)地方做文章！
        end
    end
    %***************************************
    %               粒子遷移
    %***************************************   
    F_EQ = F;
    for ic=1:1:N_c-1       %碰撞前的速度向量
        ic2=ic_op(ic);     %碰撞后的速度向量，利用前一個(gè)當(dāng)索引
        temp1=F_EQ(:,:,ic);%存儲(chǔ)了每個(gè)方向上所有點(diǎn)的平衡分布函數(shù)
        %內(nèi)部區(qū)域
        for ia=1:1:lenInner
            i=iInner(ia);
            j=jInner(ia);
            i2=i+C_y(ic);        %由于det_x=1,det_t=1,所以i,j既是索引，又可以代表空間的位置坐標(biāo)
            j2=j+C_x(ic);
            i2=yi2(i2+1);        %周期性邊界條件，巧妙，服了！但是將inlet,outlet設(shè)成兩層的意義是什么？
            F(i2,j2,ic)=temp1(i,j);%將該方向上的分布函數(shù)遷移到運(yùn)動(dòng)后的位置
        end
        %邊界                    %只有邊界層的格子才有必要驗(yàn)證是否用反彈邊界條件
        for ia=1:1:lenObs
            i=iObs(ia);
            j=jObs(ia);
            i2=i+C_y(ic);
            j2=j+C_x(ic);        
            i2=yi2(i2+1);        %豎直方向仍然要考慮周期條件
            if Channel2D(i2,j2)==0
                F(i,j,ic2) =temp1(i,j);%如果碰到管壁反彈回去，反向增加分布函數(shù)
            else
                F(i2,j2,ic)=temp1(i,j);
            end
        end
    end 
    
    
    %重新計(jì)算宏觀速度和密度
    rho=sum(F,3);
    uyout= zeros(LN,WN);
    for ic=1:N_c-1;
        uyout= uyout + C_y(ic).*F(:,:,ic);
    end
    %收斂條件計(jì)算
    if mod(Cur_Iter,Check_Iter)==0
        vect=rho(ijWhole); 
        cur_density=mean(vect);
        vect=uy(ijWhole);
        av_vel_int=mean(vect);
        av_vel_tp1=av_vel_int; 
        Condition=abs(abs(av_vel_t/av_vel_tp1 )-1);%收斂系數(shù)
        Cond_path=[Cond_path,Condition];%記錄收斂系數(shù)
        density_path=[density_path,cur_density];%記錄密度
        av_vel_t=av_vel_tp1;
        if (Condition < toler)||(Cur_Iter > Max_Iter)
            StopFlag=true;
            display('Stop iteration: Convergence met or iteration exeeding the max allowed')
            display(['Current iteration: ',num2str(Cur_Iter),' Max Number of iter: ',num2str(Max_Iter)])
            break;
        end
    end
    %輸出
    if mod(Cur_Iter,Output_Every)==0
        rho=sum(F,3);
        up=2; %從較低層可視化
        figure(15);hold off;feather(ux(LN-up,:),uy(LN-up,:));
        figure(15);hold on;plot(uy_analy_profile,'r-');
        title('Analytical vs LB calculated, fluid velocity parabolic profile');
        pause(3);
    end
end 
figure;plot(Cond_path(2:end));title('convergence path');
figure;plot(uy(LN-up,:)-uy_analy_profile);title('difference : LB - Analytical solution');
toc