プログラマ専用SNS ミクプラ

コード:

#include<stdio.h>
#define _USE_MATH_DEFINES
#include<math.h>

#define N 100


void Initial_z(int n,double *zx,double *zy,double p,double z_0){//zの生成
	int i;
	for(i=0;i<n;i++){
		zx[i]=p*cos(2*M_PI*(i+1)/n)+z_0;
		zy[i]=p*sin(2*M_PI*(i+1)/n);
	}
}

void Initial_xi(int n,double *xi_x,double *xi_y,double zx[N],double zy[N],double p,double q,double z_0){//ξの生成
	int i;
	for(i=0;i<n;i++){
		if(z_0<zx[i])
			xi_x[i]=zx[i]+fabs((q-p)*cos(2*M_PI*(i+1)/n));
		else
			xi_x[i]=zx[i]-fabs((q-p)*cos(2*M_PI*(i+1)/n));
		xi_y[i]=q*zy[i];
	}
}

void Initial_GN(int n,double GN[N][N],double zx[N],double zy[N],double xi_x[N],double xi_y[N]){//係数行列の生成
	int i,j;
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			GN[i][j]=sqrt((zx[i]-xi_x[j])*(zx[i]-xi_x[j])+(zy[i]-xi_y[j])*(zy[i]-xi_y[j]));
			GN[i][j]=log(GN[i][j]);
		}
	}
}

void Initial_logz(int n,double zx[N],double zy[N],double *logz){//連立1次方程式の右辺を生成
	int i;
	for(i=0;i<n;i++)
		logz[i]=log(sqrt(zx[i]*zx[i]+zy[i]*zy[i]));
}

void Gauss_Method(int n,double G[N][N],double f[N],double *Q){//連立一次方程式を解く関数
	int i,j,k;
	double s;	
	/*前進消去*/
	for(i=0;i<n-1;i++){
		for(j=i+1;j<n;j++){
			s=G[j][i]/G[i][i];
			for(k=i+1;k<n;k++)
				G[j][k]=G[j][k]-G[i][k]*s;
			f[j]=f[j]-f[i]*s;
		}
	}
	/*後進消去*/
	for(i=n-1;i>=0;i--){
		s=f[i];
		for(j=i+1;j<n;j++)
			s=s-G[i][j]*f[j];
		f[i]=s/G[i][i];
		Q[i]=f[i];
	}
}
void Calc_G(int n,double Q[N],double zx[N],double zy[N],double xi_x[N],double xi_y[N],double *G){//Gの計算
	int i,j;
	double sum,num;
	for(i=0;i<n;i++){
		sum=0;
		for(j=0;j<n;j++){
			num=sqrt((zx[i]-xi_x[j])*(zx[i]-xi_x[j])+(zy[i]-xi_y[j])*(zy[i]-xi_y[j]));
			sum+=Q[j]*log(num);
		}
		G[i]=-1*sum;
	}
}
void Calc_H(int n,double Q[N],double zx[N],double zy[N],double xi_x[N],double xi_y[N],double *H){//Hの計算
	int i,j;
	double sum,x,y;
	for(i=0;i<n;i++){
		sum=0;
		for(j=0;j<n;j++){
			x=1-(zx[i]*xi_x[j]+zy[i]*xi_y[j])/((xi_x[j]*xi_x[j])+(xi_y[j]*xi_y[j]));
			y=(zx[i]*xi_y[j]-zy[i]*xi_x[j])/((xi_x[j]*xi_x[j])+(xi_y[j]*xi_y[j]));
			sum+=Q[j]*atan2(y,x);
		}
		H[i]=-1*sum;
	}
}

void Calc_UV(int n,double G[N],double H[N],double zx[N],double zy[N],double *U,double *V){//U,Vの計算
	int i;
	for(i=0;i<n;i++){
		U[i]=exp(log(sqrt(zx[i]*zx[i]+zy[i]*zy[i]))+G[i])*cos(atan2(zx[i],zy[i])+H[i]);
		V[i]=exp(log(sqrt(zx[i]*zx[i]+zy[i]*zy[i]))+G[i])*sin(atan2(zx[i],zy[i])+H[i]);
	}
}

int main(){
	int n=16,i;
	double zx[N],zy[N],xi_x[N],xi_y[N],GN[N][N],logz[N],Q[N],u[N],v[N],W[N];//計算上必要な配列
	double G[N],H[N],U[N],V[N];//最終的に求めたいものを格納する配列
	double p=1.0;//拘束点の半径
	double q=1.2;//電荷点の拡大率
	double z_0=-0.25;//正規化条件

	Initial_z(n,zx,zy,p,z_0);//拘束点zの生成
	Initial_xi(n,xi_x,xi_y,zx,zy,p,q,z_0);//電荷点ξの生成
	Initial_GN(n,GN,zx,zy,xi_x,xi_y);//連立1次方程式の係数行列を生成
	Initial_logz(n,zx,zy,logz);//連立1次方程式の右辺を生成
	Gauss_Method(n,GN,logz,Q);//電荷Qを求める。
	Calc_G(n,Q,zx,zy,xi_x,xi_y,G);//Gの計算
	Calc_H(n,Q,zx,zy,xi_x,xi_y,H);//Hの計算
	Calc_UV(n,G,H,zx,zy,U,V);//U,Vの生成
	return 0;
}