FmB(middb[1],middf[3],ubar0);
FmS(ubar0,middf[3],1);

FmB(middb[0],middf[4],ubar0);
FmS(ubar0,middf[4],1);

/*out<<"midd f2 和值"<<endl;
for(kn=0;kn<4;kn++)
out<<middf[1][kn]<<"   "<<middf[2][kn]<<"   "<<middf[3][kn]<<"    "<<middf[4][kn]<<endl;*/

Fma(middf[1],middf[2],middf[1],'-');
FmS(middf[1],middf[1],1/di);
Fma(middf[3],middf[4],middf[3],'-');
FmS(middf[3],middf[3],1/dj);

/*out<<"midd f3 和值"<<endl;
for(kn=0;kn<4;kn++)
out<<middf[1][kn]<<"   "<<middf[3][kn]<<endl;*/


Fma(middf[1],middf[3],middf[1],'+');
FmS(middf[1],middf[1],dt*dt*0.5);

Fma(middf[0],middf[1],ubar0,'-');
Fma(ubar[4],ubar0,ubar0,'-');

/*out<<"midd f4 和值"<<endl;
for(kn=0;kn<4;kn++)
out<<middf[1][kn]<<"   "<<middf[0][kn]<<endl;*/
//求rou,u,v,p,e

FmS(ubar0,ubar0,1/jj0);

rout[i][j]=ubar0[0];
if(rout[i][j]<=0||rou[i][j]>=10000)
{
cout<<"worng"<<endl;
cout<<endl<<endl<<i<<"    "<<j<<endl;
break;
}
else 
{
	vut[i][j]=ubar0[1]/ubar0[0];
	vvt[i][j]=ubar0[2]/ubar0[0];
	pt[i][j]=(ubar0[3]-0.5*(ubar0[1]*vut[i][j]+ubar0[2]*vvt[i][j]))*(rg-1);

}
/*cout<<endl<<vut[i][j]<<" \t"<<vvt[i][j]<<" \t"<<rout[i][j]<<" \t"<<pt[i][j]<<endl;

out<<endl<<endl<<i<<"    "<<j<<endl;
out.close();//關閉記錄文件*/

//cin>>nl;
	}

    for(i=2;i<fi;i++)
	for(j=2;j<fj;j++)
{
		vu[i][j]=w*vut[i][j]+(1-w)*vu[i][j];
		vv[i][j]=w*vvt[i][j]+(1-w)*vv[i][j];
		rou[i][j]=w*rout[i][j]+(1-w)*rou[i][j];
		p[i][j]=w*pt[i][j]+(1-w)*p[i][j];
}

//邊界確定	

