00001 program extrapolation
00002 implicit none
00003 integer :: i, k, n = 6
00004 real(8) :: xi(4),yi(4), yin(4,6), yex(4)
00005 real(8) :: v, x(3),y(3)
00006 real(8), external :: lagint_3rd
00007
00008 open(1,file='data.dat',status='old')
00009 do i = 1, 4
00010 read(1,*) xi(i), (yin(i,k), k = 1, n)
00011 end do
00012 close(1)
00013
00014 v = 0.0
00015
00016 do k = 1, n
00017 yi(:) = yin(:,k)
00018 x(1:3) = xi(2:4)
00019 y(1:3) = yi(2:4)
00020 yex(k) = lagint_3rd(x,y,v)
00021 end do
00022
00023 write(6,*) (yex(k), k = 1, n)
00024 open(1,file='output_3rd.dat',status='unknown')
00025 do i = 1, 4
00026 write(1,'(1p,9e18.10)')xi(i), &
00027 & (dabs(1.0d0 - yin(i,k)/yex(k))*100.0, k = 1, n)
00028 end do
00029 close(1)
00030 open(1,file='output.dat',status='unknown')
00031 do i = 1, 3
00032 write(1,'(1p,9e18.10)') &
00033 & dsqrt((xi(i)/xi(4))**2 - 1.0d0)*xi(4), &
00034 & (dabs(1.0d0 - yin(i,k)/yin(4,k))*100.0, k = 1, n)
00035 end do
00036 close(1)
00037
00038 end program extrapolation
00039
00040 function lagint_3rd(x,y,v)
00041 implicit none
00042 real(8) :: lagint_3rd
00043 real(8) :: x(3),y(3), v
00044 real(8) :: dx12, dx21, dx13, dx31, dx23, dx32
00045 real(8) :: xv1, xv2, xv3, wex1, wex2, wex3
00046
00047 dx12 = x(1) - x(2)
00048 dx13 = x(1) - x(3)
00049 dx23 = x(2) - x(3)
00050 dx21 = - dx12
00051 dx31 = - dx13
00052 dx32 = - dx23
00053 xv1 = v - x(1)
00054 xv2 = v - x(2)
00055 xv3 = v - x(3)
00056 wex1 = xv2*xv3/(dx12*dx13)
00057 wex2 = xv1*xv3/(dx21*dx23)
00058 wex3 = xv1*xv2/(dx31*dx32)
00059
00060 lagint_3rd = wex1*y(1) + wex2*y(2) + wex3*y(3)
00061
00062 end function lagint_3rd
