/*　フォッカープランクプログラム */

#include <stdio.h>

main(){

/*ファイル宣言*/

FILE *indata1_file ; /* 入力ファイル１宣言（データdt1の入力用）*/
FILE *indata2_file ; /* 入力ファイル２宣言（データdt2の入力用）*/
FILE *outdata_file ; /* 出力ファイル１宣言（行列とその解の出力用） */
FILE *outdata2_file ; /* 出力ファイル２宣言（データdt3の出力用）*/

/*変数宣言*/

int m,n ; /* データ数用変数宣言 */
int i,j,k ; /* ループ用変数宣言 */

double e[100][100],f[100][100],g[100][100],h[100][100] ;/* 行列要素配列宣言 */
double T[100][100] ;/* 行列要素配列宣言 */

double dt1[100][100],dt2[100][100] ; /* 実況データ配列宣言 */
double dt3[100][100] ; /* 予測データ配列宣言 */

double e2,ef,eg,eh ; /* 行列要素宣言 */
double f2,fg,fh ; /* 行列要素宣言 */
double g2,gh ; /* 行列要素宣言 */
double h2 ; /* 行列要素宣言 */

double eT,fT,gT,hT ; /* ベクトル要素宣言 */

/* データ数の読み込み*/

printf("縦横のデータ数を入力してください。\n") ;

printf("横のデータ数はいくつですか？"); scanf("%d",&m) ; /* 横データ数読み込み */
printf("縦のデータ数はいくつですか？"); scanf("%d",&n) ; /* 縦データ数読み込み */

/* 予測回数の選択 */

int sle ; /* 予測回数選択用変数宣言 */

printf("何回目の予測ですか？1to10") ; scanf("%d", &sle) ; /* 予測回数の読み込み */

switch(sle){

case 1 :

printf("予測１回目計算開始 \n") ;
indata1_file = fopen("adt1.txt","r") ; /* 読み込みファイル１宣言*/
indata2_file = fopen("adt2.txt","r") ; /* 読み込みファイル２宣言*/
outdata_file = fopen("matrixdata3.txt","w") ; /* 書き込みファイル１宣言*/
outdata2_file = fopen("f4dt3.txt","w") ; /* 書き込みファイル２宣言*/

break ;

case 2 :

printf("予測２回目計算開始 \n") ;
indata1_file = fopen("adt2.txt","r") ;
indata2_file = fopen("f4dt3.txt","r") ;
outdata_file = fopen("matrixdata4.txt","w") ;
outdata2_file = fopen("f4dt4.txt","w") ;

break ;

/* 中略 */

/*

予測計算を(X)回行うためには、この部分に次の行を足してください。
張り付けた後は、(X)の部分（括弧部分も含む）を予測する回数の数字（半角）に直してください。

(X)＝1,2,3,4,5・・・ｎ

-------------以下から-------------

case (X) :

printf("予測(X)回目計算開始 \n") ;
indata1_file = fopen("adt(X).txt","r") ;
indata2_file = fopen("f4dt(X+1).txt","r") ;
outdata_file = fopen("matrixdata(X+2).txt","w") ;
outdata2_file = fopen("f4dt(X+2).txt","w") ;

break ;

-------------以上まで-------------

*/

case 10 :

printf("予測１０回目計算開始 \n") ;
indata1_file = fopen("f4dt10.txt","r") ;
indata2_file = fopen("f4dt11.txt","r") ;
outdata_file = fopen("matrixdata12.txt","w") ;
outdata2_file = fopen("f4dt12.txt","w") ;

break ;

default :

printf("予測１回目計算開始 \n") ;
indata1_file = fopen("adt1.txt","r") ;
indata2_file = fopen("adt2.txt","r") ;
outdata_file = fopen("matrixdata3.txt","w") ;
outdata2_file = fopen("f4dt3.txt","w") ;

break ;

}

printf("計算開始。\n") ; /* 計算中のコメント */

/* データ読み込み */

printf("データ読み込み中・・・\n") ; /* データ読み込みのコメント */

for(j=0; j<=n-1; j++){ /* 列データ */

for(i=0; i<=m-1; i++){ /* 行データ */

fscanf(indata1_file,"%lf",&dt1[i][j]) ;/* 実況データ１読み込み */
fscanf(indata2_file,"%lf",&dt2[i][j]) ;/* 実況データ２読み込み */

}

}

/* 最小自乗法用要素計算 */

printf("最小自乗法−行列要素計算中・・・\n") ; /* 計算中のコメント */

for(j=1; j<=n-2; j++){ /*列データ*/

for(i=1; i<=m-2; i++){ /*行データ*/

e[i][j] = dt1[i-1][j] - dt1[i][j] ; /* 要素 e 計算 */
f[i][j] = dt1[i+1][j] - dt1[i][j] ; /* 要素 f 計算 */
g[i][j] = dt1[i][j-1] - dt1[i][j] ; /* 要素 g 計算 */
h[i][j] = dt1[i][j+1] - dt1[i][j] ; /* 要素 h 計算 */

T[i][j] = dt2[i][j] - dt1[i][j] ; /* 要素 T 計算 */

/* 各行列要素の積算 */

e2 += e[i][j]*e[i][j] ;
ef += e[i][j]*f[i][j] ;
eg += e[i][j]*g[i][j] ;
eh += e[i][j]*h[i][j] ;

f2 += f[i][j]*f[i][j] ;
fg += f[i][j]*g[i][j] ;
fh += f[i][j]*h[i][j] ;

g2 += g[i][j]*g[i][j] ;
gh += g[i][j]*h[i][j] ;

h2 += h[i][j]*h[i][j] ;

eT += e[i][j]*T[i][j] ;
fT += f[i][j]*T[i][j] ;
gT += g[i][j]*T[i][j] ;
hT += h[i][j]*T[i][j] ;

}

}

/*
* 4*4行列
*
* ｜e2　ef　eg　eh ｜｜a｜＝｜eT｜
* ｜ 　　f2　fg　fh　｜｜b｜＝｜fT｜
* ｜　 　　　g2　gh ｜｜c｜＝｜gT｜
* ｜ * 　　 　　　h2 ｜｜d｜＝｜hT｜
*
* をGauss消去法によって解く
*
*/

printf("行列の計算中・・・\n") ; /* 計算中のコメント */　

/* 行列、ベクトル要素宣言 */

double matrixA[4][4],vectorB[4] ; /* 行列要素、ベクトル要素宣言 */

double valueX[4] ; /* 解ベクトル要素宣言 */

/* 行列A 各要素代入 */

matrixA[0][0] = e2 ;
matrixA[0][1] = ef ;
matrixA[0][2] = eg ;
matrixA[0][3] = eh ;

matrixA[1][1] = f2 ;
matrixA[1][2] = fg ;
matrixA[1][3] = fh ;

matrixA[2][2] = g2 ;
matrixA[2][3] = gh ;

matrixA[3][3] = h2 ;

/* ベクトルB 要素代入*/

vectorB[0] = eT ;
vectorB[1] = fT ;
vectorB[2] = gT ;
vectorB[3] = hT ;

/* 元の行列、ベクトルを出力データファイル１に保存する */

for(i=0; i<=3; i++){

for(j=i; j<=3; j++){ /* 対角行列のためｊはｉから始まる*/

fprintf(outdata_file,"a(%d,%d)=%lf ",i,j,matrixA[i][j]) ;
/* 行列Aの書き込み */

}

fprintf(outdata_file,"b(%d)=%lf\n",i,vectorB[i]) ; /* ベクトルBの書き込み */

}

/* Gauss消去実行*/

for(i=0; i<=2; i++){

for(j=i+1; j<=3; j++){ /* 対角行列のためｊはｉ＋１から始まる */

for(k=j; k<=3; k++){ /* 対角行列のためｋはｊから始まる */

matrixA[j][k] = matrixA[i][i] * matrixA[j][k] - matrixA[i][j] * matrixA[i][k] ;
/* 各行列要素a(j,k)を消去 */

}

vectorB[j] = matrixA[i][i] * vectorB[j] - matrixA[i][j] * vectorB[i] ;
/* 各ベクトル要素b(j)を消去 */

}

}

/* 積算用変数宣言 */

double nsum=0.0,nsumX=0.0 ;

/* 求める解Xを計算 */

for(i=3; i>=0; i--){ /* 後退しながら解Xを計算 */

nsum = 0.0 ; /* nsumの初期化 */

for(j=i+1; j<=3; j++){

nsum += matrixA[i][j] * valueX[j] ;
/* 前の解Xと要素aの乗を積算　そのときjはi+1から始まる（後退）*/

}

valueX[i] = (vectorB[i] - nsum) / matrixA[i][i]; /* 解Xを計算 */

}

/* 書き込み中のコメント */

printf("解をファイル“matrixdata.txt”に書き込み中・・・\n") ;

/* 解Xを出力データ１に書き込む */

for(i=0; i<=3; i++){

fprintf(outdata_file,"x[%d]=%lf\n",i,valueX[i]) ; /* 解Xを書き込み */

nsumX += valueX[i] ; /* 解Xを積算 */

}

nsumX = 1.0 - nsumX ; /* 積算した解Xを１から引く（計算を簡単化するため）*/

/* dt3予測計算実行 */

printf("予測値の計算中・・・\n") ; /* 計算中のコメント */

/* (i,j)=(1,1)〜(m-2,n-2)を計算 */

for(j=1; j<=n-2; j++){

for(i=1; i<=m-2; i++){

dt3[i][j] = nsumX * dt2[i][j] ;
dt3[i][j] += valueX[0] * dt2[i-1][j] ;
dt3[i][j] += valueX[1] * dt2[i+1][j] ;
dt3[i][j] += valueX[2] * dt2[i][j-1] ;
dt3[i][j] += valueX[3] * dt2[i][j+1] ;

}

}

/* ４辺の計算 */

/* (i,j)=(1,0)〜(m-2,0)と(1,n-1)〜(m-2,n-1)を計算 */

for(i=1; i<=m-2; i++){

dt3[i][0] = dt3[i][1] ;
dt3[i][n-1] = dt3[i][n-2] ;

}

/* (i,j)=(0,1)〜(0,n-2)と(m-1,1)〜(m-2,n-2)を計算 */

for(j=1; j<=n-2; j++){

dt3[0][j] = dt3[1][j] ;
dt3[m-1][j] = dt3[m-2][j] ;

}

/* ４隅の計算 */

/* dt3(0,0)の計算 */

dt3[0][0] = (dt3[1][0] + dt3[0][1]) / 2 ;

/* dt3(0,n-1)の計算 */

dt3[0][n-1] = (dt3[1][n-1] + dt3[0][n-2]) / 2 ;

/* dt3(m-1,0)の計算 */

dt3[m-1][0] = (dt3[m-2][0] + dt3[m-1][1]) / 2 ;

/* dt3(m-1,n-1)の計算 */

dt3[m-1][n-1] = (dt3[m-2][n-1] + dt3[m-1][n-2]) / 2 ;

/* dt3(i､j)の予測結果出力 */
/* 書き込み中のコメント */

printf("予測結果をファイル“f4dtn.txt”に書き込み中・・・\n") ;

for(j=0; j<=n-1; j++){

for(i=0; i<=m-1; i++){

fprintf(outdata2_file,"%lf ",dt3[i][j]) ;

}

fprintf(outdata2_file,"\n") ;

}

printf("計算終了\n");

/* ファイルクローズ */

fclose(indata1_file) ;
fclose(indata2_file) ;
fclose(outdata_file) ;
fclose(outdata2_file) ;　

return(0) ;

}