//外推法估算-特征線矯正
//壁面
for(i=2;i<fi;i++)
{
	
	vu[i][1]=2*vu[i][2]-vu[i][3];
    vv[i][1]=2*vv[i][2]-vv[i][3];
	p[i][1]=2*p[i][2]-p[i][3];
    rou[i][1]=2*rou[i][2]-rou[i][3];
fadao(i,1,visx, visy,nkk);//求向導數
uvp(vu[i][1],vv[i][1],rou[i][1], p[i][1], nkk,ubar0);//邊界處理
vu[i][1]=ubar0[0];
vv[i][1]=ubar0[1];
rou[i][1]=ubar0[2];
p[i][1]=ubar0[3];
    
	vu[i][fj]=2*vu[i][fj-1]-vu[i][fj-2];
    vv[i][fj]=2*vv[i][fj-1]-vv[i][fj-2];
    p[i][fj]=2*p[i][fj-1]-p[i][fj-2];
    rou[i][fj]=2*rou[i][fj-1]-rou[i][fj-2];
    fadao(i,fj,visx, visy,nkk);//求向導數
//nkk[0]=-nkk[0];
uvp(vu[i][fj],vv[i][fj],rou[i][fj], p[i][fj], nkk,ubar0);//邊界處理

vu[i][fj]=ubar0[0];
vv[i][fj]=ubar0[1];
rou[i][fj]=ubar0[2];
p[i][fj]=ubar0[3];



	vu[i][0]=2*vu[i][1]-vu[i][2];
    vv[i][0]=2*vv[i][1]-vv[i][2];
    p[i][0]=2*p[i][1]-p[i][2];
    rou[i][0]=2*rou[i][1]-rou[i][2];

	vu[i][fj+1]=2*vu[i][fj]-vu[i][fj-1];
    vv[i][fj+1]=2*vv[i][fj]-vv[i][fj-1];
	p[i][fj+1]=2*p[i][fj]-p[i][fj-1];
	rou[i][fj+1]=2*rou[i][fj]-rou[i][fj-1];

}
//外場
for(j=2;j<fj;j++)
{   
	//進口處理
    
	vu[1][j]=2*vu[2][j]-vu[3][j];
    vv[1][j]=2*vv[2][j]-vv[3][j];
	p[1][j]=2*p[2][j]-p[3][j];
	rou[1][j]=2*rou[2][j]-rou[3][j];
fadao(1,j,visx, visy,nkk);//求向導數
    fadao(0,j,visx, visy,n0);//求向導數
   jin(vu[0][j],vv[0][j],rou[0][j], p[1][j],vu[1][j],vv[1][j],rou[1][j], p[1][j],n0, nkk,ubar0);
   vu[1][j]=ubar0[0];
    vv[1][j]=ubar0[1];
	p[1][j]=ubar0[3];
	rou[1][j]=ubar0[2];

 //出口處理
	vu[fi][j]=2*vu[fi-1][j]-vu[fi-2][j];
    vv[fi][j]=2*vv[fi-1][j]-vv[fi-2][j];
    p[fi][j]=2*p[fi-1][j]-p[fi-2][j];
	rou[fi][j]=2*rou[fi-1][j]-rou[fi-2][j]; 
	fadao(fi,j,visx, visy,nkk);//求向導數
    fadao(fi+1,j,visx, visy,n0);//求向導數
    chukou(vu[0][j],vv[0][j],rou[0][j], p[fi][j],vu[fi][j],vv[fi][j],rou[fi][j], p[fi][j],n0, nkk,ubar0);
    vu[fi][j]=ubar0[0];
    vv[fi][j]=ubar0[1];
	p[fi][j]=ubar0[3];
	rou[fi][j]=ubar0[2];

vu[fi+1][j]=2*vu[fi][j]-vu[fi-1][j];
    vv[fi+1][j]=2*vv[fi][j]-vv[fi-1][j];
    p[fi+1][j]=2*p[fi][j]-p[fi-1][j];
	rou[fi+1][j]=2*rou[fi][j]-rou[fi-1][j]; 
	//cin>>nl;
	
}
//四個頂點
vu[1][1]=2*vu[1][2]-vu[1][3];//[1][1]
vv[1][1]=2*vv[1][2]-vv[1][3];
p[1][1]=2*p[1][2]-p[1][3];
rou[1][1]=2*rou[1][2]-rou[1][3];

vu[1][fj]=2*vu[1][fj-1]-vu[1][fj-2];//[1][fj]
vv[1][fj]=2*vv[1][fj-1]-vv[1][fj-2];
p[1][fj]=2*p[1][fj-1]-p[1][fj-2];
rou[1][fj]=2*rou[1][fj-1]-rou[1][fj-2];

vu[fi][1]=2*vu[fi][2]-vu[fi][3];//[fi][1]
vv[fi][1]=2*vv[fi][2]-vv[fi][3];
p[fi][1]=2*p[fi][2]-p[fi][3];
rou[fi][1]=2*rou[fi][2]-rou[fi][3];

vu[fi][fj]=2*vu[fi][fj-1]-vu[fi][fj-2];//[fi][fj]
vv[fi][fj]=2*vv[fi][fj-1]-vv[fi][fj-2];
p[fi][fj]=2*p[fi][fj-1]-p[fi][fj-2];
rou[fi][fj]=2*rou[fi][fj-1]-rou[fi][fj-2];
//輸出步驟
//cin>>nl;
cout<<nl<<endl;
nl=nl-1;
}

for(i=0;i<fi;i++)
for(j=0;j<fj;j++)
{
	vu0[i][j]=vu[i+1][j+1];
    vv0[i][j]=vv[i+1][j+1];
	pp[i][j]=p[i+1][j+1];
	rour[i][j]=rou[i+1][j+1];
}

out.open("chang.txt");//將x,y坐標點寫進文件
out<<"title=\"sample mesh\"\nvariables=\"u\",\"v\",\"vu0\",\"vv0\",\"pp\",\"rour\"\nzone i="<<fi<<" j="<<fj<<" f=point\n"<<endl;
for(j=0;j<fj;j++)
for(i=0;i<fi;i++)
out<<u[i][j]<<" \t"<<v[i][j]<<" \t"<<vu0[i][j]<<" \t"<<vv0[i][j]<<" \t"<<pp[i][j]<<" \t"<<rour[i][j]<<endl;
out.close();




