00001 function dfdx_4th(x,f,v)
00002 use phys_constant, only : long
00003 implicit none
00004 real(long) :: dfdx_4th
00005 real(long) :: x(5),f(5), v
00006 real(long) :: dx01, dx02, dx03, dx04
00007 real(long) :: dx10, dx12, dx13, dx14
00008 real(long) :: dx20, dx21, dx23, dx24
00009 real(long) :: dx30, dx31, dx32, dx34
00010 real(long) :: dx40, dx41, dx42, dx43
00011 real(long) :: xv0, xv1, xv2, xv3, xv4
00012 real(long) :: wex0, wex1, wex2, wex3, wex4
00013 integer :: ir0, ir1, ir2, ir3, ir4
00014
00015 ir0 = 1
00016 ir1 = 2
00017 ir2 = 3
00018 ir3 = 4
00019 ir4 = 5
00020 dx01 = x(ir0) - x(ir1)
00021 dx02 = x(ir0) - x(ir2)
00022 dx03 = x(ir0) - x(ir3)
00023 dx04 = x(ir0) - x(ir4)
00024 dx12 = x(ir1) - x(ir2)
00025 dx13 = x(ir1) - x(ir3)
00026 dx14 = x(ir1) - x(ir4)
00027 dx23 = x(ir2) - x(ir3)
00028 dx24 = x(ir2) - x(ir4)
00029 dx34 = x(ir3) - x(ir4)
00030 dx10 = - dx01
00031 dx20 = - dx02
00032 dx21 = - dx12
00033 dx30 = - dx03
00034 dx31 = - dx13
00035 dx32 = - dx23
00036 dx40 = - dx04
00037 dx41 = - dx14
00038 dx42 = - dx24
00039 dx43 = - dx34
00040 xv0 = v - x(ir0)
00041 xv1 = v - x(ir1)
00042 xv2 = v - x(ir2)
00043 xv3 = v - x(ir3)
00044 xv4 = v - x(ir4)
00045
00046 wex0 = (xv2*xv3*xv4 + xv1*xv3*xv4 + &
00047 & xv1*xv2*xv4 + xv1*xv2*xv3)/(dx01*dx02*dx03*dx04)
00048 wex1 = (xv2*xv3*xv4 + xv0*xv3*xv4 + &
00049 & xv0*xv2*xv4 + xv0*xv2*xv3)/(dx10*dx12*dx13*dx14)
00050 wex2 = (xv1*xv3*xv4 + xv0*xv3*xv4 + &
00051 & xv0*xv1*xv4 + xv0*xv1*xv3)/(dx20*dx21*dx23*dx24)
00052 wex3 = (xv1*xv2*xv4 + xv0*xv2*xv4 + &
00053 & xv0*xv1*xv4 + xv0*xv1*xv2)/(dx30*dx31*dx32*dx34)
00054 wex4 = (xv1*xv2*xv3 + xv0*xv2*xv3 + &
00055 & xv0*xv1*xv3 + xv0*xv1*xv2)/(dx40*dx41*dx42*dx43)
00056 dfdx_4th = wex0*f(ir0) + wex1*f(ir1) + wex2*f(ir2) &
00057 & + wex3*f(ir3) + wex4*f(ir4)
00058
00059 end function dfdx_4th
