! -----------------------------------------------|
! FDTD2d-mur2.f90 |
! |
! FDTD法による直交座標系(x,y)における |
! パラボラミラーによる電磁波の反射の数値計算 |
! |
! |
! プログラム作成：駒田寛 2010年4月 |
! -----------------------------------------------|
!=== 1. setup ==================================================================
! 宇野先生の本で「 setup 」となっている部分とほぼ同じである
!---1. setup 目次-------------------------------------------------------
! 1-0. 配列宣言
! 1-1．空間、時間に対して定数であるものの設定
! 1-2．空間に対して変数、時間に対して定数であるものの設定
! 1-3．時間に対して変数であるものの初期値設定
! 1-4. メモデータ作成（本にはない）
!---1. setup 目次 ここまで---------------------------------------------
!--- 1-0.配列宣言------------------------------------------
implicit none
REAL CEZ(1001,1001),CEZLX(1001,1001),CEZLY(1001,1001)
REAL CHXLY(1001,1001),CHYLX(1001,1001)
REAL MEDIUM(1001,1001)
REAL C(1001,1001) !導電率
REAL EZDUM(1001,1001) !仮座標でのEZ
REAL DA(1001,1001)
INTEGER NX
INTEGER NY
REAL PAI
REAL E
REAL MU
REAL F
REAL W
REAL DT
REAL DX,DY
REAL THETA
REAL THETA1,THETA2,THETA3,THETA4
REAL CV0 !真空中の光速(シーブイゼロ)
REAL CV(1001,1001) !媒質中の光速
INTEGER K !時間ステップ数格納用
INTEGER J !Y方向ステップ数格納用
INTEGER AI !X方向ステップ数格納用
INTEGER JDUM !仮座標用変数
INTEGER AIDUM !仮座標用変数
REAL HX(1001,1001),HY(1001,1001)
REAL EZ(1001,1001)
REAL EZ_MON_1(10000)
REAL EZ_MON_2(10000)
REAL EZX1(4,1001),EZX2(4,1001)
REAL EZY1(1001,4),EZY2(1001,4)
INTEGER KNOWN(1001,1001)
INTEGER MIRROR(1001)
INTEGER HORN_UP(1001)
INTEGER HORN_DOWN(1001)
INTEGER IT
INTEGER ITMAX
REAL T
REAL C1,C2,C3
!本のCXD,CXX,CXFYDをそれぞれC1,C2,C3と置いた。
!２次元ではCXFZDはない。
REAL AMPEZ(1001,1001)
!--- 1-0.配列宣言 ここまで-----------------------------------------------
!--- 1-1．空間、時間に対して定数（物理定数も含む）であるものの設定----------
!変数 ITMAX 以外。
! x,y方向の分割数+1の値
NX=1001 ! 吸収境界の座標はX,Yともに1,1001とする
NY=1001
!円周率
PAI = 3.141592654e0
!真空中の誘電率 [F/m]
E = 8.8541878e-12
!真空中の透磁率 [H/m]
MU = 4 * PAI * 1.0e-7
!調べる電磁波の周波数［Hz］
F = 15.0e+9
!電磁波の角周波数 [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およびメモ帳にある
!座標回転用角度[degree]（実際 マイナス ダミー）
THETA = 225.0e0
THETA1 = cos(THETA*PAI/180)
THETA2 = -sin(THETA*PAI/180)
THETA3 = sin(THETA*PAI/180)
THETA4 = cos(THETA*PAI/180)
CV0 = 1/(E*MU)**0.5e0 !真空中の光速
C1 = (CV0*DT-DX)/(CV0*DT+DX)
C2 = 2*DX/(CV0*DT+DX)
C3 = DX*( (CV0*DT)**2 ) / ( 2* ((DY)**2) *(CV0*DT+DX) )
!変数は用意してあるが本プログラムでは常にDX=DYとする。
ITMAX = 5000; !全ステップ回数
!画面での確認、ファイルへの出力
OPEN(1101,FILE='20100421_01_2ramda_SET_DATA_CONST.txt',STATUS='unknown')
WRITE(6,'(a)')'以下をSET_DATA_1ファイルに出力します'
WRITE(6,'(a30,e10.5)')'X方向のセルサイズ [m]=',DX
WRITE(1101,'(a30,e10.5)')'X方向のセルサイズ [m]=',DX
WRITE(6,'(a30,e10.5)')'Y方向のセルサイズ [m]=',DY
WRITE(1101,'(a30,e10.5)')'X方向のセルサイズ [m]=',DX
WRITE(6,'(a30,e10.5)')'時間ステップ［s］=',DT
WRITE(1101,'(a30,e10.5)')'時間ステップ［s］=',DT
WRITE(6,'(a30,e10.5)')'１秒当たりのステップ回数 =',1/DT
WRITE(1101,'(a30,e10.5)')'１秒当たりのステップ回数 =',1/DT
WRITE(6,'(a30,e10.5)')'真空中の誘電率 [F/m]=',E
WRITE(1101,'(a30,e10.5)')'真空中の誘電率 [F/m]=',E
WRITE(6,'(a30,e10.5)')'調べる電磁波の周波数［Hz］=',F
WRITE(1101,'(a30,e10.5)')'調べる電磁波の周波数［Hz］=',F
WRITE(6,'(a30,e10.5)')'１周期当たりのステップ回数 =',1/DT/F
WRITE(1101,'(a30,e10.5)')'１周期当たりのステップ回数 =',1/DT/F
WRITE(6,'(a30,e10.5)')'真空中の透磁率 [H/m]=',MU
WRITE(1101,'(a30,e10.5)')'真空中の透磁率 [H/m]=',MU
WRITE(6,'(a30,i5)')'X方向の分割数＋１の値 =',NX
WRITE(1101,'(a30,i5)')'X方向の分割数＋１の値 =',NX
WRITE(6,'(a30,i5)')'Y方向の分割数＋１の値 =',NY
WRITE(1101,'(a30,i5)')'Y方向の分割数＋１の値 =',NY
WRITE(6,'(a30,e10.5)')'円周率 =',PAI
WRITE(1101,'(a30,e10.5)')'円周率 =',PAI
WRITE(6,'(a30,e10.5)')'電磁波の角周波数 [rad/s] =',W
WRITE(1101,'(a30,e10.5)')'電磁波の角周波数 [rad/s] =',W
WRITE(6,'(a30,e10.5)')'真空中の光速[m/s] =',CV0
WRITE(1101,'(a30,e10.5)')'真空中の光速[m/s] =',CV0
WRITE(6,'(a30,e10.5)')'全ステップ回数 =',ITMAX
WRITE(1101,'(a30,e10.5)')'全ステップ回数 =',ITMAX
CLOSE(1101)
WRITE(6,'(a)')'何かキーを押してください'
pause
!--- 1-1．空間、時間に対して定数であるものの設定 ここまで-------------------
!--- 1-2．空間的には変数、時間的には定数であるものの設定--------------------
! 1-1に該当するものもあるが、揺れ動くものはすべてこちらにした。
! KNOWN(AI,J)を除く。
!--- 1-2-1. 仮設定 -----------------------------------------------------
DO 6010 J=1 , NY
DO 6010 AI=1 , NX
CEZLX(AI,J)=0.0e0 !EZを求める式に使用する係数
CEZLY(AI,J)=0.0e0 !EZを求める式に使用する係数
CHXLY(AI,J)=0.0e0 !HXを求める式に使用する係数
CHYLX(AI,J)=0.0e0 !HYを求める式に使用する係数
CEZ(AI,J)=0.0e0 !EZを求める式に使用する係数
C(AI,J)=0.0e0 !導電率
!金属部分については電界＝0として別の場所で表現する。
!電界を隣の値から求めるときにCを使うが
!金属部分についてはそれをやらないので値は関係ない
!(0のままでも良い)。
EZDUM(AI,J)=1.0e0 !仮座標でのEZ
6010 CONTINUE
!ホーンの形を決める配列HORN_UP(AI)の初期化
DO AI=1,NX
HORN_UP(AI)=0
ENDDO
!ホーンの形を決める配列HORN_DOWN(AI)の初期化
DO AI=1,NX
HORN_DOWN(AI)=0
ENDDO
!ホーンの形を決める配列HORN_UP(AI)を計算で決定する
DO AI=451,551
HORN_UP(AI)=INT(0.2*AI+410.8)
ENDDO
!ホーンの形を決める配列HORN_DOWN(AI)を計算で決定する
DO AI=451,551
HORN_DOWN(AI)=INT(-0.2*AI+591.2)
ENDDO
!仮座標でのホーン部EZの設定（設定領域は 251<=AI,J<=751 とする）
DO AI=451,551
EZDUM(AI,HORN_UP(AI))=0.0e0
EZDUM(AI,HORN_DOWN(AI))=0.0e0
ENDDO
!--- 1-2-1. 仮設定 ここまで--------------------------------------------
!--- 1-2-2. 本設定 ----------------------------------------------------
!その他のパラメータの設定
DO 6020 J=1, NY
DO 6020 AI=1, NX
MEDIUM(AI,J)=1.0e0 !比誘電率（今回は誘電体はない）
CV(AI,J) = 1/(MEDIUM(AI,J)*E*MU)**0.5e0 !媒質中の光速 [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 !同上
6020 CONTINUE
!ミラーの形を決める配列MIRROR(AI)の初期化
DO AI=1,NX
MIRROR(AI)=0
ENDDO
!ミラーの形を決める配列MIRROR(AI)を計算で決定する
DO AI=201,401
MIRROR(AI)=INT( ((AI-501)**2)/800+301 )
ENDDO
!--- 1-2-2. 本設定 ここまで ------------------------------------------
!--- 1-2．空間的には変数、時間的には定数であるものの設定 ここまで-----------
!--- 1-3．時間に対して変数であるものの初期値設定----------------------------
!--- 1-3-1.電磁場の初期値の設定 ----------------------------------------
!現時点での値
DO J=1 , NY
DO AI=1 , NX
HX(AI,J)=0.0e0 !磁場のＸ成分（0とする、以下同様）
HY(AI,J)=0.0e0 !磁場のＹ成分
EZ(AI,J)=0.0e0 !電場のＺ成分
ENDDO
ENDDO
!過去の値１
!AI=1,NXの辺とその隣の辺（計4本の線）の1ターン、２ターン前のEZ
DO AI=1 , 4
DO J=1 , NY
EZX1(AI,J)=0.0e0
EZX2(AI,J)=0.0e0
ENDDO
ENDDO
!過去の値２
!J=1,NYの辺とその隣の辺（計4本の線）の1ターン、２ターン前のEZ
DO J=1 , 4
DO AI=1 , NX
EZY1(AI,J)=0.0e0
EZY2(AI,J)=0.0e0
ENDDO
ENDDO
!--- 1-3-1.電磁場の初期値の設定 ここまで--------------------------------
!--- 1-3-2.１ポイントモニタ用配列 EZ_MON の初期化-----------------------
DO K=1,10000
EZ_MON_1(K)=0.0e-0
EZ_MON_2(K)=0.0e-0
ENDDO
!--- 1-3-2.１ポイントモニタ用配列 EZ_MON の初期化 ここまで-------------
!--- 1-3-3. 経過時間 T の初期化 ----------------------------------------
T=0
!--- 1-3-3. 経過時間 T の初期化 ここまで-------------------------------
!--- 1-3．時間に対して変数であるものの初期値設定 ここまで------------------
!--- 1-4. メモデータを作る(本にはない)--------------------------------------
OPEN(1421,FILE='20100421_01_2ramda_SET_DATA_C.txt',STATUS='unknown')
WRITE(6,'(a)')'導電率C（AI,J)ファイルを作成します '
! WRITE(1421,'(a)')'導電率C（AI,J) ' matlabの関係で書かない
DO AI=1,NX,1
WRITE(1421,501) (C(AI,J),J=1,NY,1)
501 FORMAT( e10.2, 1000(2X, e10.2) )
ENDDO
CLOSE(1421)
WRITE(6,'(a)')'何かキーを押してください'
pause
!--- 1-4. メモデータを作る ここまで--------------------------------------
!=== 1. setup ここまで ======================================================
!=== 2. メインループ ========================================================
!--- 2-1. メインループの始まり--------------------------------------------
DO 6130 IT=1,ITMAX ! 計算のスタート
! 経過時間＝(IT)＊DT
!--- 2-1. メインループの始まり ここまで------------------------------------
!--- 2-2. 電磁場のEZの計算(本では 「 CALL E-field 」)---------------------
!(P15にある)
!AI,Jともに内側を全てスキャン
!---2-2-1. 本体の計算 --------------------------------------------------
DO J=2,NY-1
DO AI=2,NX-1
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))
ENDDO
ENDDO
!---2-2-1. 本体の計算 ここまで------------------------------------------
!--- 2-2-2. 境界（内部）の設定(本には明記されていない)------------------
!ホーン部(ダミー図形を(501,501)を中心としてTHETA度回転した図形を描く)
DO AIDUM=451,551
AI = INT( (AIDUM-501)*THETA1 +(HORN_UP(AIDUM)-501)*THETA2 +501 )
J = INT( (AIDUM-501)*THETA3 +(HORN_UP(AIDUM)-501)*THETA4 +501 )
EZ(AI,J)=0.0e0
AI = INT( (AIDUM-501)*THETA1 +(HORN_UP(AIDUM)+1-501)*THETA2 +501 )
J = INT( (AIDUM-501)*THETA3 +(HORN_UP(AIDUM)+1-501)*THETA4 +501 )
EZ(AI,J)=0.0e0
AI = INT( (AIDUM-501)*THETA1 +(HORN_DOWN(AIDUM)-501)*THETA2 +501 )
J = INT( (AIDUM-501)*THETA3 +(HORN_DOWN(AIDUM)-501)*THETA4 +501 )
EZ(AI,J)=0.0e0
AI = INT( (AIDUM-501)*THETA1 +(HORN_DOWN(AIDUM)-1-501)*THETA2 +501 )
J = INT( (AIDUM-501)*THETA3 +(HORN_DOWN(AIDUM)-1-501)*THETA4 +501 )
EZ(AI,J)=0.0e0
ENDDO
!波源（(461,501)を(501,501)を中心としてTHETA度回転した点にとる）
AI = INT( (461-501)*THETA1 +(501-501)*THETA2 +501 )
J = INT( (461-501)*THETA3 +(501-501)*THETA4 +501 )
EZ(AI,J)=1.0e0*sin(W*T) ! 各ターンでSIN波波源を設定している
!パラボラ
! DO AI=201,401
! EZ(AI,MIRROR(AI))=0.0e0
! ! 電場が0となるように各ターンで修正する
! ENDDO
!--- 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))
DO J=3,NY-2
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))
ENDDO
!過去の値の更新
DO J=2,NY-1
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)
ENDDO
!---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))
DO AI=3,NX-2
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))
ENDDO
!過去の値の更新
DO AI=2,NX-1
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)
ENDDO
!--- 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 」)----------------
!--- 2-5-1. 境界（内部）の設定（本には明示していない）------------------
! HX(AI,300)=-EZ(AI,300)/120/PAI
!真空中の特性インピ-ダンスを用いて磁場も決定している(?)
!波源部分の磁界
! 進行方向のわかった波では付けると精度は上がるはずだが
! 点波源ではいらないはずである
!--- 2-5-1. 境界（内部）の設定 ここまで--------------------------------
!--- 2-5-2. 本体の計算 -------------------------------------------------
!HXの計算（P15)
DO J=1,NY-1
DO AI=2,NX-1
HX(AI,J)=HX(AI,J)-CHXLY(AI,J)*(EZ(AI,J+1)-EZ(AI,J))
ENDDO
ENDDO
!HYの計算
DO J=2,NY-1
DO AI=1,NX-1
HY(AI,J)=HY(AI,J)+CHYLX(AI,J)*(EZ(AI+1,J)-EZ(AI,J))
ENDDO
ENDDO
!--- 2-5-2. 本体の計算 ここまで----------------------------------------
!--- 2-5. 電磁場の磁場成分HX,HYの計算(「 CALL H-field 」) ここまで-------
WRITE(6,*)'IT=',IT ! ターン数の画面表示
!--- 2-6. 振幅データの作成(本にはない)--------------------------------------
IF((ITMAX-IT).LT.100) THEN
DO 6120 J=1,NY
DO 6120 AI=1,NX
! ラスト100ターンでの最大値を取る
IF(ABS(EZ(AI,J)).GE.AMPEZ(AI,J)) THEN
AMPEZ(AI,J)=ABS(EZ(AI,J)) ! AMPEZ(AI,J)=EZの振幅
ENDIF
6120 CONTINUE
ENDIF
!--- 2-6. 振幅データの作成 ここまで----------------------------------------
!--- 2-7. 半時間ステップだけ時間を進める------------------------------------
T=T+DT/2
!--- 2-7. 半時間ステップだけ時間を進める ここまで---------------------------
!--- 2-8. 各ステップ終了時のEZ(AI,J)を 1000 3000 ファイルに出力 ----------
!オプション
!計算結果出力ファイル 1000 3000 を開く
!（各ステップの結果を全ポイント記録するために開いておく）
!（３次元配列は事実上ムリ）
IF(IT==1000)then
OPEN(14,FILE='20100421_01_2ramda_1000_DATA.txt',STATUS='unknown')
! WRITE(14,'(a,2x,i4)')'IT=',IT matlabの関係で書かない
DO AI=1,NX,1 !間引きせずに書き込む
WRITE(14,504) (EZ(AI,J),J=1,NY,1)
504 FORMAT( e10.2, 1000(2X, e10.2) )
ENDDO
CLOSE(14)
ENDIF
IF(IT==3000)then
OPEN(15,FILE='20100421_01_2ramda_3000_DATA.txt',STATUS='unknown')
! WRITE(15,'(a,2x,i4)')'IT=',IT
DO AI=1,NX,1 !間引きせずに書き込む
WRITE(15,504) (EZ(AI,J),J=1,NY,1)
ENDDO
CLOSE(15)
ENDIF
!--- 2-8. 各ステップ終了時のEZ(AI,J)を 1000 3000 ファイルに出力
! ここまで --------------------------------------
!---2-9. 1ポイントモニター ------------------------------------------------
!１点でのEZ（AI,J）の値をすべての時間ステップで記録する
!適切な全体のステップ数を求めるのが目的
!記録(800,500)でのEZをモニターする
EZ_MON_1(IT)=EZ(800,500)
! EZ_MON_2(IT)=EZ(499,499)
!---2-9. 1ポイントモニター ここまで---------------------------------------
!--- 2-10. メインループの終わり-------------------------------------------
6130 CONTINUE
!--- 2-10. メインループの終わり ここまで----------------------------------
!=== 2. メインループ ここまで ================================================
!=== 3. ファイル出力(本にはない) ==============================================
!--- 3-1. AMPEZ(AI,J)の出力 -----------------------------------------------
!計算結果出力ファイル OUTPUT_AMPEZ_DATA を開く
OPEN(16,FILE='20100421_01_2ramda_OUTPUT_AMPEZ_DATA.txt',STATUS='unknown')
!算出されたEZ振幅マップAMPEZ(AI,J)を出力する
WRITE(6,'(a)')'AMPEZ(AI,J)を出力します'
! WRITE(16,'(a)')'AMPEZ(AI,J)' matlab処理の関係で書かない
DO AI=1,NX,1 !間引きせずに書き込む
WRITE(16,505) (AMPEZ(AI,J),J=1,NY,1)
505 FORMAT( e10.2, 1000(2X, e10.2) )
ENDDO
CLOSE(16) ! OUTPUT_AMPEZ_DATA ファイルを閉じる
!--- 3-1. AMPEZ(AI,J)の出力 ここまで---------------------------------------
!--- 3-2. 最後のEZ(AI,J)の出力 -----------------------------------------------
!計算結果出力ファイル OUTPUT_AMPEZ_DATA を開く
OPEN(17,FILE='20100421_01_2ramda_OUTPUT_LAST_EZ_DATA.txt',STATUS='unknown')
!算出されたEZ振幅マップAMPEZ(AI,J)を出力する
WRITE(6,'(a)')'最後のEZ(AI,J)を出力します'
! WRITE(16,'(a)')'AMPEZ(AI,J)' matlab処理の関係で書かない
DO AI=1,NX,1 !間引きせずに書き込む
WRITE(17,508) (EZ(AI,J),J=1,NY,1)
508 FORMAT( e10.2, 1000(2X, e10.2) )
ENDDO
CLOSE(17) ! LAST_EZ_DATA ファイルを閉じる
!--- 3-2. 最後のEZ(AI,J)の出力 ここまで---------------------------------------
!--- 3-3. EZ_MON_1(10000)の出力 ---------------------------------------------
!計算結果出力ファイル OUTPUT_EZ_MON_1_DATA を開く
OPEN(18,FILE='20100421_01_2ramda_OUTPUT_EZ_MON_1_DATA.txt',STATUS='unknown')
!EZ_MON(10000)を出力する
WRITE(6,'(a)')'EZ_MON_1(10000)を出力します'
! WRITE(18,'(a)')'EZ_MON_1(10000)' matlab処理の関係で書かない
DO K=1,ITMAX
WRITE(18,'(i4)') K
WRITE(18,506) EZ_MON_1(K)
506 FORMAT( e10.2 )
ENDDO
CLOSE(18) ! OUTPUT_EZ_MON_1_DATA ファイルを閉じる
!--- 3-3. EZ_MON_1(10000)の出力 ここまで-------------------------------------
!=== 3. AMPEZ(AI,J)のファイル出力(本にはない) ここまで ======================
Pause
STOP
END
! 全プログラム終了