function [lbm,part,sim] = LBM_main_module (lbm,part,sim)
%% parameters
ly = lbm.ly;
lx = lbm.lx;
rho1 = lbm.rho1;
A=lbm.A_pre;
tao = lbm.tao;
B = (sim.mesh_solid_fra.*(tao-0.5))./((1-sim.mesh_solid_fra)+(tao-0.5));
B_1 = 1-B;
uTotX2 = lbm.ux;
uTotY2 = lbm.uy;
%% 
uTotX2Us = sim.mesh_v_conv(:,:,2);
uTotY2Us = sim.mesh_v_conv(:,:,1);
tempNS1        = 3*(lbm.cxNS(1)*uTotX2+lbm.cyNS(1)*uTotY2);
tempNS2        = 3*(lbm.cxNS(2)*uTotX2+lbm.cyNS(2)*uTotY2);
tempNS3        = 3*(lbm.cxNS(3)*uTotX2+lbm.cyNS(3)*uTotY2);
tempNS4        = 3*(lbm.cxNS(4)*uTotX2+lbm.cyNS(4)*uTotY2);
tempNS5        = 3*(lbm.cxNS(5)*uTotX2+lbm.cyNS(5)*uTotY2);
tempNS6        = 3*(lbm.cxNS(6)*uTotX2+lbm.cyNS(6)*uTotY2);
tempNS7        = 3*(lbm.cxNS(7)*uTotX2+lbm.cyNS(7)*uTotY2);
tempNS8        = 3*(lbm.cxNS(8)*uTotX2+lbm.cyNS(8)*uTotY2);
tempNS9        = 3*(lbm.cxNS(9)*uTotX2+lbm.cyNS(9)*uTotY2);
temp1 = 1.5*(uTotX2.^2+uTotY2.^2);
lbm.fEq1(:,:)   = rho1.*lbm.tNS(1).*( 1 + tempNS1 + 0.5*(tempNS1.*tempNS1) - temp1 );
lbm.fEq2(:,:)   = rho1.*lbm.tNS(2).*( 1 + tempNS2 + 0.5*(tempNS2.*tempNS2) - temp1 );
lbm.fEq3(:,:)   = rho1.*lbm.tNS(3).*( 1 + tempNS3 + 0.5*(tempNS3.*tempNS3) - temp1 );
lbm.fEq4(:,:)   = rho1.*lbm.tNS(4).*( 1 + tempNS4 + 0.5*(tempNS4.*tempNS4) - temp1 );
lbm.fEq5(:,:)   = rho1.*lbm.tNS(5).*( 1 + tempNS5 + 0.5*(tempNS5.*tempNS5) - temp1 );
lbm.fEq6(:,:)   = rho1.*lbm.tNS(6).*( 1 + tempNS6 + 0.5*(tempNS6.*tempNS6) - temp1 );
lbm.fEq7(:,:)   = rho1.*lbm.tNS(7).*( 1 + tempNS7 + 0.5*(tempNS7.*tempNS7) - temp1 );
lbm.fEq8(:,:)   = rho1.*lbm.tNS(8).*( 1 + tempNS8 + 0.5*(tempNS8.*tempNS8) - temp1 );
lbm.fEq9(:,:)   = rho1.*lbm.tNS(9).*( 1 + tempNS9 + 0.5*(tempNS9.*tempNS9) - temp1 );
cuNSUs1        = 3*(lbm.cxNS(1)*uTotX2Us+lbm.cyNS(1)*uTotY2Us);
cuNSUs2        = 3*(lbm.cxNS(2)*uTotX2Us+lbm.cyNS(2)*uTotY2Us);
cuNSUs3        = 3*(lbm.cxNS(3)*uTotX2Us+lbm.cyNS(3)*uTotY2Us);
cuNSUs4        = 3*(lbm.cxNS(4)*uTotX2Us+lbm.cyNS(4)*uTotY2Us);
cuNSUs5        = 3*(lbm.cxNS(5)*uTotX2Us+lbm.cyNS(5)*uTotY2Us);
cuNSUs6        = 3*(lbm.cxNS(6)*uTotX2Us+lbm.cyNS(6)*uTotY2Us);
cuNSUs7        = 3*(lbm.cxNS(7)*uTotX2Us+lbm.cyNS(7)*uTotY2Us);
cuNSUs8        = 3*(lbm.cxNS(8)*uTotX2Us+lbm.cyNS(8)*uTotY2Us);
cuNSUs9        = 3*(lbm.cxNS(9)*uTotX2Us+lbm.cyNS(9)*uTotY2Us);
temp1 = 1.5*(uTotX2Us.^2+uTotY2Us.^2);
fEq1Us(:,:)   = rho1.*lbm.tNS(1).*( 1 + cuNSUs1 + 0.5*(cuNSUs1.*cuNSUs1) - temp1 );
fEq2Us(:,:)   = rho1.*lbm.tNS(2).*( 1 + cuNSUs2 + 0.5*(cuNSUs2.*cuNSUs2) - temp1 );
fEq3Us(:,:)   = rho1.*lbm.tNS(3).*( 1 + cuNSUs3 + 0.5*(cuNSUs3.*cuNSUs3) - temp1 );
fEq4Us(:,:)   = rho1.*lbm.tNS(4).*( 1 + cuNSUs4 + 0.5*(cuNSUs4.*cuNSUs4) - temp1 );
fEq5Us(:,:)   = rho1.*lbm.tNS(5).*( 1 + cuNSUs5 + 0.5*(cuNSUs5.*cuNSUs5) - temp1 );
fEq6Us(:,:)   = rho1.*lbm.tNS(6).*( 1 + cuNSUs6 + 0.5*(cuNSUs6.*cuNSUs6) - temp1 );
fEq7Us(:,:)   = rho1.*lbm.tNS(7).*( 1 + cuNSUs7 + 0.5*(cuNSUs7.*cuNSUs7) - temp1 );
fEq8Us(:,:)   = rho1.*lbm.tNS(8).*( 1 + cuNSUs8 + 0.5*(cuNSUs8.*cuNSUs8) - temp1 );
fEq9Us(:,:)   = rho1.*lbm.tNS(9).*( 1 + cuNSUs9 + 0.5*(cuNSUs9.*cuNSUs9) - temp1 );
colli1 = lbm.fIn1-lbm.fIn1+fEq1Us-lbm.fEq1;
colli2 = lbm.fIn4-lbm.fIn2+fEq2Us-lbm.fEq4;
colli3 = lbm.fIn5-lbm.fIn3+fEq3Us-lbm.fEq5;
colli4 = lbm.fIn2-lbm.fIn4+fEq4Us-lbm.fEq2;
colli5 = lbm.fIn3-lbm.fIn5+fEq5Us-lbm.fEq3;
colli6 = lbm.fIn8-lbm.fIn6+fEq6Us-lbm.fEq8;
colli7 = lbm.fIn9-lbm.fIn7+fEq7Us-lbm.fEq9;
colli8 = lbm.fIn6-lbm.fIn8+fEq8Us-lbm.fEq6;
colli9 = lbm.fIn7-lbm.fIn9+fEq9Us-lbm.fEq7;
f_partx = B.*(lbm.cxNS(1)*colli1+lbm.cxNS(2)*colli2+lbm.cxNS(3)*colli3+lbm.cxNS(4)*colli4+lbm.cxNS(5)*colli5+...
    lbm.cxNS(6)*colli6+lbm.cxNS(7)*colli7+lbm.cxNS(8)*colli8+lbm.cxNS(9)*colli9);
