戻る
エラーデータを線形内挿で埋める
(Fortran90)
プログラムのダウンロード
研究室の皆さん
玉川です。エラーデータの処理をするサブルーチンを作ってみました。
特に、豊吉君にお勧めです。
以下につけたようなもので
最初の program 〜 end までは、テスト用のメインプログラムです。
その下以降を自分のソフトに組み込んで使って下さい。
使い方は、
データの入った配列 x1 ( real x1(n) )
不良データを埋めたデータをいれる配列 x2 ( real x2(n) )
データ数 n integer
不良データ数 errnum integer
を用意し、
call interpol( x1, x2, n, errnum )
と呼ぶだけです。
埋める処理は、ごく簡単に、不良データの前後のデータ１つだけを見て
線形内挿しています。 (端データが不良の場合は、残った側の２つを見て)
全部不良データだったらギブアップして、何もしません。
参考にして下さい。３０分と時間を掛けてないので、間違っているかも知れません。
********************************************************
間違っていました。訂正しておきます。2007.12.5 夜
間違っていました。訂正しておきます。2008.1.7 夜
********************************************************
一応動いているところを書いておきます。左の出力が元データです。
--- 動作例 -----
tama@tama[139] !./a
./a.out < SAMPLEDATA
number of error data = 3
-9999.000 1.000000
2.000000 2.000000
3.000000 3.000000
4.000000 4.000000
-9999.000 5.000000
-9999.000 6.000000
7.000000 7.000000
8.000000 8.000000
9.000000 9.000000
10.00000 10.00000
--- その２ (２回目の修正までだめだったパターン) ----
tama@tama[107] ./a.out < SAMPLEDATA
number of error data = 8
-9999.000 1.000000
2.000000 2.000000
-9999.000 3.000000
4.000000 4.000000
-9999.000 5.000000
-9999.000 6.000000
7.000000 7.000000
-9999.000 8.000000
9.000000 9.000000
-9999.000 10.00000
--------------------------------------------------------------------------------
玉川一郎 (e-mail: tama@green.gifu-u.ac.jp)
--------------------------------------------------------------------------------
!
! 2008.1.7 debuged I. Tamagawa
!
program test
implicit none
integer, parameter :: n=10
real :: x1(n), x2(n)
integer :: i, errnum
do i=1, n
read(*,*) x1(i)
end do
call interpol( x1, x2, n, errnum )
write(*,*) 'number of error data = ', errnum
do i=1, n
write(*,*) x1(i), x2(i)
end do
end
logical function ok( x, count )
implicit none
real :: x
integer :: count
if( x < -9900 ) then
ok = (5<0)
count = count + 1
else
ok = (5>0)
endif
end function
subroutine interpol( x1, x2, n, errnum )
implicit none
real, intent(in) :: x1(n)
real, intent(out) :: x2(n)
integer, intent(in) :: n
integer, intent(out) :: errnum
! x1 ==> x2 (error data are linearly interpoleted)
integer :: i, err1st, errlast, i1, i2
integer :: count
logical, external :: ok
count = 0
errnum = 0
do i=1, n
if( ok(x1(i), count) ) then
x2(i) = x1(i)
if( count > 0 ) then
errlast = i - 1
call linear_interpol( x1, x2, err1st, errlast, n )
errnum = errnum + (errlast-err1st+1)
count = 0
end if
else
if( count == 1 ) then
err1st = i
end if
end if
end do
if( count > 0 ) then
errlast = n
call linear_interpol( x1, x2, err1st, errlast, n )
errnum = errnum + count
end if
end subroutine
subroutine find_next_ok( x, n, errlast, i2 )
! find next good data in x(i) i > errlast+1
implicit none
integer, intent(in) :: n
real, intent(in) :: x(n)
integer :: errlast, i2
integer :: i, dum
logical, external :: ok
!write(*,*) 'find next ok data'
dum=0
i = errlast+2
do while( (.NOT. ok(x(i), dum)) .AND. (i<=n) )
!write(*,*) 'i in while loop', i
i=i+1
end do
if( i == n ) then
if( ok(x(i), dum) ) then
i2 = i
else
i2 = errlast + 1
end if
else
i2 = i
end if
!write(*,*) ' i2 = ', i2
end
subroutine find_last_ok( x, n, err1st, i1 )
! find last good data in x(i) i < err1st-1
implicit none
integer, intent(in) :: n
real, intent(in) :: x(n)
integer :: err1st, i1
integer :: i, dum
logical, external :: ok
!write(*,*) 'find before ok data'
dum=0
i = err1st-2
do while( (.NOT. ok(x(i), dum)) .AND. (i>0) )
!write(*,*) 'i in while loop', i
i=i-1
end do
if( i <= 0 ) then
i1 = err1st - 1
else
i1 = i
end if
!write(*,*) ' i1 = ', i1
end
subroutine linear_interpol( x1, x2, err1st, errlast, n )
implicit none
integer, intent(in) :: n
real, intent(in) :: x1(n)
real :: x2(n)
integer :: err1st, errlast
integer :: i1, i2, i
real :: a, b
! write(*,*) ' interpol (err1st and errlast)', err1st, errlast
if( err1st == 1 .AND. errlast == n ) then ! NO GOOD DATA
write(0,*) 'Give up to interporate error data, because all', &
'data are error data'
do i=1,n
x2(i)=x1(i)
end do
return
end if
if( err1st > 1 .AND. errlast < n) then
i1 = err1st - 1; i2 =errlast +1
else if( err1st == 1 .AND. errlast < n-1 ) then
i1 = errlast +1
call find_next_ok( x1, n, errlast, i2 )
errlast = i2 - 1
else if( errlast == n .AND. err1st > 2 ) then
i2 = err1st -1
call find_last_ok( x1, n, err1st, i1 )
err1st = i1+1
end if
if( i2 > i1 ) then
a = (x1(i2)-x1(i1))/(i2-i1)
do i=err1st, errlast
x2(i) = x1(i1) + a*(i-i1)
end do
else ! only one datum is ok (maybe)
write(0,*) 'Only one datum is OK. Number of error data is wrong!'
do i=1, n
x2(i) = x1(i1)
end do
end if
end subroutine