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