out.open("changzong.txt");//將x,y坐標點寫進文件
out<<"title=\"sample mesh\"\nvariables=\"u\",\"v\",\"vu0\",\"vv0\",\"pp\",\"rour\"\nzone j="<<fj<<" i="<<fi<<" f=point\n"<<endl;
for(i=0;i<fi;i++)
for(j=0;j<fj;j++)
out<<u[i][j]<<" \t"<<v[i][j]<<" \t"<<vu0[i][j]<<" \t"<<vv0[i][j]<<" \t"<<pp[i][j]<<" \t"<<rour[i][j]<<endl;
out.close();


getchar();
}





double bar(double x1,double x2,double d)
{
	double bar;
	bar=(x1-x2)/d;
	return(bar);
}
double abr(double x1,double x2,double y1,double y2)
{
	double abr;
	abr=x1*x2+y1*y2;
	return(abr);
}
double jj(double xe,double xn,double ye,double yn)
{
	double jj;
	jj=xe*yn-xn*ye;
	return(jj);
}

void kcc(int i,int j,double visx[fi+2][fj+2],double visy[fi+2][fj+2],double u,double v,double p,double rou,char k/*小寫字母*/,double (*kcc)[4])
{
	double cit,fai,en;
	double xe,xn,ye,yn;
	double jj1,kx,ky;

xe=(visx[i+1][j]-visx[i-1][j])*0.5/di;
xn=(visx[i][j+1]-visx[i][j-1])*0.5/dj;
ye=(visy[i+1][j]-visy[i-1][j])*0.5/di;
yn=(visy[i][j+1]-visy[i][j-1])*0.5/dj;
jj1=jj(xe, xn, ye, yn);
xe=xe/jj1;
xn=xn/jj1;
ye=ye/jj1;
yn=yn/jj1;
if(k=='a')
{
	kx=yn;ky=xn;
}
else{kx=ye;ky=xe;}
	en=p/(rg-1)+0.5*rou*(pow(u,2)+pow(v,2));
	en=en/rou;
fai=((rg-1)/2)*(u*u+v*v);
cit=kx*u+ky*v;
	kcc[0][0]=0;
	kcc[0][1]=kx;
	kcc[0][2]=ky;
	kcc[0][3]=0;

	kcc[1][0]=kx*fai-u*cit;
	kcc[1][1]=kx*(2-rg)*u+cit;
	kcc[1][2]=ky*u-kx*(rg-1)*v;
	kcc[1][3]=kx*(rg-1);

	kcc[2][0]=ky*fai-v*cit;
	kcc[2][1]=kx*v-ky*(rg-1)*u;
	kcc[2][2]=ky*(2-rg)*v+cit;
	kcc[2][3]=ky*(rg-1);

	kcc[3][0]=(2*fai-rg*en)*cit;
	kcc[3][1]=kx*(rg*en-fai)-(rg-1)*u*cit;
	kcc[3][2]=ky*(rg*en-fai)-(rg-1)*v*cit;
	kcc[3][3]=rg*cit;
	//cout<<endl<<endl<<"kx  ky  值"<<endl<<kx<<"  "<<ky<<"   "<<jj1<<endl;
	//cout<<xe<<"    "<<xn<<"   "<<ye<<"  "<<yn<<endl<<endl;
}

