複素数の計算についての質問

Question

前提

複素数の計算についての質問

実現したいこと

下のコードにあるよう、行列の中の要素を実部、虚部に分けて別々の行列の中に入れたい。

発生している問題・エラーメッセージ

行列のA、Bの奇数行目の値がすべて0になる
ほかの複素数で確かめても0になることはなく、もともとの複素数も0でない

該当のソースコード

//ヤコビ法による固有値と固有ベクトル計算//

#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>

#define N  4			//次数設定
#define L  8
#define EPS 0.0001  	//収束範囲

int n = 100000;
int M = 2;                           //電子数＝2M
int c = 1;
double complex t_1, t_2, t_3, t_4, t_5, t_6,T_1,T_2;		//方向
double b_pi;
double a[L][L];
double  D[L];



int main(int argc, char** argv){

	b_pi = M_PI;
	double complex a_cube = 0.246 * pow(10, -9); //格子定数[nm]
	double complex k_x = -b_pi / a_cube; //↓のk_xの値を変えるのも忘れずに！
	double complex k_y = b_pi / M * c;

	FILE* fp;

	if ((fp = fopen("zig N=2 q=1.txt", "w")) == NULL)
	{
		printf("Cannot open the file\n");
		exit(1);
	}

	int v;
	for (v = 0;v < n; v++) {
		k_x = k_x + b_pi * 2 / a_cube / n;
			t_4 = cos(-1 * a_cube / pow(3, 0.5) * k_y) + I * sin(-1 * a_cube / pow(3, 0.5) * k_y);	//cexp(I * (-1* a_cube / pow(3, 0.5) * k_y));
			t_5 = cexp(I * (a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y)); //cos(a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y) + I * sin((a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y));	//cexp(I * (a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y));  //expの中の符号を計算したあらかじめ
			t_6 = cexp(I * (a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y)); //cos(a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y) + I * sin(a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y);	//cexp(I * (a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y));  //
			t_1 = cos(a_cube / pow(3, 0.5) * k_y) + I * sin(a_cube / pow(3, 0.5) * k_y);	//cexp(I * a_cube / pow(3, 0.5) * k_y);
			t_2 = cos(-a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y) + I * sin(-a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y);	//cexp(I * (-a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y));
			t_3 = cos(-a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y) + I * sin(-a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y);	//cexp(I * (-a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y));
			T_1= t_5 + t_6;
			T_2= t_2 + t_3;

			double complex C[N][N] = {
								 {0,0,T_1,0},
								 {0,0,t_4,T_1},
								 {T_2,t_1, 0, 0},
								 {0,T_2, 0, 0} };		//係数行列
		
			double A[N][N], B_1[N][N], B_2[N][N];
			int k, h;
			for (k = 0;k < N;k++) {
				for (h = 0;h < N;h++) {
					A[k][h] = creal(C[k][h]);
					B_1[k][h] = -1 * cimag(C[k][h]);
					B_2[k][h] = cimag(C[k][h]);
				}
			}
			/*for (k = 0;k < N;k++) {
				for (h = 0;h < N;h++) {
					printf("%7.4lf\n", A[k][h]);
				}
			}*/
			printf("%2f  %lf %lf %lf %lf\n", k_x, B_1[1][0], B_1[1][1], B_1[1][2], B_1[1][3]);
			printf("%2f  %lf %lf %lf %lf\n", k_x, A[1][0], A[1][1], A[1][2], A[1][3]);

			/*= {			{A[0][0],A[0][1],A[0][2] ,B_1[0][0],B_1[0][1],B_1[0][2]},
							{A[1][0],A[1][1],A[1][2] ,B_1[1][0],B_1[1][1],B_1[1][2]},
							{A[2][0],A[2][1],A[2][2] ,B_1[2][0],B_1[2][1],B_1[2][2]},
							{B_2[0][0],B_2[0][1],B_2[0][2] ,A[0][0],A[0][1],A[0][2]},
							{B_2[1][0],B_2[1][1],B_2[1][2] ,A[1][0],A[1][1],A[1][2]},
							{B_2[2][0],B_2[2][1],B_2[2][2] ,A[2][0],A[2][1],A[2][2]}};*/
			for (k = 0;k < L;k++) {
				for (h = 0;h < L;h++) {
					if (h < N) {
						if (k < N)
							a[h][k] = A[h % N][k % N];
						else
							a[h][k] = B_1[h % N][k % N];
					}
					else {
						if (k < N)
							a[h][k] = B_2[h % N][k % N];
						else
							a[h][k] = A[h % N][k % N];
					}
				}
			}

			printf("%f%+fi\n", creal(T_2), cimag(T_2));
			//printf("%2f  %lf %lf %lf %lf %lf %lf %lf %lf\n", k_x, a[1][0], a[1][1], a[1][2], a[1][3], a[1][4], a[1][5], a[1][6], a[1][7]);
			/*/for (k = 0;k < L;k++) {
				for (h = 0;h < L;h++)
					printf("%7.4lf\n", a[k][h]);
				}*/

			//double u[L][L];					//単位行列
			double alpha, beta, gamma;
			double s, c, w;
			double wa, wb, wc;
			double max;
			int i, j, p, q, x, y;

			/*for (y = 0;y < N;y++)
			for (x = 0;x < N;x++)
				u[y][x] = 0.0;

		for (i = 0;i < N;i++)
			u[i][i] = 1.0;*/

		while (1) {
			//最大要素の行と列を検索
			max = 0.0;
			for (i = 0;i < L - 1;i++) {
				for (j = i + 1;j < L;j++)
					if (fabs(a[i][j]) > max) {
						p = i;
						q = j;
						max = fabs(a[i][j]);
					}
			}
			//収束したら解答打ち出し
			if (max < EPS) break;

			//sin, cos　計算
			wa = a[p][p];
			wb = a[p][q];
			wc = a[q][q];
			alpha = -wb;
			beta = 0.5 * (wa - wc);
			gamma = fabs(beta) / sqrt(alpha * alpha + beta * beta);
			s = sqrt(0.5 * (1.0 - gamma));
			if (alpha * beta < 0) s = -s;
			c = sqrt(1.0 - s * s);
			//直行変換
			for (j = 0;j < L;j++) {
				w = a[p][j] * c - a[q][j] * s;
				a[q][j] = a[p][j] * s + a[q][j] * c;
				a[p][j] = w;
			}

			for (j = 0;j < L;j++) {
				a[j][p] = a[p][j];
				a[j][q] = a[q][j];
			}
			w = 2.0 * wb * s * c;
			a[p][p] = wa * c * c + wc * s * s - w;
			a[q][q] = wa * s * s + wc * c * c + w;
			a[p][q] = 0;
			a[q][p] = 0;

			D[0] = a[0][0];
			D[1] = a[1][1];
			D[2] = a[2][2];
			D[3] = a[3][3];
			D[4] = a[4][4];
			D[5] = a[5][5];
			D[6] = a[6][6];
			D[7] = a[7][7];

			//漸化式計算
			/*for (i = 0;i < N;i++) {
				w = u[i][p] * c - u[i][q] * s;
				u[i][q] = u[i][p] * s + u[i][q] * c;
				u[i][p] = w;
			}*/

			/*for (k = 0;k < L;k++) {
				for (h = 0;h < L; h++) {
					printf("%7.4lf\n", a[k]);
				}*/
		}
		
		fprintf(fp, "%2f  %lf %lf %lf %lf %lf %lf %lf %lf\n", k_x, D[0], D[1], D[2], D[3], D[4], D[5], D[6], D[7]);
	}

	/*int i;
	printf("固有値\n");
	for (i = 0;i < L;i++) {
		printf("%7.4lf\n", a[i][i]);			//対角成分　こいつが固有値
		}*/
		
	fclose(fp);
	return 0;
}

試したこと

行列の値をいちいち確かめたり、複素数自体の値も調べた
前回のミスもこの点が問題であると思う。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Answer 1

勘違いしていることにすぐ気づいた