前提
状態密度を計算したい。
実現したいこと
そのために、eとk_xの変化ごとにsum_densityを計算して、最後にeごとにsum_dosの値をファイルに書きたい。
発生している問題・エラーメッセージ
現状、eが100回増加を繰り返した後、k_xは100の二乗回増加されています。これをeが一回増加ごとにk_xが100回増加してそれごとにsum_densityを計算して。。。みたいにしたいです。
C
1#include <stdio.h> 2#include <math.h> 3#include <complex.h> 4#include <stdlib.h> 5 6#define M 5 //単位法の数 7#define N 10 //次数設定 8#define L 20 9 10#define EPS 0.000001 //収束範囲 11 12int n = 100; //100000 //電子数=2N 13int c = 1; 14double complex t_1, t_2, t_3, t_4, t_5, t_6, t_7, t_8, T_1, T_2, T_3, T_4; //方向 15double b_pi; 16double a[L][L]; 17double D[N][N]; 18double complex C[N][N]; 19double A[L], B_1[N][N], B_2[N][N]; 20double E[L]; 21 22 23 24 25int main(int argc, char** argv) { 26 27 b_pi = M_PI; 28 double a_cube = 1;//= 0.246 * pow(10, -9); //格子定数[nm] 29 double k_x = -1.5*b_pi / a_cube; //-2*b_pi/a_cube/pow(3,0.5); //↓のk_xの値を変えるのも忘れずに! K+,K-で値が変わるよ! 30 double k_y = 0; 31 double e=-3.5; 32 double sum_dos; 33 double sum_density[n]; 34 35 for(int i=0;i<n;i++){ 36 sum_dos=0.0; 37 sum_density[i]=0.0; 38 } 39 40 FILE* fp; 41 42 if ((fp = fopen("aaa N=5.txt", "w")) == NULL) 43 { 44 printf("Cannot open the file\n"); 45 exit(1); 46 } 47 48 for(int h=0;h<n;h++){ 49 e=e+7.0/n; 50 51 52 for (int r = 0;r < n; r++) { 53 k_x = k_x + b_pi * 3/ a_cube / n; 54 t_4 = cexp(-1 * I * (-1 * a_cube / pow(3, 0.5) * k_y)); 55 t_5 = cexp(-1 * I * (-a_cube / 2 * k_x - a_cube / 2 / pow(3, 0.5) * 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の中の符号を計算したあらかじめ 56 t_6 = cexp(-1 * I * (a_cube / 2 * k_x - a_cube / 2 / pow(3, 0.5) * 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)); // 57 t_1 = cexp(I * a_cube / pow(3, 0.5) * k_y); 58 t_2 = cexp(I * (-a_cube / 2 * k_x - a_cube / 2 / pow(3, 0.5) * k_y)); 59 t_3 = cexp(I * (a_cube / 2 * k_x - a_cube / 2 / pow(3, 0.5) * k_y)); 60 t_7 = -1 * t_4; 61 t_8 = -1 * t_1; 62 T_1 = t_5 + t_6; 63 T_2 = t_2 + t_3; 64 T_3 = -1 * T_1; 65 T_4 = -1 * T_2; 66 67 for (int k = 0;k < N;k++) { 68 for (int h = 0;h < N;h++) { 69 for (int i = 0;i < M - 2;i++) { 70 if (h < M) { 71 if (k < M) { 72 C[h][k] = 0; 73 } 74 else { 75 C[0][M] = T_1; 76 C[M - 1][N - 1] = T_1; 77 C[M - 1][N - 2] = t_4; 78 C[i + 1][M + i] = t_4; 79 C[i + 1][M + i + 1] = T_1; 80 } 81 } 82 else { 83 if (k < M) { 84 C[M][0] = T_2; 85 C[M][1] = t_1; 86 C[N - 1][M - 1] = T_2; 87 C[M + i + 1][i + 1] = T_2; 88 C[M + i + 1][i + 2] = t_1; 89 90 } 91 else { 92 C[h][k] = 0; 93 } 94 } 95 } 96 } 97 } 98 99 int k, h; 100 for (k = 0;k < N;k++) { 101 for (h = 0;h < N;h++) { 102 D[h][k] = creal(C[h][k]); 103 B_1[h][k] = -1 * cimag(C[h][k]); 104 B_2[h][k] = cimag(C[h][k]); 105 } 106 } 107 108 109 for (k = 0;k < L;k++) { 110 for (h = 0;h < L;h++) { 111 if (h < N) { 112 if (k < N) { 113 a[h][k] = D[h % N][k % N]; 114 } 115 else { 116 a[h][k] = B_1[h % N][k % N]; 117 } 118 } 119 else { 120 if (k < N) { 121 a[h][k] = B_2[h % N][k % N]; 122 } 123 else { 124 a[h][k] = D[h % N][k % N]; 125 } 126 } 127 } 128 } 129 130 while (1) { 131 //double u[L][L]; //単位行列 132 double alpha, beta, gamma; 133 double s, c, w; 134 double wa, wb, wc; 135 double max; 136 int i, j, p, q, x, y; 137 138 //最大要素の行と列を検索 139 max = 0.0; 140 for (i = 0;i < L - 1;i++) { 141 for (j = i + 1;j < L;j++) 142 if (fabs(a[i][j]) > max) { 143 p = i; 144 q = j; 145 max = fabs(a[i][j]); 146 } 147 } 148 149 //収束したら解答打ち出し 150 if (max < EPS) break; 151 152 //sin, cos 計算 153 wa = a[p][p]; 154 wb = a[p][q]; 155 wc = a[q][q]; 156 alpha = -wb; 157 beta = 0.5 * (wa - wc); 158 gamma = fabs(beta) / sqrt(alpha * alpha + beta * beta); 159 s = sqrt(0.5 * (1.0 - gamma)); 160 if (alpha * beta < 0) s = -s; 161 c = sqrt(1.0 - s * s); 162 //直行変換 163 for (j = 0;j < L;j++) { 164 w = a[p][j] * c - a[q][j] * s; 165 a[q][j] = a[p][j] * s + a[q][j] * c; 166 a[p][j] = w; 167 } 168 169 for (j = 0;j < L;j++) { 170 a[j][p] = a[p][j]; 171 a[j][q] = a[q][j]; 172 } 173 174 175 176 w = 2.0 * wb * s * c; 177 a[p][p] = wa * c * c + wc * s * s - w; 178 a[q][q] = wa * s * s + wc * c * c + w; 179 a[p][q] = 0; 180 a[q][p] = 0; 181 182 for (int i = 0;i < L;i++) 183 { 184 A[i] = a[i][i]; 185 } 186 187 } 188 189 double tmp; 190 for (int i=0; i<L; i++) { 191 for (int j=i+1; j<L ;j++) { 192 if (A[i] > A[j]) { 193 tmp = A[i]; 194 A[i] = A[j]; 195 A[j] = tmp; 196 } 197 } 198 } 199 double ave, squ, v, s; 200 int data_num; 201 double sum = 0.0, sumv = 0.0; 202 203 ave = 0.0; 204 squ = 0.0; 205 v = 0.0; 206 s = 0.0; 207 data_num = L; 208 sum = 0.0; 209 sumv = 0.0; 210 for(int i = 0; i < data_num; i++){ 211 sum += A[i]; // データの総和 212 sumv += A[i] * A[i]; // 二乗和 213 } 214 215 // 平均・平方和・分散・標準偏差を求める 216 ave = sum / data_num; // 平均 217 squ = sumv - (sum * sum / data_num); // 平方和 218 v = squ / (data_num - 1); // 分散 219 s = sqrt(v); // 標準偏差 220 221 double g=2.75*pow(10,-3); 222 for (int i=0;i<L;i++){ 223 E[i]=1/b_pi/pow(b_pi,0.5)/s*exp(-pow((e-A[i]),2)/2/v); 224 sum_density[r] += E[i]; 225 } 226 227 //printf("%7.4lf\n",sum_density); 228 229 /*for (int i = 0;i < L;i++) { 230 fprintf(fp, "%2f %lf\n", k_x, A[i]); 231 }*/ 232 233 fprintf(fp, "%2f %lf %lf \n",k_x,e,sum_density[r]); 234 235 236 } 237 238 double f=0.0; 239 f+=sum_density[h]; 240 for(int i=0;i<n;i++){ 241 sum_dos+=sum_density[i]; 242 } 243 244 245 //printf("%lf\n",sum_dos); 246 //fprintf(fp, "%2f %lf %lf %lf %lf\n",k_x,e,sum_density[h],sum_dos,f); 247 } 248 249 fclose(fp); 250 return 0; 251} 252
試したこと
whileにしたり、外のfor文にしたりしましたがあまりうまく行かなかったです。
また物理の状態密度を求めるソースコードとかがあれば教えていただきたいです。自分でコードを書くことにはあまり固執していないです。
補足情報(FW/ツールのバージョンなど)
ここにより詳細な情報を記載してください。
>eが100回増加を繰り返した後、k_xは100の二乗回増加されています
:
>eが一回増加ごとにk_xが100回増加してそれごとにsum_densityを計算
:
>whileにしたり、外のfor文にしたりしましたが
:
テレビを叩けば直るみたいなやり方で計算結果が良くなるはずはありません。
コンピュータは石頭です。書いたコードの通りにしか動きません。
現状がそのように動いているのであればそのようなコードを書いているからで、希望する動作があるのであればそのようなコードに書き直せば良いだけです。
どこで計算が思惑通りにならなくなっているのか、各変数を表示して机上の計算と合わせて間違いを探し出す等を行ってください。
> これをeが一回増加ごとにk_xが100回増加してそれごとにsum_densityを計算して。。。
何を希望しているのか不明瞭です。eの増加以降のコードを別の関数にしたい、ということでしょうか?
何にせよ、希望だけではなく質問を書きましょう。
回答1件
あなたの回答
tips
プレビュー