void FmB(double (*kcca)[4],double *kccg,double *kccm)
{
	int i;
for(i=0;i<4;i++)
{
kccm[i]=kcca[i][0]*kccg[0]+kcca[i][1]*kccg[1]+kcca[i][2]*kccg[2]+kcca[i][3]*kccg[3];
//cout<<endl<<"A,B"<<kccm[i]<<endl;
}

}
double temp(double p,double rou)
{
	double tem;
	tem=p/Rmg/rou;
	return(tem);
}
void FmS(double *k,double *kbar,double cita)
{
	int i;
	for(i=0;i<4;i++)
	{
		kbar[i]=k[i]*cita;
	}
}
void Fma(double *k,double *m,double *kdd,char kk)
{
	int i;
	for(i=0;i<4;i++)
	{
		if(kk=='+')
kdd[i]=k[i]+m[i];
		else
			kdd[i]=k[i]-m[i];
	}
}
void ffu(double rou,double vu,double vv,double p,double *ubar)
{
	ubar[0]=rou;
	ubar[1]=rou*vu;
	ubar[2]=rou*vv;
	ubar[3]=p/(rg-1)+0.5*rou*(pow(vu,2)+pow(vv,2));
}
void fff(double rou,double vu,double vv,double p,double *ubar)
{   
	double e;
	e=p/(rg-1)+0.5*rou*(pow(vu,2)+pow(vv,2));
    ubar[0]=rou*vu;
	ubar[1]=rou*pow(vu,2)+p;
	ubar[2]=rou*vv*vu;
	ubar[3]=(e+p)*vu;
}
void ffg(double rou,double vu,double vv,double p,double *ubar)
{
	double e;
    e=p/(rg-1)+0.5*rou*(pow(vu,2)+pow(vv,2));
    ubar[0]=rou*vv;
	ubar[1]=rou*vv*vu;
	ubar[2]=rou*pow(vv,2)+p;
	ubar[3]=(e+p)*vv;
}
void midd(double (*kcc1)[4],double (*kcc2)[4],double (*kcc)[4])
{
	int i,j;
	for(i=0;i<4;i++)
		for(j=0;j<4;j++)
		{
			kcc[i][j]=0.5*(kcc1[i][j]+kcc2[i][j]);
		}
}
void jubar(int i,int j,double visx[fi+2][fj+2],double visy[fi+2][fj+2],double rou[fi+2][fj+2],double vu[fi+2][fj+2],double vv[fi+2][fj+2],double p[fi+2][fj+2],double *ubar,double *fbar,double *gbar)
{
double xe,xn,ye,yn,jj1;
double uk[4],fk[4],gk[4];
double ubar0[4];
ffu(rou[i][j],vu[i][j], vv[i][j],p[i][j],uk);
fff(rou[i][j],vu[i][j], vv[i][j],p[i][j],fk);
ffg(rou[i][j],vu[i][j], vv[i][j],p[i][j],gk);
//求轉換系數J
xe=(visx[i+1][j]-visx[i-1][j])*0.5/di;
xn=(visx[i][j+1]-visx[i][j-1])*0.5/dj;
ye=(visy[i+1][j]-visy[i-1][j])*0.5/di;
yn=(visy[i][j+1]-visy[i][j-1])*0.5/dj;
jj1=jj(xe, xn, ye, yn);

//求轉換坐標u,f,g;
FmS(uk,ubar,jj1);

FmS(fk,fbar,yn);
FmS(gk,ubar0,-xn);
Fma(fbar,ubar0,fbar,'+');

FmS(gk,gbar,xe);
FmS(fk,ubar0,-ye);
Fma(gbar,ubar0,gbar,'+');
//cout<<endl<<"xe,xe,ye,yn"<<endl;
//cout<<jj1<<endl<<endl;
//cout<<xe<<"    "<<xn<<"   "<<ye<<"  "<<yn<<endl<<endl;
}
void fadao(int i,int j,double visx[fi+2][fj+2],double visy[fi+2][fj+2],double *n)//求向導數
 {
	 double xe,ye,xn,yn;
	 double jj0,jj1;
	 double nx,ny;
	 if(i>=1&&j>=1&&i<=fi&&j<=fj)
	 {
xe=(visx[i+1][j]-visx[i-1][j])*0.5/di;
xn=(visx[i][j+1]-visx[i][j-1])*0.5/dj;
ye=(visy[i+1][j]-visy[i-1][j])*0.5/di;
yn=(visy[i][j+1]-visy[i][j-1])*0.5/dj;}
	 else
	 {
if(i==0||j==0)
{
xe=(visx[i+1][j]-visx[i][j])/di;
xn=(visx[i][j+1]-visx[i][j])/dj;
ye=(visy[i+1][j]-visy[i][j])/di;
yn=(visy[i][j+1]-visy[i][j])/dj;
}
else
{
xe=(visx[i][j]-visx[i-1][j])/di;
xn=(visx[i][j]-visx[i][j-1])/dj;
ye=(visy[i][j]-visy[i-1][j])/di;
yn=(visy[i][j]-visy[i][j-1])/dj;
}
	 }
jj1=jj(xe, xn, ye, yn);
jj0=jj1;
nx=-ye/jj0;
ny=xe/jj0;
jj1=nx*nx+ny*ny;
jj1=sqrt(jj1);
n[0]=nx/jj1;
n[1]=ny/jj1;
 }
 void uvp(double vu,double vv,double rou,double p,double *n,double *ubar)//壁面邊界處理
{
	double a,vn;
	a=sqrt(rg*p/rou);
vn=vu*n[0]+vv*n[1];
ubar[0]=vu*pow(n[1],2)-vv*n[0]*n[1];
ubar[1]=-vu*n[0]*n[1]+vv*pow(n[0],2);
ubar[2]=pow((rg-1),2)*pow((vn-2*a/(rg-1)),2)*pow(rou,rg)/(4*rg*p);
ubar[2]=pow(ubar[2],1/(rg-1));
ubar[3]=p*pow(ubar[2],rg)/(rou,rg);
//cout<<endl<<n[0]<<"    "<<n[1]<<endl;
//cout<<endl<<vn<<"    "<<vu<<"   "<<vv<<"   "<<rou<<"   "<<p<<"    "<<a<<endl;
//cout<<endl<<ubar[0]<<"    "<<ubar[1]<<"   "<<ubar[2]<<"   "<<ubar[3]<<endl;
}
 void jin(double vu0,double vv0,double rou0,double p0,double vu,double vv,double rou,double p,double *n0,double *n,double *ubar)
 {
	 double a,vn,vt;
	 double a0,vn0,vt0;
    a=sqrt(rg*p/rou);
    a0=sqrt(rg*p0/rou0);
	vn0=vu0*n0[1]-vv0*n0[0];
	vt0=vu0*n0[0]+vv0*n0[1];
    vn=vu*n[1]-vv*n[0];
	//cout<<endl<<vn<<"    "<<vu<<"   "<<vv<<"   "<<rou<<"   "<<p<<endl;
	if(vn>a)
	{
ubar[0]=vu0;
ubar[1]=vv0;
ubar[2]=rou0;
ubar[3]=p0;
	}
	else
	{
	vn=vu*n[1]-vv*n[0];
    vn=(vn0+vn)/2+(a-a0)/(rg-1);
	vt=vt0;
	//nx,ny可能為
	if(n[0]==0)
	{
		vu=vn;
		vv=vt;
	}
	else
	{
     if(n[1]==0)
	 {
      vu=vt;
	  vv=-vn;
	 }
	 else
	 {
		 vu=(vn/n[0]+vt/n[1])/(n[1]/n[0]+n[0]/n[1]);
		 vv=vt/n[1]-vu*n[0]/n[1];
	 }
	}
ubar[0]=vu;
ubar[1]=vv;
ubar[2]=0.5*(vn-vn0)+0.5*(a0+a)/(rg-1);
ubar[2]=pow((rg-1),2)*pow(ubar[2],2)*pow(rou0,rg)/(4*rg*p0);
ubar[2]=pow(ubar[2],1/(rg-1));
ubar[3]=p0*pow(ubar[2],rg)/(rou0,rg);
	}
//cout<<endl<<vn<<"    "<<vu<<"   "<<vv<<"   "<<rou<<"   "<<p<<"    "<<a<<"  "<<a0<<endl;
//cout<<endl<<ubar[0]<<"    "<<ubar[1]<<"   "<<ubar[2]<<"   "<<ubar[3]<<n[0]<<"   "<<n[1]<<"    "<<n0[0]<<"  "<<n[1]<<endl;
 }