f_party = B.*(lbm.cyNS(1)*colli1+lbm.cyNS(2)*colli2+lbm.cyNS(3)*colli3+lbm.cyNS(4)*colli4+lbm.cyNS(5)*colli5+...
    lbm.cyNS(6)*colli6+lbm.cyNS(7)*colli7+lbm.cyNS(8)*colli8+lbm.cyNS(9)*colli9);
%% collision step
fOut1  = lbm.fIn1 - A.*B_1.*(lbm.fIn1-lbm.fEq1)+B.*colli1;%+vForce1;
fOut2  = lbm.fIn2 - A.*B_1.*(lbm.fIn2-lbm.fEq2)+B.*colli2;%+vForce2;
fOut3  = lbm.fIn3 - A.*B_1.*(lbm.fIn3-lbm.fEq3)+B.*colli3;%+vForce3;
fOut4  = lbm.fIn4 - A.*B_1.*(lbm.fIn4-lbm.fEq4)+B.*colli4;%+vForce4;
fOut5  = lbm.fIn5 - A.*B_1.*(lbm.fIn5-lbm.fEq5)+B.*colli5;%+vForce5;
fOut6  = lbm.fIn6 - A.*B_1.*(lbm.fIn6-lbm.fEq6)+B.*colli6;%+vForce6;
fOut7  = lbm.fIn7 - A.*B_1.*(lbm.fIn7-lbm.fEq7)+B.*colli7;%+vForce7;
fOut8  = lbm.fIn8 - A.*B_1.*(lbm.fIn8-lbm.fEq8)+B.*colli8;%+vForce8;
fOut9  = lbm.fIn9 - A.*B_1.*(lbm.fIn9-lbm.fEq9)+B.*colli9;%+vForce9;
%% STREAMING STEP
lbm.fIn1(:,:) = fOut1(:,:);
lbm.fIn2(:,:) = circshift(fOut2(:,:), [-lbm.cxNS(2),lbm.cyNS(2)]);
lbm.fIn6(:,:) = circshift(fOut6(:,:), [-lbm.cxNS(6),lbm.cyNS(6)]);
lbm.fIn9(:,:) = circshift(fOut9(:,:), [-lbm.cxNS(9),lbm.cyNS(9)]);
lbm.fIn4(:,:) = circshift(fOut4(:,:), [-lbm.cxNS(4),lbm.cyNS(4)]);
lbm.fIn7(:,:) = circshift(fOut7(:,:), [-lbm.cxNS(7),lbm.cyNS(7)]);
lbm.fIn8(:,:) = circshift(fOut8(:,:), [-lbm.cxNS(8),lbm.cyNS(8)]);
lbm.fIn3(:,:) = circshift(fOut3(:,:), [-lbm.cxNS(3),lbm.cyNS(3)]);
lbm.fIn5(:,:) = circshift(fOut5(:,:), [-lbm.cxNS(5),lbm.cyNS(5)]);
%% pressure inlet and outlet
rho_in= lbm.inlet_rho;
tempNS1        = 3*(lbm.cxNS(1)*((uTotX2(2,:)+uTotX2(3,:))/2)+lbm.cyNS(1)*((uTotY2(2,:)+uTotY2(3,:))/2));
tempNS2        = 3*(lbm.cxNS(2)*((uTotX2(2,:)+uTotX2(3,:))/2)+lbm.cyNS(2)*((uTotY2(2,:)+uTotY2(3,:))/2));
tempNS3        = 3*(lbm.cxNS(3)*((uTotX2(2,:)+uTotX2(3,:))/2)+lbm.cyNS(3)*((uTotY2(2,:)+uTotY2(3,:))/2));
tempNS4        = 3*(lbm.cxNS(4)*((uTotX2(2,:)+uTotX2(3,:))/2)+lbm.cyNS(4)*((uTotY2(2,:)+uTotY2(3,:))/2));
tempNS5        = 3*(lbm.cxNS(5)*((uTotX2(2,:)+uTotX2(3,:))/2)+lbm.cyNS(5)*((uTotY2(2,:)+uTotY2(3,:))/2));
tempNS6        = 3*(lbm.cxNS(6)*((uTotX2(2,:)+uTotX2(3,:))/2)+lbm.cyNS(6)*((uTotY2(2,:)+uTotY2(3,:))/2));
tempNS7        = 3*(lbm.cxNS(7)*((uTotX2(2,:)+uTotX2(3,:))/2)+lbm.cyNS(7)*((uTotY2(2,:)+uTotY2(3,:))/2));
tempNS8        = 3*(lbm.cxNS(8)*((uTotX2(2,:)+uTotX2(3,:))/2)+lbm.cyNS(8)*((uTotY2(2,:)+uTotY2(3,:))/2));
tempNS9        = 3*(lbm.cxNS(9)*((uTotX2(2,:)+uTotX2(3,:))/2)+lbm.cyNS(9)*((uTotY2(2,:)+uTotY2(3,:))/2));
temp1 = 1.5*(((uTotX2(2,:)+uTotX2(3,:))/2).^2+((uTotY2(2,:)+uTotY2(3,:))/2).^2);
lbm.fIn1(1,:)   = rho_in.*lbm.tNS(1).*( 1 + tempNS1 + 0.5*(tempNS1.*tempNS1) - temp1 );
lbm.fIn2(1,:)   = rho_in.*lbm.tNS(2).*( 1 + tempNS2 + 0.5*(tempNS2.*tempNS2) - temp1 );
lbm.fIn3(1,:)   = rho_in.*lbm.tNS(3).*( 1 + tempNS3 + 0.5*(tempNS3.*tempNS3) - temp1 );
lbm.fIn4(1,:)   = rho_in.*lbm.tNS(4).*( 1 + tempNS4 + 0.5*(tempNS4.*tempNS4) - temp1 );
lbm.fIn5(1,:)   = rho_in.*lbm.tNS(5).*( 1 + tempNS5 + 0.5*(tempNS5.*tempNS5) - temp1 );
lbm.fIn6(1,:)   = rho_in.*lbm.tNS(6).*( 1 + tempNS6 + 0.5*(tempNS6.*tempNS6) - temp1 );
lbm.fIn7(1,:)   = rho_in.*lbm.tNS(7).*( 1 + tempNS7 + 0.5*(tempNS7.*tempNS7) - temp1 );
lbm.fIn8(1,:)   = rho_in.*lbm.tNS(8).*( 1 + tempNS8 + 0.5*(tempNS8.*tempNS8) - temp1 );
lbm.fIn9(1,:)   = rho_in.*lbm.tNS(9).*( 1 + tempNS9 + 0.5*(tempNS9.*tempNS9) - temp1 );
rho_out= lbm.outlet_rho;
tempNS1        = 3*(lbm.cxNS(1)*uTotX2(lx-1,:)+lbm.cyNS(1)*uTotY2(lx-1,:));
tempNS2        = 3*(lbm.cxNS(2)*uTotX2(lx-1,:)+lbm.cyNS(2)*uTotY2(lx-1,:));
tempNS3        = 3*(lbm.cxNS(3)*uTotX2(lx-1,:)+lbm.cyNS(3)*uTotY2(lx-1,:));
tempNS4        = 3*(lbm.cxNS(4)*uTotX2(lx-1,:)+lbm.cyNS(4)*uTotY2(lx-1,:));
tempNS5        = 3*(lbm.cxNS(5)*uTotX2(lx-1,:)+lbm.cyNS(5)*uTotY2(lx-1,:));
tempNS6        = 3*(lbm.cxNS(6)*uTotX2(lx-1,:)+lbm.cyNS(6)*uTotY2(lx-1,:));
tempNS7        = 3*(lbm.cxNS(7)*uTotX2(lx-1,:)+lbm.cyNS(7)*uTotY2(lx-1,:));
tempNS8        = 3*(lbm.cxNS(8)*uTotX2(lx-1,:)+lbm.cyNS(8)*uTotY2(lx-1,:));
tempNS9        = 3*(lbm.cxNS(9)*uTotX2(lx-1,:)+lbm.cyNS(9)*uTotY2(lx-1,:));
temp1 = 1.5*(uTotX2(lx-1,:).^2+uTotY2(lx-1,:).^2);
lbm.fIn1(lx,:)   = rho_out.*lbm.tNS(1).*( 1 + tempNS1 + 0.5*(tempNS1.*tempNS1) - temp1 );
lbm.fIn2(lx,:)   = rho_out.*lbm.tNS(2).*( 1 + tempNS2 + 0.5*(tempNS2.*tempNS2) - temp1 );
lbm.fIn3(lx,:)   = rho_out.*lbm.tNS(3).*( 1 + tempNS3 + 0.5*(tempNS3.*tempNS3) - temp1 );
lbm.fIn4(lx,:)   = rho_out.*lbm.tNS(4).*( 1 + tempNS4 + 0.5*(tempNS4.*tempNS4) - temp1 );
lbm.fIn5(lx,:)   = rho_out.*lbm.tNS(5).*( 1 + tempNS5 + 0.5*(tempNS5.*tempNS5) - temp1 );
lbm.fIn6(lx,:)   = rho_out.*lbm.tNS(6).*( 1 + tempNS6 + 0.5*(tempNS6.*tempNS6) - temp1 );
lbm.fIn7(lx,:)   = rho_out.*lbm.tNS(7).*( 1 + tempNS7 + 0.5*(tempNS7.*tempNS7) - temp1 );
lbm.fIn8(lx,:)   = rho_out.*lbm.tNS(8).*( 1 + tempNS8 + 0.5*(tempNS8.*tempNS8) - temp1 );
lbm.fIn9(lx,:)   = rho_out.*lbm.tNS(9).*( 1 + tempNS9 + 0.5*(tempNS9.*tempNS9) - temp1 );
%%
rho1 = lbm.fIn1+lbm.fIn2+lbm.fIn3+lbm.fIn4+lbm.fIn5+lbm.fIn6+lbm.fIn7+lbm.fIn8+lbm.fIn9;
jx1  = lbm.cxNS(1)*lbm.fIn1+lbm.cxNS(2)*lbm.fIn2+lbm.cxNS(3)*lbm.fIn3+lbm.cxNS(4)*lbm.fIn4...
    +lbm.cxNS(5)*lbm.fIn5+lbm.cxNS(6)*lbm.fIn6+lbm.cxNS(7)*lbm.fIn7+lbm.cxNS(8)*lbm.fIn8+lbm.cxNS(9)*lbm.fIn9;
jy1  = lbm.cyNS(1)*lbm.fIn1+lbm.cyNS(2)*lbm.fIn2+lbm.cyNS(3)*lbm.fIn3+lbm.cyNS(4)*lbm.fIn4...
    +lbm.cyNS(5)*lbm.fIn5+lbm.cyNS(6)*lbm.fIn6+lbm.cyNS(7)*lbm.fIn7+lbm.cyNS(8)*lbm.fIn8+lbm.cyNS(9)*lbm.fIn9;
uTotX2 = jx1 ./ rho1;
uTotY2 = jy1 ./ rho1;
lbm.ux = uTotX2;
lbm.uy = uTotY2;
lbm.rho1 = rho1;
lbm.p = lbm.rho1/3;
lbm.f_partx = f_partx;
lbm.f_party = f_party;
end
