Fortran実装関係

　実装した関数などをまとめていこうかなと

線形代数系

三点の座標から三角形の面積

三点の座標から三角形の面積

!三点の座標から三角形の面積を出す関数
real(8) function areatri(x1,y1,x2,y2,x3,y3)

        implicit none
        
        real(8),intent(in)::x1,y1,x2,y2,x3,y3
        real(8)::a12,a23,a31,s,area

        !各辺の長さ(a12は点1,2の間の辺を表す)
        a12 = sqrt((x1-x2)**2 + (y1-y2)**2)
        a23 = sqrt((x2-x3)**2 + (y2-y3)**2)
        a31 = sqrt((x3-x1)**2 + (y3-y1)**2)

        !ヘロンの公式を用いる
        s = (a12+a23+a31)/2.0
        area = sqrt(s*(s-a12)*(s-a23)*(s-a31))

        !return
        areatri = area
        
    end function

大気分光において等価吸収幅を求める際、スペクトルの形をおおよそ三角形なんじゃね？と近似してしまうのに用いる。

統計関係

最小二乗法（直線）

最小二乗法（直線）

!x(:),y(:)のデータからy=ax+bを求める
!s_a,s_bはa,bの標準不確かさ
subroutine least_sq(x,y,a,b,s_a,s_b)

        implicit none
        
        real(8), allocatable , intent(in):: x(:),y(:)
        real(8) :: a,b,ave_x,ave_y,ave_x2,ave_xy,s_y_p2,s_a_p2,s_b_p2,s_a,s_b
        integer :: n,i

        n = size(x)
        ave_x = sum(x)/real(n,8)
        ave_y = sum(y)/real(n,8)

        ave_x2 = 0
        do i =lbound(x,dim=1),ubound(x,dim=1)
            ave_x2 = ave_x2 + x(i)**2
        end do
        ave_x2 = ave_x2/(real(n,8))

        ave_xy = 0
        do i =lbound(y,dim=1),ubound(y,dim=1)
            ave_xy = ave_xy + x(i)*y(i)
        end do
        ave_xy = ave_xy/(real(n,8))

        a = (ave_xy-ave_x * ave_y)/(ave_x2 - ave_x**2)
        b = ave_y - a*ave_x

        s_y_p2 = 0
        do i = lbound(x,dim=1),ubound(x,dim=1)
            s_y_p2 = s_y_p2 + (y(i)-a*x(i)-b)**2
        end do
        s_y_p2 = s_y_p2/(real(n-2,8))

        s_a_p2 = (s_y_p2)/(real(n,8)*(ave_x2-ave_x**2))
        s_a = sqrt(s_a_p2)
        
        s_b_p2 = (ave_x2*s_y_p2)/(real(n,8)*(ave_x2-ave_x**2))
        s_b = sqrt(s_b_p2)


    end subroutine 

　みんな大好き最小二乗法。a,bの標準不確かさも同時に求めてくれるので重宝している。n>3でしか使えない。系統誤差がある時は、

    s_y_p2 = s_y_p2 + !系統不確かさ^2

をs_a_p2の前にぶっ込んでしまったりする。