void chukou(double vu0,double vv0,double rou0,double p0,double vu,double vv,double rou,double p,double *n0,double *n,double *ubar)
{
	 double a,vn,vt;
	 double a0,vn0,vt0;
    a=sqrt(rg*p/rou);
    a0=sqrt(rg*p0/rou0);
	vn0=vu0*n0[1]-vv0*n0[0];
	vt0=vu0*n0[0]+vv0*n0[1];
    vn=vu*n[1]-vv*n[0];
	if(vn>a)
	{
ubar[0]=vu;
ubar[1]=vv;
ubar[2]=rou;
ubar[3]=p;
	}
	else
	{
	vn=vu*n[1]-vv*n[0];
    vn=(vn0+vn)/2+(a-a0)/(rg-1);
	vt=vt0;
	//nx,ny可能為
	if(n[0]==0)
	{
		vv=vt;
		vu=vn;
	}
	else
	{
     if(n[1]==0)
	 {
      vu=vt;
	  vv=-vn;
	 }
	 else
	 {
		 vu=(vn/n[0]+vt/n[1])/(n[1]/n[0]+n[0]/n[1]);
		 vv=vt/n[1]-vu*n[0]/n[1];
	 }
	}
ubar[0]=vu;
ubar[1]=vv;
ubar[2]=0.5*(vn-vn0)+0.5*(a0+a)/(rg-1);
ubar[2]=pow((rg-1),2)*pow(ubar[2],2)*pow(rou0,rg)/(4*rg*p0);
ubar[2]=pow(ubar[2],1/(rg-1));
ubar[3]=p0*pow(ubar[2],rg)/(rou0,rg);
	}
 }