[Fortran] 纯文本查看 复制代码
Program www_fcode_cn
use SplineMod
Implicit None
Integer , parameter :: N = 10
Integer , parameter :: M = 37
Real :: x(N) , y(N) , xx(M) , yy(M)
integer :: i
Open( 12 , File = 'a.txt' )
Open( 13 , File = 'b.txt' )
Do i = 1 , N
Read( 12 , * ) x(i) , y(i)
End Do
call SplineMod_Init( N, x , y )
Do i = 1 , M
xx(i) = (i-1) * 0.5
yy(i) = SplineMod_Interp( N , x , y , xx(i) )
End Do
Do i = 1 , M
Write( 13 , * ) xx(i) , yy(i)
End Do
End Program www_fcode_cn
Module SplineMod
Implicit None
Integer , private :: i , N
Real , private , Allocatable :: S2(:) , Dely(:) , S3(:)
Contains
Subroutine SplineMod_Init( iN , x , y )
Integer , Intent( IN ) :: iN
Real , Intent( IN ) :: x( iN ) , y( iN )
Real :: B( iN ) , C( iN ) , H( iN ) , H1( iN ) , Delsqy
Integer jj , N1
N = iN
N1 = N - 1
Allocate( S2( N ) )
Allocate( Dely( N ) )
Allocate( S3( N ) )
Do i = 1 , N1
H( i ) = x( i+1 ) - x( i )
if ( abs(H(i)) < 1./3600 ) H(i) = 1./3600
Dely( i ) = ( y(i+1) - y(i) ) / H( i )
End Do
Do i = 2 , N1
H1( i ) = H( i-1 ) + H( i )
B( i ) = 0.5 * H(i-1) / H1( i )
Delsqy = ( Dely(i) - Dely(i-1) ) / H1( i )
S2( i ) = 2.0 * Delsqy
C( i ) = 3.0 * Delsqy
End Do
S2( 1 ) = 0.0
S2( N ) = 0.0
Do jj = 1 , 26
Do i = 2 , N1
S2(i) = (C(i)-B(i)*S2(i-1)-(0.5-B(i))*S2(i+1)-S2(i))*1.0717968+S2(i)
End Do
End Do
Do i = 1 , N1
S3( i ) = ( S2(i+1) - S2(i) ) / H( i )
End Do
End Subroutine SplineMod_Init
Real Function SplineMod_Interp( iN , x , y , T )
Integer , Intent( IN ) :: iN
Real , Intent( IN ) :: T
Real , Intent( IN ) :: x( iN ) , y( iN )
Integer i
Real :: ht1 , ht2 , Delsqs
i = 1
if( ( T - x(i) ) <= 0.0 ) goto 17
if( ( T - x(N) ) < 0.0 ) goto 57
goto 59
56 if( ( T - x(i) ) < 0.0 ) goto 60
if( ( T - x(i) ) == 0.0 ) goto 17
57 i = i + 1
GOTO 56
59 i = N
60 i = i - 1
17 HT1 = T - x(i)
HT2 = T - x(i+1)
Delsqs = ( 2.0 * S2(i) + S2(i+1) + HT1 * S3(i) ) / 6.0
SplineMod_Interp = y(i) + HT1 * Dely( i ) + HT1 * HT2 * Delsqs
End Function SplineMod_Interp
Subroutine SplineMod_UnInit()
DeAllocate( S2 )
DeAllocate( Dely )
DeAllocate( S3 )
End Subroutine SplineMod_UnInit
End Module SplineMod
