#define WX (uint)128
#define WY (uint)128
#define ipn (uint)4
#define PRESS_HYOJI_XY (uint)20

typedef struct{
	float x;
	float y;
	float z;
} Vec3;

typedef struct{
	float x;
	float y;
	float z;
	float w;
} Vec4;


void mulqarray(Vec4 q1, Vec4 q2,Vec4 *q3)
{
	q3->x=q1.x*q2.x-q1.y*q2.y-q1.z*q2.z-q1.w*q2.w;
	q3->y=q1.x*q2.y+q2.x*q1.y +q1.z*q2.w-q1.w*q2.z;
	q3->z=q1.x*q2.z+q2.x*q1.z +q1.w*q2.y-q1.y*q2.w;
	q3->w=q1.x*q2.w+q2.x*q1.w +q1.y*q2.z-q1.z*q2.y;
}


//流体にUFO壁を適応するためにまずUFO回転
__kernel void UfoRote(__global Vec3 *dCLid,__global uchar *ufom,float qsx,float qsy,float rad2)
{
	float q20=cos(rad2);
	float q22=sin(rad2);
	int psx=((int)qsx)-6;
	int psy=((int)qsy)-6;
	int icx=get_global_id(0);

	Vec3 vert=dCLid[icx];
	Vec4 P;
	P.x=0.0;
	P.y=vert.x;
	P.z=vert.y;
	P.w=vert.z;
	Vec4 R;
	R.x=q20;
	R.y=-0.0;
	R.z=-q22;
	R.w=-0.0;
	Vec4 Q;
	Q.x=q20;
	Q.y=0.0;
	Q.z=q22;
	Q.w=0.0;
	Vec4 tmp;
	mulqarray(R,P,&tmp);
	mulqarray(tmp,Q,&P);
	vert.x=P.y+qsx;
	vert.y=P.z;
	vert.z=P.w+qsy;
//	vnuid[icx]=vert;
//	Vec3 vert=vnuid[icx];
	float Px1=vert.x-psx;
	float Pz1=vert.z-psy;
	if (Px1<  0.0){Px1=0.0;}
	if (Px1>=12.0){Px1=11.99999;}
	if (Pz1<  0.0){Pz1=0.0;}
	if (Pz1>=12.0){Pz1=11.99999;}
	uint i1=(int)Px1;
	uint i2=(int)Pz1;
	uint tmpp=mad24(i2,(uint)12,i1);
	if (tmpp<12*12){
		ufom[tmpp]=0;
	}
}


//UFOの個体のところを壁にする
//もともと壁=0、流体=255のところを、UFOの壁適応で壁=1にする。あとで解除するので絶対0にはしてはいけない
//UFOの速度を壁に適応。ここで回転速度も考慮する
__kernel void UfoMap(__global uchar *ufom,__global uchar *kabeP,__global uchar *kabeX,__global uchar *kabeY,int psx,int psy,__global float *YUN,__global float *YVN,__global float *GXU,__global float *GYU,__global float *GXV,__global float *GYV,float spdx,float spdy,float ufoQuoSpd)
{
	int gid=get_global_id(0);
	int icy=gid/12;
	int icx=gid-icy*12;
	uchar kabeda=ufom[gid];
	float radSpdX=(float)icx-6.0;
	float radSpdY=(float)icy-6.0;
	float len=sqrt(radSpdX*radSpdX+radSpdY*radSpdY);

	float radTmp=atan2(radSpdY,radSpdX)+3.14159265359/2.0;

	float yajix=cos(radTmp)*len*ufoQuoSpd;
	float yajiy=sin(radTmp)*len*ufoQuoSpd;

	icx+=psx;
	icy+=psy;

	uint gid2=icx+icy*WX;
	if (gid2<(WX*WY-WX-1)){
		if (kabeda==0){
			if (kabeP[gid2   ]==255){kabeP[gid2   ]=1;}
			if (kabeX[gid2   ]==255){kabeX[gid2   ]=1;YUN[gid2   ]=spdx+yajix;GXU[gid2   ]=0.0f;GYU[gid2   ]=0.0f;}
			if (kabeX[gid2+1 ]==255){kabeX[gid2+1 ]=1;YUN[gid2+1 ]=spdx+yajix;GXU[gid2+1 ]=0.0f;GYU[gid2+1 ]=0.0f;}
			if (kabeY[gid2   ]==255){kabeY[gid2   ]=1;YVN[gid2   ]=spdy+yajiy;GXV[gid2   ]=0.0f;GYV[gid2   ]=0.0f;}
			if (kabeY[gid2+WX]==255){kabeY[gid2+WX]=1;YVN[gid2+WX]=spdy+yajiy;GXV[gid2+WX]=0.0f;GYV[gid2+WX]=0.0f;}
		}
	}
}

