Hello everyone, Please I have a table of values of 33 rows and 15790 columns(x, f1(x),f2(x),...,f15789(x)). The first contains the values of x and the rest of those of f(X) ( 15789 f(x)). I want to read in this table, do the interpolation and then write the fitted values in an array of 128 grids points.
Code
module interpolation ! start of interpolation module.
implicit none
real(kind = 8),dimension(:),allocatable::x,y
real(kind = 8),allocatable::a(0:n), b(0:n), c(0:n), d(0:n)
integer::n
end module interpolation ! end of interpolation module.
program main ! start of main program.
use interpolation
implicit none
character*(*),parameter::inputdatafile = 'Tableau.dat'
real(kind = 8)::sj, xin, f
integer:: m, j
! count data set.
call count_data(inputdatafile,m)
! allocate x and y arrays.
n = m - 1
allocate(x(0:n),y(0:n))
! load data set.
call load_data(inputdatafile,n,x,y)
x(0:32) = tab1(1:33,1)
do i = 1, 15789
a(0:32) = tab1(1:33, i+1)
! b = 0.0, c = 0.0, d = 0.0
! allocate a, b, c, and d arrays.
allocate(a(0:32), b(0:32), c(0:32), d(0:32))
!
!a = y;
!
! call natural cubic spline subroutine.
call natural_cubic_spline(n,x,a,b,c,d)
mat(1,1:33,i) = b(0:32)
mat(2,1:33,i) = c(0:32)
mat(3,1:33,i) = d(0:32)
enddo
! Interpolate data and write to file.
!open(100, file = 'interpolate1.dat', action = 'write')
xin = x(0)
do while(xin .le. x(n))
write(*,*))xin,f(xin)
! xin = xin + 0.04
enddo
!close(100)
do i = 1, 15789
b(0:32) = mat(1,1:33,i)
c(0:32) = mat(2,1:33,i)
d(0:32) = mat(3,1:33,i)
open(j+90, file = 'interpolate.dat', action = 'write')
do i=1,128
xin=4.0+dble(i-1)*(12.0-4.0)/(128-1)
write(j+90,77) f(xin)
77 format '(128e25.14)'
enddo
enddo
!
deallocate(x,y)
deallocate(a, b, c, d)
!
stop
end program main ! end of main program.
real(kind = 8) function f(xin) ! start of interpolant function f.
use interpolation
implicit none
real(kind = 8),intent(in)::xin
integer:: j, i
!
do i = 0,n-1
if((x(i) .le. xin).and.(xin .le. x(i+1)))then
j = i
exit
endif
enddo
!
f = a(j) + b(j)*(xin - x(j)) + c(j)*(xin - x(j))**2 + d(j)*(xin - x(j))**3
!
return
end function f ! end of interpolant function f.
subroutine load_data(filename,n,x,y) ! start of load_data subroutine.
! load data from filename.
implicit none
character*(*),intent(in)::filename
integer,intent(in)::n
real(kind = 8),dimension(0:n),intent(out)::x,y
integer::i,ios
!
open(100, file = filename, action = 'read')
ios = 0; i = 0
do while((ios .ge. 0).and.(i .le. n))
read(100,*,iostat = ios)x(i),y(i)
if(ios == 0)then
i = i + 1
endif
enddo
close(100)
!
return
end subroutine load_data ! end of load_data subroutine.
subroutine count_data(filename,n) ! start of count_data subroutine.
! This subroutine counts the number of data set in filename.
implicit none
character*(*),intent(in):: filename
integer,intent(out)::n
real(kind = 8)::x,y
integer::ios
!
open(100, file = filename, action = 'read')
ios = 0; n = 0
do while(ios .ge. 0)
read(100,*,iostat = ios)x,y
if(ios == 0)then
n = n + 1
endif
enddo
close(100)
!
return
end subroutine count_data ! end of count_data subroutine.
subroutine natural_cubic_spline(n,x,a,b,c,d) ! start of natural_cubic_spline subroutine.
implicit none
integer,intent(in)::n
real(kind = 8),intent(in)::x(0:n)
real(kind = 8),intent(inout):: a(0:n)
real(kind = 8),intent(out):: b(0:n), c(0:n), d(0:n)
integer:: i, j
real(kind = 8):: h(0:n), alpha(0:n), l(0:n), mu(0:n), z(0:n)
!
h = 0.0; alpha = 0.0
!
! step 1:
do i = 0,n-1
h(i) = x(i+1) - x(i)
enddo
!
! step 2:
do i = 1,n-1
alpha(i) = (3.0/h(i))*(a(i+1) - a(i)) - (3.0/h(i-1))*(a(i) - a(i-1))
enddo
!
! step 3:
l(0) = 1.0; mu(0) = 0.0; z(0) = 0.0;
!
! step 4:
do i = 1,n-1
l(i) = 2.0*(x(i+1) - x(i-1)) - h(i-1)*mu(i-1)
mu(i) = h(i)/l(i)
z(i) = (alpha(i) - h(i-1)*z(i-1))/l(i)
enddo
!
! step 5:
l(n) = 1.0; z(n) = 0.0; c(n) = 0.0
!
! step 6:
do j = n-1,0,-1
c(j) = z(j) - mu(j)*c(j+1)
b(j) = (a(j+1) - a(j))/h(j) - h(j)*(c(j+1) + 2.0*c(j))/3.0
d(j) = (c(j+1) - c(j))/(3.0*h(j))
enddo
!
return
end subroutine natural_cubic_spline ! end of natural_cubic_spline subroutine.
Sample data
x f1(x) f2(x),..., f1579(x)
4.000 0.55977887068214E-01
4.250 0.33301830925674E-01
4.500 0.19357033058574E-01
4.750 0.10851672673603E-01
5.000 0.57601109719013E-02
5.250 0.27793714306045E-02
5.500 0.10812628240276E-02
5.750 0.15097129824141E-03
6.000 -0.32616649818477E-03
6.250 -0.54065650568596E-03
6.500 -0.60723892717492E-03
6.750 -0.59463219747272E-03
7.000 -0.54352144592126E-03
7.250 -0.47743610751912E-03
7.500 -0.40939905472145E-03
7.750 -0.34605270924249E-03
8.000 -0.29025099403574E-03
8.250 -0.24269312486307E-03
8.500 -0.20294374360102E-03
8.750 -0.17005393420182E-03
9.000 -0.14292238222802E-03
9.250 -0.12048981961726E-03
9.500 -0.10182978503448E-03
9.750 -0.86178058019242E-04
10.000 -0.72928575326720E-04
10.250 -0.61613399268902E-04
10.500 -0.51877242408502E-04
10.750 -0.43452338487951E-04
11.000 -0.36136444117556E-04
11.250 -0.29774951350710E-04
11.500 -0.24247094860570E-04
11.750 -0.19455757052863E-04
12.000 -0.15320203216670E-04
results
real(kind = 8),allocatable::a(0:n), b(0:n), c(0:n), d(0:n)
1
Error: Explicit shaped array with nonconstant bounds at (1)
main1.f90:9:5:
use interpolation
1
Fatal Error: Can't open module file ‘interpolation.mod’ for reading at (1): Aucun fichier ou dossier de ce type
compilation terminated.
