module integral
  !
  !== 台形公式を用いた数値積分モジュール
  !
  ! Authors::   森川 靖大
  ! Version::   0.3 2008-03-09 morikawa
  ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
  !
  !--
  ! (#-- から #++ までの部分を RDoc は解釈しません.)
  ! "=", "==", ""===" は見出しを表します.
  !
  != 見出しレベル1
  !== 見出しレベル2
  !=== 見出しレベル3
  !++
  !
  ! trapezoid サブルーチンに積分区間に関する情報を与えると,
  ! 計算を行い値を返します. 
  ! 
  !--
  ! モジュール名やメソッド名はそのままモジュールやメソッドへのリンクに
  ! 変換されます.
  !++
  ! 
  !=== 参考資料
  ! 
  ! * http://www.gfd-dennou.org/library/dcmodel
  !   1. {RDoc Fortran 90/95 ソースコード解析機能強化版}[http://www.gfd-dennou.org/library/dcmodel/rdoc-f95]
  !   2. {地球流体電脳davis/rubyワークショップ}[http://www.gfd-dennou.org/library/davis/workshop/2008-03-10/]
  !
  !--
  !== リストの表示に関して
  !
  ! リストは以下のような記号が付いたパラグラフです. 
  !
  !  - '*' もしくは '-' で普通のリスト
  !  - 数字+ピリオドで番号付きリスト
  !  - アルファベット+ピリオドでアルファベットリスト
  !
  !== リンクに関して
  !
  ! http:, mailto:, ftp:, www. で始まるテキストはウェブへのリンクだと
  ! 判別されます.
  !
  ! label[url] の形式でもハイパーリンクが張れます. この場合は lavel が表
  ! 示され, url がリンク先となります. label が複数の単語を含んでいる場合 
  ! (日本語の場合はこっちを使ってください), 中括弧を使い, <em>{multi word
  ! label}[</em>url<em>]</em>としてください. 
  !++
  ! 
  !=== 開発履歴
  ! 
  ! * 0.3  Ruby 用 RDoc チュートリアルを元に作ってみた. 
  ! 

  private
  public trapezoid
contains
  function trapezoid(a,b,n) result(s)
    !
    ! 台形公式を用いた数値積分を行います. 
    !
    ! $ f(x) = x^2 - 2x + 2 $ の積分値を台形公式によって求めます.
    !
    ! * 積分区間は $ [a, b] $ とします.
    ! * 実際に計算する式は以下の通りです.
    !   ここで $ h $ は $ h = \frac{b-a}{n} $ です.
    !
    ! \[
    !    \int_{a}^{b} f(x) dx =
    !      \frac{h}{2} \left\{f(a) + f(b)\right\} + h\sum_{i=1}^{n-1} f(a+ih)
    ! \]
    !
    ! このサンプルプログラムは http://www.gfd-dennou.org/library/dcmodel
    ! にて入手できます. 
    !
    !--
    ! 別のモジュール内のサブルーチンなどへリンクする場合は
    ! "<i>モジュール名</i>#<i>サブルーチン名</i>" と指定します. 
    !++

    real,    intent(in) :: a ! 積分区間の下限
    real,    intent(in) :: b ! 積分区間の上限
    integer, intent(in), optional :: n ! 積分区間の分割数
    real :: s  ! 積分値
    real :: h
    integer :: i, div = 10
    if (present(n)) div = n
    h = (b-a) / div
    s = (f(a)+f(b)) / 2.0
    do i = 1, n-1
      s = s + f(a + i*h)
    end do
    s = h*s
  end function trapezoid

  function f(x) result(res)
    !
    ! $ f(x) = x^2 - 2x + 2 $ の値を返す内部サブルーチンです. 
    !
    real, intent(in) :: x
    real :: res
    res = x**2 - 2*x + 2
  end function f
end module integral