//Ufoマッピングで壁になったところを解除
//ufomはUFOマスク。0だと壁、255だと流体
__kernel void UfoUnMap(__global uchar *kabeP,__global uchar *kabeX,__global uchar *kabeY,int psx,int psy)
{
	int gid=get_global_id(0);
	int icy=gid/12;
	int icx=gid-icy*12;
	icx+=psx;
	icy+=psy;
	uint gid2=icx+icy*WX;
	if (gid2<(WX*WY-WX-1)){
		if (kabeP[gid2   ]==1){kabeP[gid2]=255;}
		if (kabeX[gid2   ]==1){kabeX[gid2]=255;}
		if (kabeX[gid2+1 ]==1){kabeX[gid2+1]=255;}
		if (kabeY[gid2   ]==1){kabeY[gid2]=255;}
		if (kabeY[gid2+WX]==1){kabeY[gid2+WX]=255;}
	}
}

//噴射口1ドットの湧出設定
__kernel void UfoHun(int psxy,float rada,__global uchar *kabeP,__global uchar *kabeX,__global uchar *kabeY,__global float *YUN,__global float *YVN,__global float *GXU,__global float *GYU,__global float *GXV,__global float *GYV)
{
	uint posxy=psxy;
	float hunspdx=0.9*sin(rada);
//	float sgnx=+1.0;
//	if (hunspdx<0.0){sgnx=+1.0;}
	float hunspdy=0.9*cos(rada);
//	float sgny=+2.0;
//	if (hunspdy<0.0){sgny=+2.0;}

	if (posxy<(WX*WY-WX-1)){
		if (kabeP[posxy   ]==255){kabeP[posxy   ]=1;}
		if (kabeX[posxy   ]==255){kabeX[posxy   ]=1;YUN[posxy   ]=hunspdx;}//GXU[posxy    ]=-1.0;GYU[posxy    ]= 0.0;
		if (kabeX[posxy+1 ]==255){kabeX[posxy+1 ]=1;YUN[posxy+1 ]=hunspdx;}//GXU[posxy+1  ]=-1.0;GYU[posxy+1  ]= 0.0;
		if (kabeY[posxy   ]==255){kabeY[posxy   ]=1;YVN[posxy   ]=hunspdy;}//GXV[posxy    ]= 0.0;GYV[posxy    ]=-1.0;
		if (kabeY[posxy+WX]==255){kabeY[posxy+WX]=1;YVN[posxy+WX]=hunspdy;}//GXV[posxy+128]= 0.0;GYV[posxy+128]=-1.0;
	}
}


//噴射口1ドットの湧出設定解除
__kernel void UfoUnHun(int psxy,__global uchar *kabeP,__global uchar *kabeX,__global uchar *kabeY)
{
	uint posxy=psxy;
	if (posxy<(WX*WY-WX-1)){
		if (kabeP[posxy   ]==1){kabeP[posxy   ]=255;}
		if (kabeX[posxy   ]==1){kabeX[posxy   ]=255;}
		if (kabeX[posxy+1 ]==1){kabeX[posxy+1 ]=255;}
		if (kabeY[posxy   ]==1){kabeY[posxy   ]=255;}
		if (kabeY[posxy+WX]==1){kabeY[posxy+WX]=255;}
	}
}


