#include
#include
#define N 5 /* 行列のサイズ */
#define M 10 /*点の個数*/
#define eb 3.0
#define ea 1.0
#define a 7.0
#define w 2.0
#define PI2 6.28318530717959
main(void)
{
int i, j;
float A[N][N], B[N][N];
double P[N][N] = {{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,1,0},{0,0,0,0,1}};
double P2[N][N] = {{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,1,0},{0,0,0,0,1}};
double x,k,c;
int m,q,l,n,b;
for(c = 0.0; c < M; c++){
k=c/M*PI2/(2*a);
printf("k = %f,",k);
/*行列の初期値を代入する*/
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
if (i == j)
A[i][i] = (1/eb+w/a*(1/eb-1/ea))*(k+PI2/a*(j-2))*(k+PI2/a*(i-2));
else
A[i][j] = (1/ea-1/eb)*w/a*sin*1/*2*(k+PI2/a*(j-2))*(k+PI2/a*(i-2));
}
}
for( b=0; b<10; b++){
/*非対角成分で最大の成分を見つける*/
/*01-02-03-04-12-13-14-23-24-34の順に成分を比較する*/
q = 0, l = 1,m = 0,n = 2,x = 0.0;
while (m < 4){
if (A[q][l] < A[m][n]) {
q = m;
l = n;
}
if (n < 4)
n++;
else {
m++;
n = m+n-3;
}
}
/*ユニタリ行列PとPの転置行列をつくる*/
x = 1/2*atan(2*A[q][l]/(A[q][q]-A[l][l]));
P[q][q] = P[l][l] = P2[q][q] = P2[l][l] = cos(x);
P[q][l] = P2[l][q] = -sin(x);
P[l][q] = P2[q][l] = sin(x);
/*行列の計算を行う*/
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++){
B[i][j] = 0;
for (m = 0; m < N; m++){
B[i][j] += A[i][m]*P[m][j];
}
}
}
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++){
A[i][j] = 0;
for (m = 0; m< N; m++){
A[i][j] += P2[i][m]*B[m][j];
}
}
}
A[q][l] = A[l][q] = 0;
}
for(i=0;i
