くだすれFORTRAN(超初心者用) at TECH
[2ch|▼Menu]
1:デフォルトの名無しさん
06/01/24 09:48:23
このスレッドは、他のスレッドでは書き込めない超低レベル、 
もしくは質問者自身何が何だが分からない質問を勇気を持って書き込むスレッドです。 
FORTRAN使いが優しくコメントを返しますが、 
お礼はFORTRANの布教と初心者の救済をお願いします。 


2:漂泊の2ゲッター
06/01/24 17:19:15
       ローリング!!    ∧∧  
               (゚Д゚,,)  
               ⊂⊂,,ヽ  
                (_ (  )ノ
     クルン       
              /⌒⌒ヽノ  )))
              (   )て )  
          (((   ∨∨⊂ノ
ズサギコ!! 
                    (´´
   ∧∧  )  ≡≡≡≡≡(´⌒(´≡≡ 
 ⊂(゚Д゚⊂⌒つ ≡≡≡≡(´⌒;;≡ 
           (´⌒(´⌒;;
2げっとーーーーーーーーー!!!

3:毛の生えたフサギコ
06/01/25 20:48:31
なにがなんだが ×
なにがなんだか ○

すみません、私が悪かったです。コピペするとき
添削してくれたらよかったのですが・・・。

4:お願いします;;
06/01/27 01:28:03
任意の生年月日と今日の日付を入力すれば、
産まれてから今日まで生きてきた日数を出力するプログラムを作れ。
さらに、任意の整数nを入力して、産まれてからn日目が、
いつかを出力させるプログラムを作れ。
ちなみに、閏年は以下のようにして決められる。

1)基本的に西暦が4で割り切れたら閏年。
2)しかし、西暦が100で割り切れたら閏年ではない。
3)しかし、西暦が400で割り切れたら閏年とする。

レポート作成
予備的考察、プログラム、プログラムの説明、計算結果、考察

フォートラン使用・・・


一応書き込んでみました。

少しあきらめられない自分がいました・・・。

5:デフォルトの名無しさん
06/01/27 01:41:50
>>4
Rosenmaidenを見ていたおれ様が来ましたよw
まあ少し待たれよw

6:デフォルトの名無しさん
06/01/27 02:39:51
MODULE mod_subs
IMPLICIT NONE
INTEGER, PARAMETER :: ndays_of_month(12, 2) = &
RESHAPE( (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, &
31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /), (/12, 2/) )
CONTAINS
LOGICAL FUNCTION qleapyear(iyear)
IMPLICIT NONE
INTEGER, INTENT(IN) :: iyear
qleapyear = .FALSE.
IF ( MOD(iyear, 4) == 0 ) qleapyear = .TRUE.
IF ( MOD(iyear, 100) == 0 ) qleapyear = .FALSE.
IF ( MOD(iyear, 1000) == 0 ) qleapyear = .TRUE.
RETURN
END FUNCTION qleapyear
!
INTEGER FUNCTION idays_from_jan1(iyear, imonth, iday)
IMPLICIT NONE
INTEGER, INTENT(IN) :: iyear, imonth, iday
INTEGER :: i, ileap, n
ileap = 1
IF ( qleapyear(iyear) ) ileap = 2
n = 0
DO i = 1, imonth - 1
n = n + ndays_of_month(i, ileap)
END DO
idays_from_jan1 = n + iday
RETURN
END FUNCTION idays_from_jan1
END MODULE mod_subs



7:デフォルトの名無しさん
06/01/27 02:41:02
PROGRAM AliceGame
USE mod_subs
IMPLICIT NONE
INTEGER :: iyear0, imonth0, iday0, iyear1, imonth1, iday1
INTEGER :: iyear, ndays, ndays0, ndays1
!
WRITE(*, *) "input birth date : Year, Month, Day"
READ(*, *) iyear0, imonth0, iday0
!
WRITE(*, *) "input today's date : Year, Month, Day"
READ(*, *) iyear1, imonth1, iday1
!
ndays0 = idays_from_jan1(iyear0, imonth0, iday0)
ndays1 = idays_from_jan1(iyear1, imonth1, iday1)
ndays = -ndays0
DO iyear = iyear0, iyear1 - 1
IF ( qleapyear(iyear) ) THEN
ndays = ndays + 366
ELSE
ndays = ndays + 365
END IF
END DO
ndays = ndays + ndays1
WRITE(*, *) "The number of days from birth =", ndays
STOP
END PROGRAM AliceGame

専ブラ インデント付き用
>>6-7

8:デフォルトの名無しさん
06/01/27 02:48:07
>>4
結果を余りチェックしていないので、間違ってたらごめんwwww謝罪も賠償もしないww

とりあえず、課題の前半だ。後半は眠いから作らないw

紀元前の日付には対応していないが、入力値はチェックしていないので危ない。

出力例:
input birth date : Year, Month, Day
1900 1 1
input today's date : Year, Month, Day
2000 1 1
The number of days from birth = 36524
Press any key to continue

1世紀で1年の平均日数が365.24日になっているのでグレゴリオ暦としては正しいと
思うが、間違ってたらごめんw

9:デフォルトの名無しさん
06/01/27 02:54:18
>>4
日数計算のアルゴリズムはこんな感じ↓
5/10の10日先は、10+10=20なので、5/20。
5/10の20日先は、10+20=30なので、5/30。
5/10の30日先は、10+30=40なので、5/40、
しかし、5月は31日までなので、日は40-31=9、月は5+1=6。よって、6/9。
これを繰り返せば、任意の日付のx日後を求められる。
3/23 + 100 = 3/123 = 4/92 = 5/62 = 6/31 = 7/1 といった要領で。

閏年の判定は
LOGICAL FUNCTION CheckLeap(year)
とかいう関数を作って、
入力:西暦(整数型year)、出力:フラグ(trueなら閏年)
みたいにしておけばいいと思う。

やる気があれば、これくらいのヒントで十分作れるだろう。
あとは頑張れ。


10:デフォルトの名無しさん
06/01/27 02:54:58
プログラムの解説だ。
関数 qleapyear はうるう年なら .TRUE. を、普通の年なら .FALSE. を返す。
関数 idays_from_jan1 は、1月1日から入力日までの日数を返す。
  ( 1月1日を入れた場合は0日、1月2日を入れたら1日と数える。 )
主プログラムのアルゴリズム
最初の年は、365(または366)−誕生日までの日数。
最後の年は、今日までの日数。
途中の年は、うるう年を考慮しつつ365(366)を足してゆく。

料金の100万円はマカオ銀行に金正日あてで振り込んでおいてくれ。

11:9
06/01/27 02:56:16
ありゃあ作っちゃってるよ。
投稿するの遅かったな。
でも、module文とか使っていいの?
Fortran77じゃ動かないけど…。


12:5
06/01/27 03:06:30
>>11
ごめw 
練習のため90で書いてしまったw 
ともかく後半を任せた もう寝るw



13:デフォルトの名無しさん
06/01/27 09:53:36
次の宿題もってこい!w

14:助けて下さい(:_;)
06/01/27 20:05:54
誤差関数E(x)を次のように定義する。
E(x)=2/√π∫(範囲:0からx)e(-yの2乗)乗dy
この時、E(x)=2x(2)乗-x(3)乗
の実数解を「すべて」求めよ。(解はひとつとは限らないので解の個数についても議論せよ。)

このプログラムがどうしてもできません…fortranで作ってほしいですm(--)m
if文とループがうまく使えなくてエラーばっかりでてしまいます↓時間がないので助けてほしいですM(--)M×∞


15:5
06/01/27 21:58:04
>>4 後半。だるいので適当w  マイナス入れると暴走、入力チェックしてないw

MODULE mod_subs
IMPLICIT NONE
INTEGER, PARAMETER :: ndays_of_month(12, 2) = &
RESHAPE( (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, &
31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /), (/12, 2/) )
CONTAINS
LOGICAL FUNCTION qleapyear(iyear)
IMPLICIT NONE
INTEGER, INTENT(IN) :: iyear
qleapyear = .FALSE.
IF ( MOD(iyear, 4) == 0 ) qleapyear = .TRUE.
IF ( MOD(iyear, 100) == 0 ) qleapyear = .FALSE.
IF ( MOD(iyear, 1000) == 0 ) qleapyear = .TRUE.
RETURN
END FUNCTION qleapyear
END MODULE mod_subs


16:5
06/01/27 21:59:02
PROGRAM AliceGame2
USE mod_subs
IMPLICIT NONE
INTEGER :: iyear0, imonth0, iday0, iyear, imonth, iday
INTEGER :: ileap, ndays, idays
!
WRITE(*, *) "input birth date : Year, Month, Day"
READ(*, *) iyear0, imonth0, iday0
!
WRITE(*, *) "input days "
READ(*, *) ndays
!
ndays = ndays + iday0
iyear = iyear0
DO
IF ( qleapyear(iyear) ) THEN
ileap = 2
idays = 366
ELSE
ileap = 1
idays = 365
END IF
IF ( ndays - idays > 0 ) THEN
iyear = iyear + 1
ndays = ndays - idays
ELSE
EXIT
END IF
END DO


17:5
06/01/27 22:00:19
!
imonth = imonth0
DO
IF ( ndays - ndays_of_month(imonth, ileap) > 0 ) THEN
imonth = imonth + 1
ndays = ndays - ndays_of_month(imonth, ileap)
ELSE
EXIT
END IF
END DO
iday = ndays
WRITE(*, *) "The date is ", iyear, imonth, iday
STOP
END PROGRAM AliceGame2

センブラインデント用
>>15-17

結果チェックしてないw 自分でチェックして書き直すですよ!

18:5
06/01/27 22:05:32
ごめw バグってたwwww

19:5
06/01/27 22:08:35
IF ( ndays - ndays_of_month(imonth, ileap) > 0 ) THEN
 imonth = imonth + 1
 ndays = ndays - ndays_of_month(imonth, ileap)
ELSE

IF ( ndays - ndays_of_month(imonth, ileap) > 0 ) THEN
 ndays = ndays - ndays_of_month(imonth, ileap)
 imonth = imonth + 1
ELSE

順番入れ替えてくれw

20:デフォルトの名無しさん
06/01/27 23:32:37
>>14
f(x)=ERF(x) - 2x**2 + x**3 
とおく。求めたいxはf(x)=0を満たすx。

1階及び2階の微係数を求めると、
f'(x)=exp(-x**2) - 4x + 3x**2
f"(x)=-2x exp(-2x**2) - 4 + 6x

また、それぞれのx=0とx→∞の値は、
f(0) = 0, lim f(x) >0 
f'(0) > 0, lim f'(x) >0 
f"(0) < 0, lim f"(x) >0

ここでf"(x)の概形を調べてみると
f"(x) = 2x( 3 - exp(-2x**2) ) - 4 > 2x(3 - 1) - 4 = 4x - 4 〜 1次直線

これらのグラフを書けば、f'(x)はx〜1付近で最小値を持つ。
f'(1) < 0 なので、f'(0) > 0, lim f'(x) >0より f'(x) は2回0点を切ることがわかる。
つまり、f(x)は2個極値を持つ。

さらにこのことと、f(0) = 0, lim f(x) >0 また f(1) < 1 - 2 + 1 = 0 より fはx>0で
2回0点を切ることが分かる。

ゆえに自明な解x=0も含めて3個解を持つ。
(その解はx=1付近に前後に1こづづあるはずww)


21:デフォルトの名無しさん
06/01/27 23:35:42
FORTRANの宿題っていうか授業してるところってまだあるんだね。
機械とか建築関係ならまだ多いのかな?

22:デフォルトの名無しさん
06/01/27 23:41:14
理学部はFORTRANしかやらねーな。もしくはFORTRANすらやらねーな。

23:デフォルトの名無しさん
06/01/28 00:22:53
>>14
とりあえず、誤差関数ERF(x)は持っていると仮定した。
表式が解析的に表されているので、微係数の式を求めた上で、ニュートン・ラフソン法でといた。

MODULE mod_sub
IMPLICIT NONE
CONTAINS
REAL FUNCTION f(x)
IMPLICIT NONE
REAL, INTENT(IN) :: x
REAL :: ERF
EXTERNAL :: ERF
f = ERF(x) - 2.0 * x**2 + x**3
END FUNCTION f
!
REAL FUNCTION fderiv(x)
IMPLICIT NONE
REAL, INTENT(IN) :: x
REAL :: pi
pi = 4.0 * ATAN(1.0)
fderiv = 2.0 / SQRT(pi) * EXP(-x**2) - 4.0 * x + 3.0 * x**2
END FUNCTION fderiv
END MODULE mod_sub
!


24:デフォルトの名無しさん
06/01/28 00:23:38
PROGRAM erf0
USE mod_sub
IMPLICIT NONE
INTEGER :: i
REAL :: x
WRITE(*, *) ' input start value x0 :'
READ(*, *) x
DO i = 1, 100
x = x - f(x) / fderiv(x) ! Newton-Raphson method
PRINT *, x, f(x) ! debug info
IF ( ABS(f(x)) < EPSILON(x) ) THEN
WRITE(*, *) ' Root found x= ', x, ' : f(x) =', f(x)
STOP
END IF
END DO
WRITE(*, *) ' Not converged in 100 times iteration. Try another start value.'
STOP
END PROGRAM erf0


25:デフォルトの名無しさん
06/01/28 00:25:33
実行例
input start value x0 :
0.1
-2.5092300E-02 -2.9582733E-02
-1.0398688E-03 -1.1755297E-03
-1.9108115E-06 -2.1561270E-06
-6.5146001E-12 -7.3509384E-12
Root found x= -6.5146001E-12 : f(x) = -7.3509384E-12
Press any key to continue

input start value x0 :
0.6
0.7874328 -1.7304393E-02
0.7620823 -9.1377224E-05
0.7619469 1.8604648E-08
Root found x= 0.7619469 : f(x) = 1.8604648E-08
Press any key to continue

input start value x0 :
1.3
4.483723 50.93219
3.281823 14.80570
2.510043 4.213046
2.034683 1.139575
1.769606 0.2661909
1.657068 3.9245199E-02
1.633732 1.5385470E-03
1.632740 2.6429243E-06
1.632738 7.6136075E-08
Root found x= 1.632738 : f(x) = 7.6136075E-08
Press any key to continue


26:助けてください(;-;)
06/01/28 00:26:24
解いてくれてありがとうございます!!これはfortranでプログラム組む事できるんですか??

27:助けてください(
06/01/28 00:27:08
うわ!!ありがとうございます!!!

28:デフォルトの名無しさん
06/01/28 00:31:07
つまり、求める解は、
x = 0.0, x = 0.7619469, x = 1.632738
の3つ。

ニュートン法の収束条件判定はいい加減につくった。
自分で工夫してくれ。
単精度実数では、1.0e-7~8 位が数値誤差だ。




この位の宿題は自分でやるです!
間違ってても、謝罪も賠償もしないです!

29:デフォルトの名無しさん
06/01/28 00:37:08
>>26
FORTRAN90で書いたから、理解したうえでFORTRAN77に直してくれw

なお、ここではERF(x)は、IMSLのライブラリ外部関数を使っている。
その辺の特殊性がちょっとある。

fderiv は fの微分だ。

f(x)のグラフを書いて理解を深めてくれ。

もう寝るです。

30:助けてください(
06/01/28 00:38:44
ありがとうございましたmm

31:20
06/01/28 01:03:41
ごめんw ちょっと間違ってたです。

>>20
ERF(x)の微分から 2/√π が抜けてるwww
微妙に式が変わってくるが結論はあんまり変わらないから許せw



32:お願いします;;
06/01/28 16:30:46
<デフォルトの名無しさん><5>さんへ
閏年の問題ありがとうございます!!
参考になりました!!

33:デフォルトの名無しさん
06/01/28 17:56:03
>>32
おまえ、いろんなとこにマルチしまくっているな。

Yahooの方にもいやがった。
トップ > コンピュータとインターネット > プログラミング言語 > 全般 > ◆◆◆プログラミング質問箱◆◆◆

34:デフォルトの名無しさん
06/01/28 18:06:14
YahooID検索したら、ここにもいたw

トップ > コンピュータとインターネット > ソフトウェア > オペレーティングシステム > Windows > Windows XP > fortranわかる人いますか?

何故XPトピ?

35:デフォルトの名無しさん
06/01/28 19:26:52
すいません、ちょっとお力をお貸しください。

曲線y=xlog(a+x)と直線y=axとで囲まれた部分の面積の近似値を
連立方程式の a+3b=1 2a-b=2 を行列計算を用いて求めたa,bを使い答えを求めよ。

というプログラミングなのですがうまくできません。
行列計算においてdo文の使い方がいまいち…
もしよろしければどなたかご助言を。


36:デフォルトの名無しさん
06/01/28 19:32:34
>>33
向こうでの回答はFORTRAN66風味だったです。
ついでにだれかFORTRAN77風に書いてみやがれです。

>>32もマルチポストする気力があるなら、自分で解きやがれです。
サブルーチン2個とメインルーチン1個、それぞれ核心部分は10行程度なんだから、
たいしたこと無いだろです。よく学んでおくです。

37:デフォルトの名無しさん
06/01/28 19:37:37
>>35
そもそも、bが式に現れていない。問題を写し間違えていないか?
それとも単に、連立方程式を解いて出たaだけを使うという問題のなのか?

問題の背景とここまでやった自分の考えを述べなさい。

38:35
06/01/28 20:00:48
>>37
ええと、問題はこれで合っています。
bは使わないで、問題を解いて出てきたaのみで計算するようです。

行列計算ということなのでガウス消去法を用いて解こうとし
それぞれの座標に
a(1,1)=1
というように値を当てはめていったのですが、その後計算するにあたって
fortranのプログラムにおいてガウス消去法がどのようにやるのか分からず、今に至ります。

39:デフォルトの名無しさん
06/01/28 20:02:47
とりあえずフォートランに多相型っていつはいるんですか?

40:デフォルトの名無しさん
06/01/28 20:10:27
>>39
polymorphismならFORTRAN2003でそれなりに入っている。



41:デフォルトの名無しさん
06/01/28 20:10:42
>>35
まず、両方の曲線の交点が2つあって、その間の面積は解析的に求めることは簡単ですね。
で、その解析的に出た面積はaが分かれば数値的に出てくると。
そのaは二元連立方程式の解の一つであると。

二元連立方程式はプログラミングしないでも、中学レベルの頭ならすぐ解ける。

どこにプログラミングの必要性があるのか疑問。


42:デフォルトの名無しさん
06/01/28 20:17:28
>>38
いまいち問題設定が見えてこないなwww これが問題のすべてか?背景は無いのか?

ガウスの消去法は、例外処理を考えなければ、手でやるとおりに組めw
例外処理はめんどくさい。

2*2限定なら、高校で習った公式に突っ込んだほうが早い。例外がdet=0で排除できるし。
というか手でといたほうが早いわけだが。

はっきり言って2*2をガウスの消去法で解く意味が分からんw
というわけだから張り切って一般のn次行列について作れw


43:35
06/01/28 20:29:59
>>41
これをプログラミングで解いてレポートして提出しろという課題なんですよ。
一応実際に計算して解くことは出来ました。
しかしプログラミングのほうが分からず…
>>42
ガウスは使わない方がいいのでしょうか。
となると逆行列を用いての計算のほうがやりやすいのでしょうか?

44:デフォルトの名無しさん
06/01/28 20:32:57
>>43
じゃあさ、面積を求める部分をSimpson法やいろんな数値積分法で組んでしまうのはどう?
面積を求める部分をサブルーチンか関数にしてaを渡せばそれなりにプログラミングの形になるし。

45:42
06/01/28 20:49:16
>>43
ガウスの消去法はアルゴリズムが単純なので、小さめの連立方程式を解くための
練習問題としてはよく使われる。行列がsingularな時や誤差が溜まりやすい状況を、
あまり考慮しなければ、比較的楽に作れるだろう。

素朴には、対角要素が全部0で無いなら、1列目からn列目まで、順に各列の対角要素以外を
0にするように行ごとひいていけばいい。手でとくのと同じだ。
対角要素が0や0付近だと例外処理をすることになる。このへんがまんどくさい。

この辺は、問題製作者の意図がどの辺にあるのかで、こだわり度が決まってくるのだが、
この問題設定では、いまいち理解しにくい。

2*2での一般形なら単純に公式を使うほうが早いが、n*nの一般サブルーチンを求めているなら、
ガウス法なり何らかの方法が必要になる。

後半では数値積分をすることになるわけだが、問題として全く意図が読めないので、
これだけの情報では、どの方法をとるべきか答えがたい。


46:35
06/01/28 21:01:23
>>44
うーん、すいません、またtくの初心者ですのでよくわかりません…
助言は非常に嬉しいのですが…すいません。
>>45
今回はどうやら2*2行列のようです。
それでは今回は公式の方を用いて計算することにします。
まぁ、普通の行列のプログラミングもよくわからないのですが…
後半部分はあれだけです。
タダ単に前半で求めた値を代入して面積を求めるだけでして、背景は無いと思います。

47:デフォルトの名無しさん
06/01/28 22:46:30
>>35
まあ、適当にやれw

a = 1.0, b = 0.0

交点
x = 0.0, x = e - 1

積分 f(x) = ax - x log(a + x) (a = 1.0 積分範囲 [0.0, e-1])

解析的な値
-e^2 / 4 + e - 5 / 4

台形公式 短精度 分割数 n=1000 で こんな感じだ。 
誤差は( (e-1)/1000 )^2〜1e-6位なので、まぁもっともな気がする。

numerical integral 0.378981858
analytic value   0.378982186
Program Completed
Press Enter to Continue.


48:デフォルトの名無しさん
06/01/28 22:51:17
PROGRAM vip
IMPLICIT NONE
INTEGER, PARAMETER :: n = 1000
REAL :: x0, x1, s
x0 = 0.0
x1 = EXP(1.0) - 1.0
s = trapez_integ(n, x0, x1)
PRINT *, 'numerical integral', s
PRINT *, 'analytical value ', 5.0 / 4.0 + EXP(1.0)**2 / 4.0 - EXP(1.0)
STOP
CONTAINS
REAL FUNCTION trapez_integ(n, x0, x1)
INTEGER, INTENT(IN) :: n
REAL , INTENT(IN) :: x0, x1
INTEGER :: i
REAL :: s, h
h = (x1 - x0) / REAL(n)
s = 0.0
DO i = 1, n
s = s + 0.5 * ( func(x0 + (i - 1) * h) + func(x0 + i * h) )
END DO
trapez_integ = s * h
RETURN
END FUNCTION trapez_integ
!
REAL FUNCTION func(x)
REAL, INTENT(IN) :: x
REAL, PARAMETER :: a = 1.0
func = a * x - x * LOG(a + x)
RETURN
END FUNCTION func
END PROGRAM vip

49:デフォルトの名無しさん
06/01/28 22:56:38
訂正

解析的な値
-e^2 / 4 + e - 5 / 4

e^2 / 4 - e + 5 / 4

正負逆だったw

交点が 0.0 と e-1 なのは func(0.0) と func(e-1) を書かせて確かめるのもよし。

前半は、2*2の逆行列の公式でいいだろう。
Ax=b 
x=(A^-1)b
で一応行列の積も入るし。

50:デフォルトの名無しさん
06/01/29 12:21:03
ひょっとして、数値計算プログラム作るなら
Cよりフォートランの方が高速コードを生成しますか?

51:デフォルトの名無しさん
06/01/29 13:06:57
>>50
昔はそうだった。最近はそれほどで差が無い。
コンパイラ自体がフロントエンド部とバックエンド部に分けて作られるようになり、
中間言語→機械語生成のバックエンド部は共通になったから。

しかしF90以降ならどう考えてもFORTRANの方がバグが少なく、
簡潔明快なソースでやりたいことを表現できる。

52:デフォルトの名無しさん
06/01/29 13:12:43
速度は?

53:デフォルトの名無しさん
06/01/29 13:26:03
スパコンのレベルになると、事実上FORTRANしか選択肢がなかったりする

でもそこらのパソコンでやるなら、Cで書こうがFORTRANで書こうが大して変わらんね。
g77とかはほとんど内部でf2cを動かして、C言語に変換してからgccでコンパイルしてるし

54:デフォルトの名無しさん
06/01/29 13:44:16
>>52
コンパイラメーカーが同じなら、速度もたいして変わらない。
要するに、ソースを中間言語に直すところでの最適化の違いだけ。
中間言語から機械語へは共通化している。
ただ、FORTRANの方が制約がきつい分コンパイラーの判断できる最適化は
Cより高いとされている。

>>53
スパコンに関して言えば、90年代はじめごろまではベクトル化コンパイラーはFORTRAN
だけで、Cコンパイラーがあっても非ベクトル型だった。しかし、その後Cコンパイラーも
ベクトル化コードを吐くようになってきた。

55:デフォルトの名無しさん
06/01/29 14:03:46
FORTRAN90/95はコードを書いていて気持ちいい。
柔軟性と制約のバランスが取れていてぴったりくる手袋や靴下のよう。

イラつくのは宣言部が重くなりがちな点くらいか。

みなにお勧めする。

56:デフォルトの名無しさん
06/01/29 14:09:26
じゃあ、特に速度を要求する部分を、書けるならアセンブリで、書けないならフォートランで、
んで普段はC++で書くとかって感じですかね

57:デフォルトの名無しさん
06/01/29 14:26:58
MPIのルーチンはCとFortranしか提供されていないが、微妙にCより。

C++は動的バインディングのせいでCにくらべて10〜20%遅くなるという話だ。
次期FORTRANにOOP機能が付け加わったらどうなることか。

58:デフォルトの名無しさん
06/01/29 19:40:26
IMPLICIT REAL*8 (A-H,O-Z)
WRITE(*,*) 'あ'
N=3
COMPLEX*16 A
END


でコンパイルしようとして

$ g77 test3.f
test3.f: In program `MAIN__':
test3.f:2:
WRITE(*,*) 'あ'
1
test3.f:4: (continued):
COMPLEX*16 A
2
Statement at (2) invalid in context established by statement at (1)

って言われるんですけど、何がいけないんでしょう?
教えてください

59:デフォルトの名無しさん
06/01/29 20:22:25
>>58
COMLEX文は宣言文だから、最初の実行文よりも前に置かれねばならない。
このばあい、WRITE文の前。

それで怒られている。

60:デフォルトの名無しさん
06/01/30 10:41:15
2006年の1月から3月までのカレンダーを作れという課題がでました。

・・・・・・・・・・・
do 20 t=1,3
if(t.eq.2)then
b=28
else
b=31
end if
・・・・・・・・・
・・・・・・・・・
write(6,*)" ",t,"月        "
write(6,*)" 日 月 火 水 木 金 土  "
・・・・・・・・・・
・・・・・・・・・・
write(6,100)(a(i),i=6,b)
100 format(7i5)
・・・・・・・・・・
・・・・・・・・・・

・・・・の部分がわかりません。よろしくお願いします。

61:デフォルトの名無しさん
06/01/30 11:08:08
いまどきフォートラン使ってる人って、どういう立場だなんだ。

趣味? 仕事? 学校? 環境は? 作ってるアプリの種類は?

ちと興味ある。

62:デフォルトの名無しさん
06/01/30 11:24:39
>>61
失敬な!
FORTRANは現役ですよ!いまやモジュール型言語として完成。次はOOPw
INTELが直接開発・販売しているコンパイラはC++とFORTRANだけですよん。

63:デフォルトの名無しさん
06/01/30 12:37:58
いやだから何作ってるのか教えてよ。現役なのは存じてるから

64:デフォルトの名無しさん
06/01/30 15:13:53
FORTRANを科学技術計算以外で使う人なんているの?

65:デフォルトの名無しさん
06/01/30 15:35:34
ということはここに居る人みんな科学者!?

66:デフォルトの名無しさん
06/01/30 16:43:06
>>61, >>63
俺は固体の量子計算。学校の研究です。
アプリっていうとGUIを想像してしまうけど、そうじゃなくて、
コマンドラインで実行するようなプログラム。

>>65
たぶんそうだと思う。
あとガッコウの宿題ができない人々…。


67:デフォルトの名無しさん
06/01/30 17:29:24
Fortran使って構造計算してる土建屋も科学者か?

68:デフォルトの名無しさん
06/01/30 21:27:22
FORTRAN mp3 player
URLリンク(www.pbase.com)

ワロス


69:デフォルトの名無しさん
06/01/30 21:55:33
fortran mp3 encoder
URLリンク(members.at.infoseek.co.jp)



70:デフォルトの名無しさん
06/01/30 21:56:14
ところで、fortranに高階関数はいつ入るんですか?

71:デフォルトの名無しさん
06/01/30 22:25:00
>>70
高階かんすうってなんですか?おしえてかださい

72:デフォルトの名無しさん
06/01/31 01:46:45
            ,,.  -‐―- 、
          / _..-‐―- 、  .ヽ
         / /-‐―-..__..ヽ _..、ヽ
        .!三三三´    `´  : i
        |.!三三   __   ̄ _,. .;. i 暗黙の型宣言を守らない
        |ii i!"  ´ ェェ` l´rェ、 i..!      >>60は地獄へ逝く
         !i、 i!:::..... - ̄r  ヽ  .!
         !ii!.|:::::::  / `−´ヽ /
  _,,.. ===="/%;;;::::.   '∠ニニ」  /、            _ -‐‐-、
/´       | .% :::;;;;: ,,,.,,,:. ⌒: ./´  ̄==--  _,,.. -‐‐'"´ _,..-‐ ´
        |  %  ::::::;;;;;;;;;:.""%./      r" (<二 ̄ ̄>
       .|  %。゚+。  。+% /       `ヽiノ<二. ̄ ̄>
       |   %+*゚+。*゚%/         ( <二 ̄ ̄>




73:デフォルトの名無しさん
06/01/31 02:26:28
>>60 特別にFORTRAN77で書いてやったです。 感謝しやがれです。
PROGRAM calndr
C calender for Jan to Mar 2006
INTEGER ia(31)
CHARACTER*30 fmt
mon1 = 1
C Jan 1 2006 is Sunday (= 1)
DO 10 i = 1, 31
ia(i) = i
10 CONTINUE
DO 20 it = 1, 3
IF (it .eq. 2) THEN
ib=28
ELSE
ib=31
END IF
WRITE(fmt, '(a, i2, a, i1, a)')
1 '(', 4 * (mon1 - 1), 'x, ', (8 - mon1), 'i4, / (7i4))'
mon1 = MOD(ib + mon1 - 1, 7) + 1
WRITE(6, *) ' ', it, '月        '
WRITE(6, *) ' 日 月 火 水 木 金 土  '
WRITE(6, fmt) (ia(i), i = 1, ib)
WRITE(6, *)
20 CONTINUE
STOP
END

>>73


74:73
06/01/31 02:28:02
実行結果です。

1 月      
 日 月 火 水 木 金 土
1 2 3 4 5 6 7
8 9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30 31

2 月      
 日 月 火 水 木 金 土
1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28

3 月      
 日 月 火 水 木 金 土
1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30 31

Press any key to continue

75:73
06/01/31 02:31:04
画面が崩れたので、専ブラ用のアンカーを書いてやるです。 
>>73,75

FORMATの技巧が高度すぎて、>>60はもちろん>>60に出題した
どんくさいチビ人間には理解できないだろです!

76:73
06/01/31 02:31:53
専ブラ用のアンカー間違ったです。 
>>73,74


77:デフォルトの名無しさん
06/01/31 03:14:40
>>71 関数を値として引数に取れる関数だよ。

78:デフォルトの名無しさん
06/01/31 07:00:31
関数のポインタ?

79:デフォルトの名無しさん
06/01/31 09:50:36
>>77
単にexternalなら今でも出来るが、関数ポインターならFORTRAN2003で可能に。

80:デフォルトの名無しさん
06/01/31 16:21:20
今日の宿題はまだか?w

81:デフォルトの名無しさん
06/01/31 22:44:21
ところで型推論っていつFortranにはいるんですか?

82:デフォルトの名無しさん
06/01/31 23:40:14
>>81
型推論ってなんですか?

83:デフォルトの名無しさん
06/02/01 00:18:53
すみません、質問なのですが。
フォートランで作った実行ファイルとソースファイルから
そのソースでcallから呼び出されてる関数
のライブラリ(モジュール?)の場所を探す方法はありますか?

というのも、プログラムを使っていたUNIXマシンから別のUNIXマシンにプログラムを移したところ
マルチバイト文字が正しくありません。とエラーが出てしまいます。
なのでソースがあるので別のマシンでコンパイル、リンクしようと思ったのです。
そうするとソースの中にcallで呼び出されたサブルーチン記述されていない名前の関数が
あったのでリンクのときにこのモジュールを指定しなければいけないんですよね?
その場所がわからずコンパイルでエラーが出てしまいます。

長文で申し訳ないのですがよろしくお願いします。


84:デフォルトの名無しさん
06/02/01 00:56:50
>>83
エラーメッセージを見ていないので状況がいまいちつかめないが、
MAP出力とかのコンパイラ(リンカ、ローダー)オプションで、
必要な情報が得られるかもしれない。



85:83
06/02/01 20:00:59
>>84
レスありがとうです。

いろいろ悩みましたが
コマンド nm を使ってライブラリを全部grepに渡して関数名を検索したら出てきました。
新しいマシンに移してコンパイルリンクでき、実行も出来ました。


86:デフォルトの名無しさん
06/02/01 22:34:43
今日は宿題無いのか?

87:こんにちは
06/02/01 22:44:03
日付を2つ入力して,その間に何日あるかを書き出すプログラム(Nannichi)を作れ。実行例はつぎの通りである. 日付を2つ入力して下さい.
何年? 2004
何月? 7
何日? 1
何年? 2005
何月? 7
何日? 1
日数は 365日です.
これの答えを教えてください お願いします

88:デフォルトの名無しさん
06/02/01 22:53:21
>>87
>>5-6をみろ。すでにより一般的なうるう年を考慮したサブルーチンが存在している。
それをちょっといじればよい。

質問があれば、日本語でおkwwwww

89:デフォルトの名無しさん
06/02/02 07:55:46
>>82 型を宣言しなくても勝手に正しく推論してくれるしくみだよ。

いつはいるの?

90:デフォルトの名無しさん
06/02/02 10:02:57
>>89
率直に言ってそれは下らないアイデアだと思う。
書きなぐるスクリプト言語ではよくあるが、なんんというかだらしなさのあらわれだ。

どっちにしろ暗黙の型宣言があるし。


91:デフォルトの名無しさん
06/02/03 01:39:04
ROZEN MAIDEN ・・・ orz

92:デフォルトの名無しさん
06/02/03 23:41:27
ところで、Fortranにproof-carring-codeはいつ入るの?

93:デフォルトの名無しさん
06/02/03 23:49:39
>>92
proof-carring-codeってなんでちゅか?ばぶー

94:デフォルトの名無しさん
06/02/04 00:46:34
だれかこのプログラミングを教えてください。お願いします。
x=0.1 , 0.2 , 0.3 , 0.4 , 0.5
y=1.228 , 1.005 , 0.823 , 0.674 , 0.552
近似式 y=a*exp(bx)
上の関数とデータで最小自乗法で近似した式を求めよ。
この問題のプログラムがわかりません・・・。お願いします。

95:デフォルトの名無しさん
06/02/04 01:35:39
>>94
PROGRAM vip
IMPLICIT NONE
INTEGER, PARAMETER :: n = 5
REAL, PARAMETER :: x0(n) = (/0.1 , 0.2 , 0.3 , 0.4 , 0.5 /), &
y0(n) = (/1.228, 1.005, 0.823, 0.674, 0.552/)
REAL :: y(n), t(n, 2), z(2, 2), zinv(2, 2), v(2)
INTEGER :: i
y = LOG(y0)
t(:, 1) = 1.0
t(:, 2) = x0
z = MATMUL(TRANSPOSE(t), t)
v = MATMUL(TRANSPOSE(t), y)
zinv = ainverse(z)
v = MATMUL(zinv, v)
WRITE(*, *) ' least square fit a=', EXP(v(1)), ' b=', v(2)
STOP
CONTAINS
FUNCTION ainverse(x)
IMPLICIT NONE
REAL, INTENT(IN) :: x(2, 2)
REAL :: ainverse(2, 2)
REAL :: determinant
determinant = x(1, 1) * x(2, 2) - x(1, 2) * x(2, 1)
IF (ABS(determinant) < 1.0e-7) STOP ' singular matrix '
ainverse(1, 1) = x(2, 2) / determinant
ainverse(2, 2) = x(1, 1) / determinant
ainverse(1, 2) = -x(1, 2) / determinant
ainverse(2, 1) = -x(2, 1) / determinant
RETURN
END FUNCTION ainverse
END PROGRAM vip

96:デフォルトの名無しさん
06/02/04 01:36:42
>>95
least square fit a= 1.499271 b= -1.998701
1.227659 1.005252 0.8231379 0.6740159 0.5519093
1.228000 1.005000 0.8230000 0.6740000 0.5520000

97:94
06/02/04 03:35:31
>>95さん >>96さん
ほんとありがとうございます。
月曜にでも学校でやってみます。
また、わからないことがあれば質問させてください。
ありがとうございました。

98:デフォルトの名無しさん
06/02/04 18:37:30
>>95について。
y=a*exp(bx)を次の式に変形し
log(y)=log(a)+bx
yとxの各値を上の式に代入すると
次の5本の式となる。
log(y1)=log(a)*1+b*x1
log(y2)=log(a)*1+b*x2
log(y3)=log(a)*1+b*x3
log(y4)=log(a)*1+b*x4
log(y5)=log(a)*1+b*x5
ここで、
log(yi)がi番目の要素である列ベクトルをy、
i行1列を1、i行2列をxiとする行列をt、
1番目の要素をlog(a)、2番目の要素をbとする
列ベクトルをvとすると、上の5式は、
y=tv
と表せる。
tが1次独立な列から成るので、
vの最小二乗推定は、次の式
(t't)**(-1)t'y
となる定理を>>95は使用していると思います。
(注)t'はtの転置を表す。(t't)**(-1)はt'tの逆行列を表す。)

99:デフォルトの名無しさん
06/02/04 18:56:37
ところでFortranに配列のバウンダリチェックはいつ入るのでしょうか?

100:デフォルトの名無しさん
06/02/04 19:05:46
>>99
subscript check は、debug 用にoptionで入るのが普通。1960年代からある。
しかし、遅くなるから実際の実行時に使う奴はいない。
この辺はスピード優先の数値計算の価値観を反映している。

101:95
06/02/04 19:40:40
オッス!おれ極右!
>>98氏解説d

最近寒いので湯たんぽ代わりにノートパソコンを腹の上において寝ているお。
昨夜そういう寝そべった状態で布団の中で作ったのが>>95だお。
オラはVipperだから間違っていても謝罪も賠償もしないお!

最小二乗法なんて普段使わないから、きょうの昼おっきしてからNumericalRecipesで
式を確かめといたおw 方式はちょっと違うが式おkwww (15.2.4~6)あたり。
URLリンク(www.library.cornell.edu)

FORTRAN77で書くと3倍くらいの長さになるとおもうお。
MATMULやTRANSPOSEはF90にしかない行列積と転置だお。
LOGを取ったのは線形化するためだお。
普通、最小二乗法は線形の式にしか使わないと思ったから、そうしてみたおw

NumericalRecipesは、旧版のC言語版なら訳本が出ているお。
URLリンク(www.amazon.co.jp)
本のコードにはバグがちりばめられているので、本を打ち込むのはお勧めしないお。
金出して買うコードのほうはこまごまと修正されているお。

(^ω^)このくらいソースコードが短いと勉強する気にもなるだろ?

102:デフォルトの名無しさん
06/02/05 07:53:27
>>87
フォートランはわかんないからベーシック
あってるかどうかもわかんない
こんにちは

10 gosub 60
20 nannichi0=nannichi
30 gosub 60
40 print "日数は";nannichi-nannichi0
50 end
60 input "年,月,日",year,month,day
70 if month<3 then year=year-1:month=month+12
80 nannichi=year+int(year/4)-int(year/100)+int(year/400)+int((153*month-457)/5)+day
90 return

103:デフォルトの名無しさん
06/02/05 11:38:27
>>102
解説頼む

104:デフォルトの名無しさん
06/02/05 20:01:35
>>103
まちがいた。ちょと訂正。

80 nannichi=365*year+int(year/4)-int(year/100)+int(year/400)+int(306*(month+1)/10)+day-428

「これは Fairfield公式 といって、西暦1年1月1日からの経過日数を算出する式らしい。
 ただしグレゴリオ暦が採用されたのは1582年10月15日からだが。
 要するにこれを使って2つの日付の差を求めればよい。」
っておばあちゃんがいってた。


105:94
06/02/05 20:12:40
すいませんがこないだのプログラムを
DO文で作っていただけませんか。
何回もお願いして、すいません。

106:デフォルトの名無しさん
06/02/05 20:33:52
コンピューターおばあちゃん♪
コンピューターおばあちゃん♪
イエーイ イエーイ 僕は大好きさ〜♪

107:デフォルトの名無しさん
06/02/05 21:07:36
>>104
把握したw
2月の28日の例外を2月をお尻にすることで回避。
大の月と小の月の入れ替えを1ヶ月の平均日数を30.6日にすることで実現。
かなり絶妙な公式だ。オモシロスw


だが最近流行の明示的なプログラミングから離れているがwwwww

108:デフォルトの名無しさん
06/02/05 21:54:46
>>105 世話が焼けるッス! FORT77で、素朴に書いといたからw
PROGRAM vipper
   PARAMETER (n = 5)
REAL x(n), y(n)
   DATA x /0.1 , 0.2 , 0.3 , 0.4 , 0.5 /
DATA y /1.228, 1.005, 0.823, 0.674, 0.552/
   sumx = 0.0
   sumy = 0.0
   sumxx = 0.0
   sumxy = 0.0
sum = REAL(n)
DO 10 i = 1, n
    sumx = sumx + x(i)
sumy = sumy + LOG(y(i))
sumxx = sumxx + x(i) * x(i)
sumxy = sumxy + x(i) * LOG(y(i))
10 CONTINUE
det = sum * sumxx - sumx * sumx
alog = ( sumxx * sumy - sumx * sumxy ) / det
b = ( sum * sumxy - sumx * sumy ) / det
WRITE(6, *) ' a = ', EXP(alog), ' b = ', b
   END

>>108 インデント

109:やり直しw
06/02/05 22:11:27
PROGRAM vipper
PARAMETER (n = 5)
REAL x(n), y(n)
DATA x /0.1 , 0.2 , 0.3 , 0.4 , 0.5 /
DATA y /1.228, 1.005, 0.823, 0.674, 0.552/
sumx = 0.0
sumy = 0.0
sumxx = 0.0
sumxy = 0.0
sum = REAL(n)
DO 10 i = 1, n
sumx = sumx + x(i)
sumy = sumy + LOG(y(i))
sumxx = sumxx + x(i) * x(i)
sumxy = sumxy + x(i) * LOG(y(i))
10 CONTINUE
det = sum * sumxx - sumx * sumx
alog = ( sumxx * sumy - sumx * sumxy ) / det
b = ( sum * sumxy - sumx * sumy ) / det
WRITE(6, *) ' a = ', EXP(alog), ' b = ', b
C check
WRITE(6, *) ' i x y',
1 ' a * EXP(bx) (dy)^2'
DO 20 i = 1, n
yy = EXP(alog + b * x(i))
WRITE(6, *) i, x(i), y(i), yy, (y(i) - yy)**2
20 CONTINUE
STOP
END

>>109 インデント

110:デフォルトの名無しさん
06/02/05 22:12:16
a = 1.499271 b = -1.998701
i x y a * EXP(bx) (dy)^2
1 0.1000000 1.228000 1.227659 1.1664589E-07
2 0.2000000 1.005000 1.005252 6.3688631E-08
3 0.3000000 0.8230000 0.8231379 1.9023346E-08
4 0.4000000 0.6740000 0.6740159 2.5326941E-10
5 0.5000000 0.5520000 0.5519093 8.2298044E-09
Press any key to continue
>>110 インデント

111:デフォルトの名無しさん
06/02/05 22:28:29
>>94
死ぬほど易しく作ったから理解しろw
>>101の本の公式通りに作り、変数名も大体あわせておいた。
ここで誤差のσは各要素で共通一定と仮定した。この場合キャンセルしあうから式から無視してよし。

本を見れば分かるが、二乗誤差が極値をとる条件を、a,bに関してそれぞれ微分を取ることで求めている。
二乗誤差は極値を1個だけもちかつそれは最小だからこれでおkとなる。
(誤差は外せばいくらでも大きく出来るが、合わせるほうはに限界がある。最高に合いまくって0。)

プログラム解説:
10番LOOPで総和を取る。
その後、2*2の逆行列をかけるのと同じ理屈で、2元連立方程式を解いている。
20番LOOPでどのくらい近似できたか見ている。

求めていないが、一番右側の列の二乗差の総和が最小になっているはずである。

a,bの値を手動で変えて、求めた二乗差の総和が最小であることを確かめてみそ。
その結果を考察に書いておけば+10点の花丸ハンバーグだ!


112:94
06/02/06 00:15:35
みなさんほんと詳しく教えていただきありがとうございます。
ほんと助かりました。

113:デフォルトの名無しさん
06/02/06 12:36:36
ところで、Fortranにオブジェクトデータ構造はいつ入るんですか.

114:sage
06/02/06 13:32:23
>>113
FORTRAN2003。コンパイラーの実装待ち。

115:デフォルトの名無しさん
06/02/06 17:44:37
>>113
オブジェクトデータ構造ってなんですか?
オブイェークト

116:デフォルトの名無しさん
06/02/06 19:53:17
ERK=ER(K)*DCOS(GANMA2(K))-EI(K)*DSIN(GANMA2(K))
EIK=EI(K)*DCOS(GANMA2(K))+ER(K)*DSIN(GANMA2(K))

実数 ERK=ERK*DCOS(GANMA3(K))-EIK*DSIN(GANMA3(K))
虚数 EIK=EIK*DCOS(GANMA3(K))+ERK*DSIN(GANMA3(K))

ERK=ERK*DCOS(GANMA4(K))-EIK*DSIN(GANMA4(K))
EIK=EIK*DCOS(GANMA4(K))+ERK*DSIN(GANMA4(K))

この式をEXPを用いた複素数型に変換できなくこまってます。
どなたか優しい方おしえていただけませんか(><)

117:デフォルトの名無しさん
06/02/06 21:53:43
>>116
やぁ!やさしいお兄さんだよ!援助してあげるから、今日からパパと呼ぶんだよ!
定義: 
ze = er + i ei
zg = COS(Γ) + i SIN(Γ) = EXP(iΓ)

積:
ze * zg = ( er + i ei ) * ( COS(Γ) + i SIN(Γ) )
     = ( er * COS(Γ) - ei * SIN(Γ) ) + i ( ei * COS(Γ) + er * SIN(Γ) )

これより
COMPLEX*8 :: ze, zek
ze = (er, ei)
zek = ze * EXP( CMPLX(0.0d0, gamma2) )

これでOKさ!

ただ細かいことを言うと、倍精度複素数はF77標準規格に存在しないので、
完全に規格に乗っ取ろうとすると複素変数は使えないのさ!
でも、大概のコンパイラは倍精度複素数があるのでOKさ。
EXPとCMPLXは総称名なので、精度は引数に合わせてくれるはず。

でもマニュアルを調べて総称名が、ちゃんと倍精度複素数になってるかどうか
確認したほうがいいかもね。

118:117
06/02/06 22:00:01
おっと
ze = CMPLX(er, ei)
だったね。

さぁ、>>116 おっぱいうpうp
    _  ∩
  ( ゚∀゚)彡 おっぱい!おっぱい!
  (  ⊂彡
   |   | 
   し ⌒J


119:デフォルトの名無しさん
06/02/06 22:39:00
なんだかしりませんが
Fortran2003って普通にかな〜りすごそうなんですけど
ほんとに実装されるんですか?

120:デフォルトの名無しさん
06/02/06 23:17:41
>>118
URLリンク(kahouweb.hp.infoseek.co.jp)

121:デフォルトの名無しさん
06/02/06 23:29:50
>>119
FORTRANの場合スーパーコンピュータやCPUのキラーアプリとして
単独での採算を度外視されて開発されているので、多分実装されるでしょう。

INTELが旧DECのFORTRANコンパイラーチームをCompaqから買い取ったのは
そういう方向からだし。

122:デフォルトの名無しさん
06/02/07 11:50:47
次の宿題もってこい!

123:デフォルトの名無しさん
06/02/07 15:28:52
f90を使い初めてプログラムの複素数化をしているのですがうまくいきません、どなたか助言をおねがいします。
PROGRAM TTT
IMPLICIT NONE
INTEGER :: I,N,M
REAL(8) :: EI,ER(0:2047),EI(0:2047),A,B
REAL(8) :: FACTR,H,PI,RN1,RN2,EP0,C,WL,TINTY,TH,S(0:2047),TT(0:2047),GEFF,ERI,EII
!
PI=4.0*DATAN(DBLE(1.0))
C=3.0E+8
WL=1.55E-6
DO I=0,N-1
FACTR=-H*PI*RN2*RN1*EP0*C/WL
TINTY=(ABS(S(I)*S(I))+ABS(TT(I)*TT(I)))/DBLE(2.0)
TH=FACTR*TINTY
ERI=EXP(-GEFF*H)*(1*DCOS(TH)-1*DSIN(TH))
EII=EXP(-GEFF*H)*(1*DCOS(TH)+1*DSIN(TH))
ER(I)=ERI
EI(I)=EII
OPEN(20,FILE='TH.dat')
WRITE(20,*)ER(I),EI(I)
END DO
END PROGRAM TTT
上のプログラムからSINやCOSをけしてEXP等で表現
し複素数でのプログラムにするにはどうしたらいいでしょうか?

124:デフォルトの名無しさん
06/02/07 16:48:19
>>123
まずその前に・・・・・
C=3.0E+8
WL=1.55E-6
これだと単制度で定数が入るので、ここで数値が腐ってしまう。E→Dにせよ。

吸収係数かなんかの計算か?
発想が逆さまだ。もともとの公式を素直にプログラムするほうがよっぽど楽だぞw

ほいで
(1 + i) * EXP(iθ) = (1 + i) * (cosθ+ isinθ)
          = (cosθ+ isinθ) + (icosθ - sinθ)
          = (cosθ - sinθ) + i(cosθ +sinθ)
ゆえに
ZE = EXP(-geff * h) * (1.0d0, 1.0d0) * EXP(CMPLX(0.0d0, th))
こうじゃねーか?

間違っていてもシラネw


次ページ
最新レス表示
スレッドの検索
類似スレ一覧
話題のニュース
おまかせリスト
▼オプションを表示
暇つぶし2ch

4302日前に更新/404 KB
担当:undef