//UFOが流れ圧力から外力を受ける計算
__kernel void UfoGairyoku(int psxy,__global uchar *kabeP,__global float *YPN,__global float *GAIRYOKU)
{
	uint gid=get_global_id(0);
	uint gid2=gid*6;
	uint icy=gid/12;
	uint icx=gid-icy*12;
	uint tmp0=psxy+icx+icy*WX;


	GAIRYOKU[gid2  ]=0.0f;
	GAIRYOKU[gid2+1]=0.0f;
	GAIRYOKU[gid2+2]=0.0f;
	GAIRYOKU[gid2+3]=0.0f;
	GAIRYOKU[gid2+4]=0.0f;
	GAIRYOKU[gid2+5]=0.0f;
	float yokokabefg=0.0;//左が流体なら+256、右が流体なら+1
	float tatekabefg=0.0;//上が流体なら+256、下が流体なら+1
	float yokop0=0.0f;
	float tatep0=0.0f;
	float yokop1=0.0f;
	float tatep1=0.0f;
	if (tmp0>(WX+1)){if (tmp0<(WX*WY-WX-1)){
		if (kabeP[tmp0]==1){
			if (kabeP[tmp0-1]==255){yokop0+=YPN[tmp0-1];yokokabefg+=256.0;}
			if (kabeP[tmp0+1]==255){yokop1-=YPN[tmp0+1];yokokabefg+=1.0;}
			if (kabeP[tmp0-WX]==255){tatep0+=YPN[tmp0-WX];tatekabefg+=256.0;}
			if (kabeP[tmp0+WX]==255){tatep1-=YPN[tmp0+WX];tatekabefg+=1.0;}
			GAIRYOKU[gid2  ]=yokop0;
			GAIRYOKU[gid2+1]=tatep0;
			GAIRYOKU[gid2+2]=yokop1;
			GAIRYOKU[gid2+3]=tatep1;
			GAIRYOKU[gid2+4]=yokokabefg;
			GAIRYOKU[gid2+5]=tatekabefg;
		}
	}}
}

//UFOが流れから外力を受ける。圧力を合計
__kernel void UfoGairyoku2(__global float *GAIRYOKU)
{
	uint gid=get_global_id(0);
	float tmpx=0.0f;

	for(int i=gid;i<144*6;i+=6){
		tmpx+=GAIRYOKU[i];
	}

	GAIRYOKU[gid]=tmpx;
}



//爆発時にDIVに加算する関数
__kernel void divexp(__global uchar *kabeP,__global uchar *kabeX,__global uchar *kabeY,int psx,int psy,__global float *DIV,float exploF) {
	int gid=get_global_id(0);
	int icy=gid/12;
	int icx=gid-icy*12;
	icx+=psx;
	icy+=psy;
	uint gid2=icx+icy*WX;
	if (gid2<(WX*WY-WX-1)){
		if (kabeP[gid2   ]==1){kabeP[gid2]=255;DIV[gid2]-=exploF;}
		if (kabeX[gid2   ]==1){kabeX[gid2]=255;}
		if (kabeX[gid2+1 ]==1){kabeX[gid2+1]=255;}
		if (kabeY[gid2   ]==1){kabeY[gid2]=255;}
		if (kabeY[gid2+WX]==1){kabeY[gid2+128]=255;}
	}
}



