質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

Q&A

解決済

1回答

564閲覧

繰り返しを独立させたい

kisara11235

総合スコア18

C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

0グッド

0クリップ

投稿2022/11/08 02:56

前提

状態密度を計算したい。

実現したいこと

そのために、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/ツールのバージョンなど)

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

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

jimbe

2022/11/08 04:27

>eが100回増加を繰り返した後、k_xは100の二乗回増加されています : >eが一回増加ごとにk_xが100回増加してそれごとにsum_densityを計算 : >whileにしたり、外のfor文にしたりしましたが : テレビを叩けば直るみたいなやり方で計算結果が良くなるはずはありません。 コンピュータは石頭です。書いたコードの通りにしか動きません。 現状がそのように動いているのであればそのようなコードを書いているからで、希望する動作があるのであればそのようなコードに書き直せば良いだけです。 どこで計算が思惑通りにならなくなっているのか、各変数を表示して机上の計算と合わせて間違いを探し出す等を行ってください。
int32_t

2022/11/08 04:38

> これをeが一回増加ごとにk_xが100回増加してそれごとにsum_densityを計算して。。。 何を希望しているのか不明瞭です。eの増加以降のコードを別の関数にしたい、ということでしょうか? 何にせよ、希望だけではなく質問を書きましょう。
guest

回答1

0

ベストアンサー

直接的な答えではありません。
指摘事項として、同じ変数名(ループカウンタ)を使いまわさないこと。
どの同名の変数名を配列アクセスに使用しているのでかなり危ない。
紛らわしく読み間違える可能性が高い。よくエラーが出ないなと思う。

本当に意図通りに値が変化しているか見てください。

最後の方でfを算出しているが何にも使っていないので無意味です。

投稿2022/11/08 08:13

編集2022/11/08 08:16
ardin

総合スコア544

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問