00001 function dfdx_3rd(x,f,v)
00002 implicit none
00003 real(8) :: dfdx_3rd
00004 real(8) :: x(4),f(4), v
00005 real(8) :: dx01, dx02, dx03, dx10, dx12, dx13
00006 real(8) :: dx20, dx21, dx23, dx30, dx31, dx32
00007 real(8) :: xv0, xv1, xv2, xv3
00008 real(8) :: wex0, wex1, wex2, wex3
00009 integer :: ir0, ir1, ir2, ir3
00010
00011 ir0 = 1
00012 ir1 = 2
00013 ir2 = 3
00014 ir3 = 4
00015 dx01 = x(ir0) - x(ir1)
00016 dx02 = x(ir0) - x(ir2)
00017 dx03 = x(ir0) - x(ir3)
00018 dx12 = x(ir1) - x(ir2)
00019 dx13 = x(ir1) - x(ir3)
00020 dx23 = x(ir2) - x(ir3)
00021 dx10 = - dx01
00022 dx20 = - dx02
00023 dx21 = - dx12
00024 dx30 = - dx03
00025 dx31 = - dx13
00026 dx32 = - dx23
00027 xv0 = v - x(ir0)
00028 xv1 = v - x(ir1)
00029 xv2 = v - x(ir2)
00030 xv3 = v - x(ir3)
00031 wex0 = (xv2*xv3 + xv1*xv3 + xv1*xv2)/(dx01*dx02*dx03)
00032 wex1 = (xv2*xv3 + xv0*xv3 + xv0*xv2)/(dx10*dx12*dx13)
00033 wex2 = (xv1*xv3 + xv0*xv3 + xv0*xv1)/(dx20*dx21*dx23)
00034 wex3 = (xv1*xv2 + xv0*xv2 + xv0*xv1)/(dx30*dx31*dx32)
00035 dfdx_3rd = wex0*f(ir0) + wex1*f(ir1) + wex2*f(ir2) + wex3*f(ir3)
00036
00037 end function dfdx_3rd