//右の画面の圧力のやつ
//PRESS_HYOJI_XY=20という前提
__kernel void presshyoji(__global uchar *kabeP,__global float *YPN,__global uint *clmem_presshyouji_color,int psx,int psy,__global float *clmem_presshyouji_sum){
	int ipsx=psx-PRESS_HYOJI_XY/2;//圧力取得座標hの左上
	int ipsy=psy-PRESS_HYOJI_XY/2;//圧力取得座標hの左上

	int gid=get_global_id(0);
	int gidy=gid/PRESS_HYOJI_XY;
	int gidx=gid-gidy*PRESS_HYOJI_XY;
	ipsx+=gidx;
	ipsy+=gidy;
	int indx=ipsy*WX+ipsx;
	float yptmp=0.0;
	float yptmr=0.0;
	float yptmg=0.0;
	float yptmb=0.0;

	if ((ipsx>=0)&(ipsy>=0)&(ipsx<WX)&(ipsy<WY)){
		if ((uint)kabeP[indx]<4){//壁なら白にする
			clmem_presshyouji_color[gid*4  ]=4278848010;
			clmem_presshyouji_color[gid*4+1]=4278848010;
			clmem_presshyouji_color[gid*4+2]=4278848010;
			clmem_presshyouji_color[gid*4+3]=4278848010;
		}else{//流体ならちゃんと計算して色代入
			yptmp=YPN[indx]-clmem_presshyouji_sum[0];
			if (yptmp<0.0){
				yptmb=log(0.12*(-1000.0*yptmp)+1.0)*45.0;
				yptmg=-yptmp*321.9;
			}else{
				yptmr=log(0.12*(1000.0*yptmp)+1.0)*45.0;
				yptmg=yptmp*321.9;
			}//陰圧なら青、陽圧なら赤
			if (yptmr<0.0){yptmr=0.0;}
			if (yptmg<0.0){yptmg=0.0;}
			if (yptmb<0.0){yptmb=0.0;}
			if (yptmr>255.0){yptmr=255.0;}
			if (yptmg>255.0){yptmg=255.0;}
			if (yptmb>255.0){yptmb=255.0;}
			clmem_presshyouji_color[gid*4  ]=((uint)yptmr)+((uint)yptmg)*256+((uint)yptmb)*65536+(uint)4278190080;
			clmem_presshyouji_color[gid*4+1]=((uint)yptmr)+((uint)yptmg)*256+((uint)yptmb)*65536+(uint)4278190080;
			clmem_presshyouji_color[gid*4+2]=((uint)yptmr)+((uint)yptmg)*256+((uint)yptmb)*65536+(uint)4278190080;
			clmem_presshyouji_color[gid*4+3]=((uint)yptmr)+((uint)yptmg)*256+((uint)yptmb)*65536+(uint)4278190080;
		}
	}else{//範囲外なら0代入
		clmem_presshyouji_color[gid*4  ]=255*65536*256;
		clmem_presshyouji_color[gid*4+1]=255*65536*256;
		clmem_presshyouji_color[gid*4+2]=255*65536*256;
		clmem_presshyouji_color[gid*4+3]=255*65536*256;
	}
}


//圧力の合計の値を計算。
__kernel void presshyoji2(__global float *YPN,__global float *clmem_presshyouji_sum,__local float *Psum,int psx,int psy,__global uchar *kabeP) {
	int ipsxd=psx-PRESS_HYOJI_XY/2;//圧力取得座標hの左上
	int ipsyd=psy-PRESS_HYOJI_XY/2;//圧力取得座標hの左上
	int gid=get_global_id(0);
	int gidy=0;
	int gidx=0;
	int indx=0;
	int ipsx=0;
	int ipsy=0;
	float flg=0.0;
	float tmpPsum=0.0;
	//測れる部分の圧力を自分のレジスタに
	for(uint i=gid;i<PRESS_HYOJI_XY*PRESS_HYOJI_XY;i+=64){
		gidy=i/PRESS_HYOJI_XY;
		gidx=i-gidy*PRESS_HYOJI_XY;
		ipsx=ipsxd+gidx;
		ipsy=ipsyd+gidy;
		if ((ipsx>=0)&(ipsy>=0)&(ipsx<WX)&(ipsy<WY)){
			indx=ipsy*WX+ipsx;
			if (kabeP[indx]>=4){
				tmpPsum+=YPN[indx];
				flg+=1.0;
			}
		}
	}
	//共有メモリにコピー
	Psum[gid]	=tmpPsum;
	Psum[gid+64]=flg;
	barrier(CLK_LOCAL_MEM_FENCE);
	//だんだんと合計していく
	for(uint i=32;i>=1;i/=2){
		//iは32 16 8 4 2 1となっていく
		if (gid<i){
			Psum[gid]	+=Psum[gid+i];
			Psum[gid+64]+=Psum[gid+i+64];
		}
		barrier(CLK_LOCAL_MEM_FENCE);
	}

	if (gid==0){clmem_presshyouji_sum[0]=Psum[0]/Psum[64];}

}


