#include
#include
#include
#include


#define?R??????????????????????????????????18.0//initial（產生半圓液滴）
#define?X0?????????????????????????????????1
#define?Y0?????????????????????????????????（LY-30）

#define?tau0???????????????????????????????1.0//流體性質參數
#define?tau1???????????????????????????????1.0
#define?m0?????????????????????????????????1.0
#define?m1?????????????????????????????????1.0
#define?G??????????????????????????????????0.2//作用力參數
#define?GW0????????????????????????????????-0.01
#define?GW1????????????????????????????????0.01
#define?GR?????????????????????????????????-0.0004

//define
#define?LX????????????????????????????????32
#define?LY????????????????????????????????302
#define?Q?????????????????????????????????9

#define?fluid??????????????????????????????0//initial_geo
#define?solid??????????????????????????????1



#define?n0_in??????????????????????????????2.0
#define?n0_out?????????????????????????????0.01
#define?n1_in??????????????????????????????0.01
#define?n1_out?????????????????????????????2.0



#define?PI?????????????????????????????????3.1415926

#define?ERR????????????????????????????????1.0e-6//收斂判斷
#define?IT?????????????????????????????????20000
#define?max（xy）	??????????????????????（（（x）?<=?（y））???（y）?:?（x））

int?flag[LX][LY];
double?s[LX][LY]?e[Q][2]
????????f0[LX][LY][Q]f1[LX][LY][Q]
????????n0[LX][LY]n1[LX][LY]n[LX][LY]den0[LX][LY]den1[LX][LY]den[LX][LY]
		veq0x[LX][LY]veq0y[LX][LY]veq1x[LX][LY]veq1y[LX][LY]vx[LX][LY]vy[LX][LY]//平衡態速度
		g0[LX][LY][Q]g1[LX][LY][Q]
		F0x[LX][LY]F0y[LX][LY]F1x[LX][LY]F1y[LX][LY]
	????pressure[LX][LY]p[LX][LY]
		C[LX][LY];
double?rfa[2]gravity0gravity1;

double?eq（int?indexint?kdouble?ndouble?uxdouble?uy）;
void?ini_geo（void）;
void?initial（void）;
void?collision（void）;
void?stream（void）;
void?para（void）;
void?evolve（void）;
void?contact_angle（int?Tdouble?errint?N）;
//double?err（double?px[LX][LY]double?py[LX][LY]）;
double?err（double?p[LX][LY]）;
void?write（int?T）;


double?eq（int?indexint?kdouble?ndouble?uxdouble?uy）
{
	double?rf;
	rf=rfa[index];
	double?euuueq;
	eu=ux*e[k][0]+uy*e[k][1];
	uu=ux*ux+uy*uy;
	if（k==0）
		eq=rf*n-2.0/3*n*uu;
	else?if（k>=1&&k<=4）
		eq=（1-rf）*n/5+1.0/3*n*eu+1.0/2*n*eu*eu-1.0/6*n*uu;
	else
		eq=（1-rf）*n/20+1.0/12*n*eu+1.0/8*n*eu*eu-1.0/24*n*uu;

	return?eq;
}
//ini_geo
void?ini_geo（）
{
	int?ij;
	for（j=0;j		for（i=0;i		{
			if（i==0||i==LX-1）
			{
				flag[i][j]=solid;
				s[i][j]=1.0;
			}
			else
			{
				flag[i][j]=fluid;
				s[i][j]=0.0;
			}
		}
}
//initial
void?initial（）
{
	int?ijk;
	for（j=0;j		for（i=0;i		{
			if（i==0||i==LX-1）
				n0[i][j]=n1[i][j]=0;
			else
			{
				if（（double（（i-X0）*（i-X0）+（j-Y0）*（j-Y0）））<=R*R）
				{
					n0[i][j]=n0_in;
					n1[i][j]=n1_in;
				}
				else
				{
					n0[i][j]=n0