// -----------------------------------------------|
// parabora-using-sub2.cpp |
// |
// FDTD法による直交座標系(x,y)における |
// パラボラミラーによる電磁波の反射の数値計算 |
// (動作確認済み) |
// |
// プログラム作成：駒田寛 2010年6月 |
// -----------------------------------------------|
#include <vector>
#include <complex>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
// 2次元配列をフリーストアに作成するサブルーチン
// ＜実数配列用＞
double** makeArray(int row, int col)
{
double** a=new double*[row]; //全体のための配列を作る
for(int i=0; i<row; i++){
a[i]=new double[col]; //各行のための配列を作る
}
return a;
}
// 上で作った実数配列の削除
void deleteArray(double** a, int row){
for (int i=0;i<row; ++i) delete[] a[i];
delete[] a;
}
// 2次元配列をフリーストアに作成するサブルーチン
// ＜整数配列用＞
int** makeArrayInt(int row, int col){
int** a=new int*[row]; //全体のための配列を作る
for(int i=0; i<row; i++){
a[i]=new int[col]; //各行のための配列を作る
}
return a;
}
// 上で作った整数配列の削除
void deleteArrayInt(int** a, int row){
for (int i=0;i<row; ++i) delete[] a[i];
delete[] a;
}
// 電磁界が既知である部分のマップを作成する
void knownMap(char* str, char* str2, int** KNOWN, int NX, int NY){
cout << setw(30) << str << endl;
ofstream fout02(str2); //「やさしいC++」16章
cout << "ファイルをオープンしました。
";
for(int AI=1;AI<=NX;AI++){ //間引きせずに書き込む
for(int J=1;J<=NY;J++){
fout02 << setw(6) << KNOWN[AI][J];
}
fout02 << endl;
}
cout << "ファイルに書き込みました。
";
fout02.close();
cout << "ファイルをクローズしました。
";
cout << "何かキーを押してください";
getchar();
}
// 2次元実数配列の各種出力データを作成する
void kakikomi(char* str, char* str2, double** a, int NX, int NY){
cout << str << endl;
ofstream fout05(str2);
cout << " " << "ファイルをオープンしました。" << endl;
for(int AI=1;AI<=NX;AI++){ //間引きせずに書き込む
for(int J=1;J<=NY;J++){
fout05 << " " << scientific << a[AI][J];
}
fout05 << endl;
}
cout << " " << "ファイルに書き込みました。" << endl;
fout05.close();
cout << " " << "ファイルをクローズしました。" << endl;
}
//その他の出力データ
//一点における時間変化
void onePointMonitor(char* str, char* str2, double* EZ_MON_1, int ITMAX ){
cout << str << endl;
ofstream fout10(str2);
cout << " " << "ファイルをオープンしました。" << endl;
for(int K=1;K<=ITMAX;K++){
fout10 << " " << scientific << EZ_MON_1[K] << endl;
}
cout << "__ファイルに書き込みました。" << endl;
fout10.close();
cout << "__ファイルをクローズしました。" << endl;
}
//定数データ用クラス
struct ConstData{
char* myname;
double myvalue;
void show(); //メソッドの宣言
};
//メソッドの定義 定数データを画面とファイルの両方に出力する
void ConstData::show()
{
ofstream fout01("20100729_01_SET_DATA_CONST.txt",ios_base::app);
cout << setw(30) << myname << ":" << myvalue << endl;
fout01 << setw(30) << myname << ":" << myvalue << endl;
fout01.close();
}
//--- 1-0.配列宣言------------------------------------------
//要素を１行、1列それぞれ余分に宣言しておく。
//そして、添え字０の要素はデータとして使わない。
//（fortranプログラムとの兼ね合いである）
//フリーストアに格納する配列（日経BP 「VisualC++入門」）
double** CEZ = makeArray(1002,3002); //2次元配列の行数、列数を渡す
double** CEZLX = makeArray(1002,3002);
double** CEZLY = makeArray(1002,3002);
double** CHXLY = makeArray(1002,3002);
double** CHYLX = makeArray(1002,3002);
double** MEDIUM = makeArray(1002,3002);
double** C = makeArray(1002,3002); //導電率
double** DA = makeArray(1002,3002);
double** CV = makeArray(1002,3002); //媒質中の光速
int** KNOWN = makeArrayInt(1002,3002);
double** HX = makeArray(1002,3002);
// double** HX1 = makeArray(1002,3002);
double** HY = makeArray(1002,3002);
// double** HY1 = makeArray(1002,3002);
double** EZ = makeArray(1002,3002);
// double** EZ1 = makeArray(1002,3002);
double* EZ_MON_1=new double[8005]; //ここは1次元配列
double** EZX1 = makeArray(5,3002); //4+1行（後ろ参照）
double** EZX2 = makeArray(5,3002); //4+1行（後ろ参照）
double** EZY1 = makeArray(1002,5); //4+1列（後ろ参照）
double** EZY2 = makeArray(1002,5); //4+1列（後ろ参照）
double** AMPEZ= makeArray(1002,3002);
//自動メモリに格納する配列
int NX;
int NY;
double PAI;
double E;
double MU;
double F;
double W;
double DT;
double DX,DY;
double THETA;
double THETA1,THETA2,THETA3,THETA4;
double CV0; //真空中の光速(シーブイゼロ)
int K; //時間ステップ数格納用
int J; //Y方向ステップ数格納用
int AI; //X方向ステップ数格納用
int JDUM,JDUM2; //仮座標用変数
int AIDUM; //仮座標用変数
int IT;
int ITMAX;
double T,Tw;
double C1,C2,C3;
//本のCXD,CXX,CXFYDをそれぞれC1,C2,C3と置いた。
//２次元ではCXFZDはない。
//--- 1-0.配列宣言 ここまで-----------------------------------------------
/* 定数設定 */
void const_set(void){
//--- 1-1．空間、時間に対して定数（物理定数も含む）であるものの設定----------
//全ステップ回数
ITMAX = 12000; //
// x,y方向の分割数+1（EZノードの数）の値
NX=1001; // 吸収境界の座標はXは1,1001とする //いまだけ10分の１
NY=3001; // 吸収境界の座標はYは1,2001とする //いまだけ10分の１
//円周率
PAI = 3.141592654e0;
//真空中の誘電率 [F/m]
E = 8.8541878e-12;
//真空中の透磁率 [H/m]
MU = 4 * PAI * 1.0e-7;
//調べる電磁波の周波数［Hz］
F = 15.0e+9;
//2周期分の時間 [s]
Tw = 2/F;
//電磁波の角周波数 [rad/s]
W = F * 2 * PAI;
//時間ステップ［s］
DT = 1.0e-12;
// 注１ 最適時間ステップの計算がEXCELおよびメモ帳にある
//X方向のセルサイズ［m］
DX = 1.0e-3;
//Y方向のセルサイズ [m]
DY = 1.0e-3;
// 注２ 最適セルサイズの計算がEXCELおよびメモ帳にある
CV0 = 1/pow((E*MU) , 0.5); //真空中の光速
C1 = (CV0*DT-DX)/(CV0*DT+DX);
C2 = 2*DX/(CV0*DT+DX);
C3 = DX * pow((CV0*DT),2) / ( 2 * pow((DY),2) * (CV0*DT+DX) );
//変数は用意してあるが本プログラムでは常にDX=DYとする。
//画面での確認、ファイルへの出力
cout << "定数設定をSET_DATA_CONSTファイルに出力します。" << endl;
// "20100729_01_SET_DATA_CONST.txt"
ConstData aa0;
aa0.myname = "X方向のセルサイズ [m]=";
aa0.myvalue = DX;
aa0.show();
ConstData aa1;
aa1.myname = "Y方向のセルサイズ [m]=";
aa1.myvalue = DY;
aa1.show();
ConstData aa2;
aa2.myname = "時間ステップ［s］=";
aa2.myvalue = DT;
aa2.show();
ConstData aa3;
aa3.myname = "１秒当たりのステップ回数 =";
aa3.myvalue = 1/DT;
aa3.show();
ConstData aa4;
aa4.myname = "真空中の誘電率 [F/m]=";
aa4.myvalue = E;
aa4.show();
ConstData aa5;
aa5.myname = "調べる電磁波の周波数［Hz］=";
aa5.myvalue = F;
aa5.show();
ConstData aa6;
aa6.myname = "１周期当たりのステップ回数 =";
aa6.myvalue = 1/DT/F;
aa6.show();
ConstData aa7;
aa7.myname = "真空中の透磁率 [H/m]=";
aa7.myvalue = MU;
aa7.show();
ConstData aa8;
aa8.myname = "X方向の分割数＋１の値 =";
aa8.myvalue = NX;
aa8.show();
ConstData aa9;
aa9.myname = "Y方向の分割数＋１の値 =";
aa9.myvalue = NY;
aa9.show();
ConstData aa10;
aa10.myname = "円周率 =";
aa10.myvalue = PAI;
aa10.show();
ConstData aa11;
aa11.myname = "電磁波の角周波数 [rad/s] =";
aa11.myvalue = W;
aa11.show();
ConstData aa12;
aa12.myname = "真空中の光速[m/s] =";
aa12.myvalue = CV0;
aa12.show();
ConstData aa13;
aa13.myname = "全ステップ回数 =";
aa13.myvalue = ITMAX;
aa13.show();
cout << "何かキーを押してください" << endl;
getchar();
//--- 1-1．空間、時間に対して定数であるものの設定 ここまで-------------------
}
//メイン関数
int main(void)
{
//=== 1. setup ==================================================================
// 宇野先生の本で「 setup 」となっている部分とほぼ同じである
//---1. setup 目次-------------------------------------------------------
// 1-0. 配列宣言
// 1-1．空間、時間に対して定数であるものの設定
// 1-2．空間に対して変数、時間に対して定数であるものの設定
// 1-3．時間に対して変数であるものの初期値設定
// 1-4. メモデータ作成（本にはない）
//---1. setup 目次 ここまで---------------------------------------------
//--- 1-1．空間、時間に対して定数（物理定数も含む）であるものの設定----------
const_set();
//--- 1-2．空間的には変数、時間的には定数であるものの設定 -------------------
// 1-1に該当するものもあるが、揺れ動くものはすべてこちらにした。
//使わない添え字0の要素も念のため初期化する
for(J=0;J<=NY;J++){
for(AI=0;AI<=NX;AI++){
CEZLX[AI][J]=0.0; //EZを求める式に使用する係数
CEZLY[AI][J]=0.0; //EZを求める式に使用する係数
CHXLY[AI][J]=0.0; //HXを求める式に使用する係数
CHYLX[AI][J]=0.0; //HYを求める式に使用する係数
CEZ[AI][J]=0.0; //EZを求める式に使用する係数
C[AI][J]=0.0; //導電率
//金属部分については電界＝0として別の場所で表現する。
//電界を隣の値から求めるときにCを使うが
//金属部分についてはそれをやらないので値は関係ない
//(0のままでも良い)。
KNOWN[AI][J]=0;
}
}
for(J=0;J<=NY;J++){
for(AI=0;AI<=NX;AI++){
MEDIUM[AI][J]=0.0;
CV[AI][J]=0.0;
DA[AI][J]=0.0;
CEZ[AI][J]=0.0;
CEZLX[AI][J]=0.0;
CEZLY[AI][J]=0.0;
CHXLY[AI][J]=0.0;
CHYLX[AI][J]=0.0;
}
}
for(J=1;J<=NY;J++){
for(AI=1;AI<=NX;AI++){
MEDIUM[AI][J]=1.0; //比誘電率（今回は誘電体はない）
CV[AI][J] = 1/pow((MEDIUM[AI][J]*E*MU),0.5); //媒質中の光速 [m/s]
DA[AI][J]=C[AI][J]*DT/(2 * E * MEDIUM[AI][J]); // 導電率0の部分では0になる
CEZ[AI][J]=(1 -DA[AI][J])/(1+DA[AI][J]); // 導電率0の部分では1になる
CEZLX[AI][J]=DT/E/MEDIUM[AI][J]/(1+DA[AI][J])/DX;
CEZLY[AI][J]=DT/E/MEDIUM[AI][J]/(1+DA[AI][J])/DY;
CHXLY[AI][J]=DT/MU/DY;
//P13(1.31)にある。磁性体があれば比誘電率を掛けた値をMUとして用いる
CHYLX[AI][J]=DT/MU/DX; //同上
}
}
//波源の設定（既知の部分は1回も微分方程式で解かないにする)
//（ KNOWN の値 ＝ 既知で0：１、既知で時間変動がある：２、未知：０）
//図形回転用角度（実際 マイナス ダミー）[degree] の準備
THETA = 126.9e0;
THETA1 = cos(THETA*PAI/180);
THETA2 = -sin(THETA*PAI/180);
THETA3 = sin(THETA*PAI/180);
THETA4 = cos(THETA*PAI/180);
for(AIDUM=433;AIDUM<=475;AIDUM++){ //1次放射器（金属板）
JDUM=467;
//図形の回転操作（(500,467)を中心としてTHETA度回転する）
AI = floor( (AIDUM-500)*THETA1 +(JDUM-467)*THETA2 +500 );
J = floor( (AIDUM-500)*THETA3 +(JDUM-467)*THETA4 +467 );
KNOWN[AI][J]=1;
}
for(AIDUM=525;AIDUM<=567;AIDUM++){ //1次放射器（金属板）
JDUM=467;
//図形の回転操作（(500,467)を中心としてTHETA度回転する）
AI = floor( (AIDUM-500)*THETA1 +(JDUM-467)*THETA2 +500 );
J = floor( (AIDUM-500)*THETA3 +(JDUM-467)*THETA4 +467 );
KNOWN[AI][J]=1;
}
for(AIDUM=567;AIDUM<=667;AIDUM++){ //遮蔽器（金属板）
JDUM=1034-AIDUM;
//図形の回転操作（(500,467)を中心としてTHETA度回転する）
AI = floor( (AIDUM-500)*THETA1 +(JDUM-467)*THETA2 +500 );
J = floor( (AIDUM-500)*THETA3 +(JDUM-467)*THETA4 +467 );
KNOWN[AI][J]=1;
KNOWN[AI][J+1]=1;
}
for(JDUM=1;JDUM<=367;JDUM++){ //遮蔽器（金属板）
AIDUM=667;
JDUM2=JDUM+0;
//図形の回転操作（(500,467)を中心としてTHETA度回転する）
AI = floor( (AIDUM-500)*THETA1 +(JDUM2-467)*THETA2 +500 );
J = floor( (AIDUM-500)*THETA3 +(JDUM2-467)*THETA4 +467 );
KNOWN[AI][J]=1;
}
for(AI=772;AI<=996;AI++){ //追加遮蔽器（金属板）
J=880;
KNOWN[AI][J]=1;
}
for(AIDUM=485;AIDUM<=515;AIDUM++){ //1次放射器（余弦波波源部分）
JDUM=400;
//図形の回転操作（(500,467)を中心としてTHETA度回転する）
AI = floor( (AIDUM-500)*THETA1 +(JDUM-467)*THETA2 +500 );
J = floor( (AIDUM-500)*THETA3 +(JDUM-467)*THETA4 +467 );
KNOWN[AI][J]=2;
}
//パラボラミラーを設置する
for(AI=201;AI<=401;AI++){
int ss=(int)pow((double)(AI-500),2)/800+267;
KNOWN[AI][ss]=1;
//下に凸の放物線 （ 頂点（500、267）、焦点（500、467）のもの ）上にとる
}
// 心配ならこのように書く
// IF(1<=AI .and. AI<=1001 .and. 1<=J .and. J<=1001)then
// KNOWN[AI][J]=1
// ENDIF
//--- 1-2．空間的には変数、時間的には定数であるものの設定
//ここまで----------------------------------------
//--- 1-3．時間に対して変数であるものの初期値設定--------------------------
//--- 1-3-1.電磁場の初期値の設定 --------------------------------------
//現時点での値
for(J=0;J<=NY;J++){
for(AI=0;AI<=NX;AI++){
HX[AI][J]=0.0; //磁場のＸ成分（0とする、以下同様）
HY[AI][J]=0.0; //磁場のＹ成分
EZ[AI][J]=0.0; //電場のＺ成分
}
}
// //過去の値1
// //1ターン前のEZ,HX,HY
// DO AI=1,NX
// DO J=1,NY
// EZ1[AI][J]=0.0e0
// HX1[AI][J]=0.0e0
// HY1[AI][J]=0.0e0
// ENDDO
// ENDDO
//過去の値2
//AI=1,NXの辺とその隣の辺（計4本の線）の1ターン、２ターン前のEZ
//（添え字0を含めて5本）
for(AI=0;AI<=4;AI++){
for(J=0;J<=NY;J++){
EZX1[AI][J]=0.0;
EZX2[AI][J]=0.0;
}
}
//過去の値3
//J=1,NYの辺とその隣の辺（計4本の線）の1ターン、２ターン前のEZ
for(J=0;J<=4;J++){
for(AI=0;AI<=NX;AI++){
EZY1[AI][J]=0.0;
EZY2[AI][J]=0.0;
}
}
//--- 1-3-1.電磁場の初期値の設定 ここまで------------------------------
// //--- 1-3-2.１ポイントモニタ用配列 EZ_MON の初期化---------------------
//
// DO K=1,10000
// EZ_MON_1(K)=0.0e-0
// ENDDO
//
// //--- 1-3-2.１ポイントモニタ用配列 EZ_MON の初期化 ここまで------------
//--- 1-3-3. 経過時間 T の初期化 --------------------------------------
T=0.0;
//--- 1-3-3. 経過時間 T の初期化 ここまで------------------------------
//--- 1-3．時間に対して変数であるものの初期値設定 ここまで----------------
//--- 1-4. メモデータを作る(本にはない)------------------------------------
knownMap("既知部分KNOWN（AI,J)ファイルを作成します。" ,
"20100729_01_SET_DATA_KNOWN.txt" , KNOWN , NX , NY);
//--- 1-4. メモデータを作る ここまで--------------------------------------
//=== 1. setup ここまで ======================================================
//=== 2. メインループ ========================================================
//--- 2-1. メインループの始まり--------------------------------------------
for(IT=0;IT<=ITMAX;IT++){ // 計算のスタート
// 経過時間＝(IT)＊DT
//--- 2-1. メインループの始まり ここまで------------------------------------
//--- 2-2. 電磁場のEZの計算(本では 「 CALL E-field 」)---------------------
//(P15にある)
//AI,Jともに内側を全てスキャン
//---2-2-1. 余弦波波源部分-----------------------------------------
for(AIDUM=485;AIDUM<=515;AIDUM++){ //1次放射器（余弦波波源部分）
JDUM=400;
//図形の回転操作（(500,467)を中心としてTHETA度回転する）
AI = floor( (AIDUM-500)*THETA1 +(JDUM-467)*THETA2 +500 );
J = floor( (AIDUM-500)*THETA3 +(JDUM-467)*THETA4 +467 );
if(-0.1<=T && T<=Tw){
EZ[AI][J]=1.0/2.0*(1.0-cos(PAI*T/Tw))*sin(W*T);
}else{
EZ[AI][J]=3.0*cos(W*T);
}
}
// //---2-2-1. 余弦波波源部分 ここまで--------------------------------
//---2-2-2. 本体の計算 --------------------------------------------------
for(J=2;J<=NY-1;J++){ //両端以外全部
for(AI=2;AI<=NX-1;AI++){ //同上
if(KNOWN[AI][J]==0){
EZ[AI][J]=CEZ[AI][J]*EZ[AI][J]+CEZLX[AI][J]*(HY[AI][J]-HY[AI-1][J])
-CEZLY[AI][J]*(HX[AI][J]-HX[AI][J-1]);
}
}
}
//---2-2-2. 本体の計算 ここまで------------------------------------------
//--- 2-2. CALL E-field ここまで-------------------------------------------
//--- 2-3. 境界（周辺）の設定(本では「 CALL A.B.C. 」)---------------------
//吸収境界条件(P60)
//４辺のすぐ近くに誘電体はないとして式を立てる
//---2-3-1. AI=1,NXの辺での計算 ----------------------------------------
//本体の計算
J=2;
EZ[1][J]=EZX1[2][J]+C1*(EZ[2][J]-EZX1[1][J]);
EZ[NX][J]=EZX1[3][J]+C1*(EZ[NX-1][J]-EZX1[4][J]);
J=NY-1;
EZ[1][J]=EZX1[2][J]+C1*(EZ[2][J]-EZX1[1][J]);
EZ[NX][J]=EZX1[3][J]+C1*(EZ[NX-1][J]-EZX1[4][J]);
for(J=3;J<=NY-2;J++){
EZ[1][J]=-EZX2[2][J]+C1*(EZ[2][J]+EZX2[1][J])
+C2*(EZX1[1][J]+EZX1[2][J])
+C3*(EZX1[1][J+1]-2.0e-0*EZX1[1][J]+EZX1[1][J-1]
+EZX1[2][J+1]-2.0e-0*EZX1[2][J]+EZX1[2][J-1]);
EZ[NX][J]=-EZX2[3][J]+C1*(EZ[NX-1][J]+EZX2[4][J])
+C2*(EZX1[4][J]+EZX1[3][J])
+C3*(EZX1[4][J+1]-2.0e-0*EZX1[4][J]+EZX1[4][J-1]
+EZX1[3][J+1]-2.0e-0*EZX1[3][J]+EZX1[3][J-1]);
}
//過去の値の更新
for(J=2;J<=NY-1;J++){
EZX2[1][J]=EZX1[1][J];
EZX2[2][J]=EZX1[2][J];
EZX2[3][J]=EZX1[3][J];
EZX2[4][J]=EZX1[4][J];
EZX1[1][J]=EZ[1][J];
EZX1[2][J]=EZ[2][J];
EZX1[3][J]=EZ[NX-1][J];
EZX1[4][J]=EZ[NX][J];
}
//---2-3-1. AI=1,NXの辺での計算 ここまで-------------------------------
//---2-3-2. J=1,NYの辺での計算 ----------------------------------------
//本体の計算
AI=2;
EZ[AI][1]=EZY1[AI][2]+C1*(EZ[AI][2]-EZY1[AI][1]);
EZ[AI][NY]=EZY1[AI][3]+C1*(EZ[AI][NY-1]-EZY1[AI][4]);
AI=NX-1;
EZ[AI][1]=EZY1[AI][2]+C1*(EZ[AI][2]-EZY1[AI][1]);
EZ[AI][NY]=EZY1[AI][3]+C1*(EZ[AI][NY-1]-EZY1[AI][4]);
for(AI=3;AI<=NX-2;AI++){
EZ[AI][1]=-EZY2[AI][2]+C1*(EZ[AI][2]+EZY2[AI][1])
+C2*(EZY1[AI][1]+EZY1[AI][2])
+C3*(EZY1[AI+1][1]-2.0e-0*EZY1[AI][1]+EZY1[AI-1][1]
+EZY1[AI+1][2]-2.0e-0*EZY1[AI][2]+EZY1[AI-1][2]);
EZ[AI][NY]=-EZY2[AI][3]+C1*(EZ[AI][NY-1]+EZY2[AI][4])
+C2*(EZY1[AI][4]+EZY1[AI][3])
+C3*(EZY1[AI+1][4]-2.0e-0*EZY1[AI][4]+EZY1[AI-1][4]
+EZY1[AI+1][3]-2.0e-0*EZY1[AI][3]+EZY1[AI-1][3]);
}
//過去の値の更新
for(AI=2;AI<=NX-1;AI++){
EZY2[AI][1]=EZY1[AI][1];
EZY2[AI][2]=EZY1[AI][2];
EZY2[AI][3]=EZY1[AI][3];
EZY2[AI][4]=EZY1[AI][4];
EZY1[AI][1]=EZ[AI][1];
EZY1[AI][2]=EZ[AI][2];
EZY1[AI][3]=EZ[AI][NY-1];
EZY1[AI][4]=EZ[AI][NY];
}
//--- 2-3. 境界（周辺）の設定(本では「 CALL A.B.C. 」) ここまで------------
//--- 2-4. 半時間ステップだけ時間を進める------------------------------------
T=T+DT/2;
// WRITE(6,'(a,e15.7)')'T=',T //画面での確認
//--- 2-4. 半時間ステップだけ時間を進める ここまで---------------------------
// --- 2-5. 電磁場の磁場成分HX,HYの計算(「 CALL H-field 」)----------------
//HXの計算（P15)
for(J=1;J<=NY-1;J++){
for(AI=2;AI<=NX-1;AI++){
HX[AI][J]=HX[AI][J]-CHXLY[AI][J]*(EZ[AI][J+1]-EZ[AI][J]);
}
}
//HYの計算
for(J=2;J<=NY-1;J++){
for(AI=1;AI<=NX-1;AI++){
HY[AI][J]=HY[AI][J]+CHYLX[AI][J]*(EZ[AI+1][J]-EZ[AI][J]);
}
}
//--- 2-5. 電磁場の磁場成分HX,HYの計算(「 CALL H-field 」) ここまで-------
// //--- 2-6. 過去の値の更新 --------------------------------------------------
//
// for(AI=1;AI<=NX;AI++){
// for(J=1;J<=NY;J++){
// EZ1[AI][J]=EZ[AI][J];
// HX1[AI][J]=HX[AI][J];
// HY1[AI][J]=HY[AI][J];
// }
// }
//
// //--- 2-6. 過去の値の更新 ここまで------------------------------------------
cout << setw(10) << "IT=" << IT <<endl; // ターン数の画面表示
//--- 2-7. 振幅データの作成(本にはない)--------------------------------------
if((ITMAX-IT) <= 100){
for(J=1;J<=NY;J++){
for(AI=1;AI<=NX;AI++){
// ラスト100ターンでの最大値を取る
if(abs(EZ[AI][J]) >= AMPEZ[AI][J]){
AMPEZ[AI][J]=abs(EZ[AI][J]); // AMPEZ[AI][J]=EZの振幅
}
}
}
}
//--- 2-7. 振幅データの作成 ここまで----------------------------------------
//--- 2-8. 半時間ステップだけ時間を進める------------------------------------
T=T+DT/2;
//--- 2-8. 半時間ステップだけ時間を進める ここまで---------------------------
//--- 2-9. 1000ステップ毎終了時のEZ[AI][J]を 各ファイルに出力 ------
if(IT==10){
kakikomi("10ターン目をファイルに出力します。" ,
"20100729_01_10_DATA.txt" , EZ , NX , NY);
}
if(IT==8000){
kakikomi("8000ターン目をファイルに出力します。" ,
"20100729_01_8000_DATA.txt" , EZ , NX , NY);
}
if(IT==10000){
kakikomi("10000ターン目をファイルに出力します。" ,
"20100729_01_10000_DATA.txt" , EZ , NX , NY);
}
//--- 2-9. 100、300、500ステップ終了時のEZ[AI][J]を 各ファイルに出力
// ここまで --------------------------------------
//---2-10. 1ポイントモニター ----------------------------------------------
//１点でのEZ（AI,J）の値をすべての時間ステップで記録する
//適切な全体のステップ数を求めるのが目的
//記録(8,8)でのEZをモニターする
// EZ_MON_1[IT]=EZ[8][8];
//---2-10. 1ポイントモニター ここまで-------------------------------------
//--- 2-11. メインループの終わり-------------------------------------------
}
//--- 2-11. メインループの終わり ここまで----------------------------------
//=== 2. メインループ ここまで ================================================
//=== 3. ファイル出力(本にはない) ==============================================
//--- 3-1. AMPEZ[AI][J]の出力 ---------------
kakikomi("AMPEZ[AI][J]を出力します。" ,
"20100729_01_01_OUTPUT_AMPEZ_DATA.txt" , AMPEZ , NX , NY);
//--- 3-2. 最後のEZ[AI][J]の出力 ------------
kakikomi("EZ[AI][J]を出力します。" ,
"20100729_01_OUTPUT_LAST_EZ_DATA.txt" , EZ , NX , NY);
//--- 3-3. 最後のHX[AI][J]の出力 ------------
kakikomi("HX[AI][J]を出力します。" ,
"20100729_01_OUTPUT_LAST_HX_DATA.txt" , HX , NX , NY);
//--- 3-4. 最後のHY[AI][J]の出力 ------------
kakikomi("HY[AI][J]を出力します。" ,
"20100729_01_OUTPUT_LAST_HY_DATA.txt" , HY , NX , NY);
//--- 3-5. EZ_MON_1の出力 -------------------
onePointMonitor("EZ_MON_1を出力します。" ,
"20100729_01_OUTPUT_EZ_MON_1_DATA.txt" , EZ_MON_1 , ITMAX);
//=== 3. ファイル出力(本にはない) ここまで ======================
//=== 4. フリーストアの配列の削除 ================================
deleteArray(CEZ,NX+1); //配列名以外に行数だけ渡す
deleteArray(CEZLX,NX+1);
deleteArray(CEZLY,NX+1);
deleteArray(CHXLY,NX+1);
deleteArray(CHYLX,NX+1);
deleteArray(MEDIUM,NX+1);
deleteArray(C,NX+1);
deleteArray(DA,NX+1);
deleteArray(CV,NX+1);
deleteArrayInt(KNOWN,NX+1);
deleteArray(HX,NX+1);
// deleteArray(HX1,NX+1);
deleteArray(HY,NX+1);
// deleteArray(HY1,NX+1);
deleteArray(EZ,NX+1);
// deleteArray(EZ1,NX+1);
delete[] EZ_MON_1;
deleteArray(EZX1,5);
deleteArray(EZX2,5);
deleteArray(EZY1,NX+1);
deleteArray(EZY2,NX+1);
deleteArray(AMPEZ,NX+1);
//===4. フリーストアの配列の削除 ここまで =======================
cout << "何かキーを押してください";
getchar();
}
// 全プログラム終了