00001 function lagint_4th(x,y,v)
00002 use phys_constant, only : long
00003 implicit none
00004 real(long) :: lagint_4th
00005 real(long) :: x(4),y(4), v
00006 real(long) :: dx12, dx13, dx14, dx23, dx24, dx34
00007 real(long) :: dx21, dx31, dx32, dx41, dx42, dx43
00008 real(long) :: xv1, xv2, xv3, xv4, wex1, wex2, wex3, wex4
00009
00010 dx12 = x(1) - x(2)
00011 dx13 = x(1) - x(3)
00012 dx14 = x(1) - x(4)
00013 dx23 = x(2) - x(3)
00014 dx24 = x(2) - x(4)
00015 dx34 = x(3) - x(4)
00016 dx21 = - dx12
00017 dx31 = - dx13
00018 dx32 = - dx23
00019 dx41 = - dx14
00020 dx42 = - dx24
00021 dx43 = - dx34
00022 xv1 = v - x(1)
00023 xv2 = v - x(2)
00024 xv3 = v - x(3)
00025 xv4 = v - x(4)
00026 wex1 = xv2*xv3*xv4/(dx12*dx13*dx14)
00027 wex2 = xv1*xv3*xv4/(dx21*dx23*dx24)
00028 wex3 = xv1*xv2*xv4/(dx31*dx32*dx34)
00029 wex4 = xv1*xv2*xv3/(dx41*dx42*dx43)
00030
00031 lagint_4th = wex1*y(1) + wex2*y(2) + wex3*y(3) + wex4*y(4)
00032
00033 end function lagint_4th
