00001 include '../Module/phys_constant.f90'
00002 include '../Module/def_matter.f90'
00003 include '../Module/def_matter_parameter.f90'
00004 include '../Module/def_quantities.f90'
00005 include '../Module/grid_parameter.f90'
00006 include '../Module/coordinate_grav_r.f90'
00007 include '../Module/coordinate_grav_phi.f90'
00008 include '../Module/coordinate_grav_theta.f90'
00009 include '../Module/coordinate_grav_extended.f90'
00010 include '../Module/make_array_1d.f90'
00011 include '../Module/make_array_2d.f90'
00012 include '../Module/make_array_3d.f90'
00013 include '../Module/make_array_4d.f90'
00014 include '../Module/make_array_5d.f90'
00015 include '../Module/trigonometry_grav_theta.f90'
00016 include '../Module/trigonometry_grav_phi.f90'
00017
00018 include '../Module_interface/interface_grgrad_4th_gridpoint_bhex.f90'
00019 include '../Subroutine/grgrad_4th_gridpoint_bhex.f90'
00020 include '../Module_interface/interface_grgrad_4th_gridpoint.f90'
00021 include '../Subroutine/grgrad_4th_gridpoint.f90'
00022
00023 include '.././Subroutine/grgrad_4th.f90'
00024 include '.././Function/dfdx_4th.f90'
00025
00026
00027 program main
00028 implicit none
00029 call coordinate_patch_kit_grav
00030 call fncdiff
00031 end program main
00032 subroutine coordinate_patch_kit_grav
00033 use grid_parameter
00034 use coordinate_grav_r
00035 use coordinate_grav_phi
00036 use coordinate_grav_theta
00037 use coordinate_grav_extended
00038 use trigonometry_grav_theta
00039 use trigonometry_grav_phi
00040 implicit none
00041
00042 call read_parameter
00043 call grid_r
00044 call grid_theta
00045 call trig_grav_theta
00046 call grid_phi
00047 call allocate_trig_grav_mphi
00048 call trig_grav_phi
00049 call grid_extended
00050 end subroutine coordinate_patch_kit_grav
00051
00052 subroutine fncdiff
00053 use grid_parameter
00054 use coordinate_grav_r
00055 use coordinate_grav_phi
00056 use coordinate_grav_theta
00057 use coordinate_grav_extended
00058 use trigonometry_grav_theta
00059 use trigonometry_grav_phi
00060 use make_array_3d
00061 use interface_grgrad_4th_gridpoint_bhex
00062 implicit none
00063 real(8), pointer :: fnc(:,:,:)
00064 real(8) :: x, y, z, dfdx, dfdy, dfdz
00065 real(8) :: df1, df2, df3, e1, e2, e3, error, small = 1.0d-14
00066 integer :: irg, itg, ipg
00067
00068 call alloc_array3d(fnc,0,nrg,0,ntg,0,npg)
00069 do irg = 0, 30
00070 do itg = 0, ntg
00071 do ipg = 0, npg
00072 x = rg(irg)*sinthg(itg)*cosphig(ipg)
00073 y = rg(irg)*sinthg(itg)*sinphig(ipg)
00074 z = rg(irg)*costhg(itg)
00075 fnc(irg,itg,ipg) = ((x**4 + 1)*y + y**3)*dcos(z)
00076 end do
00077 end do
00078 end do
00079
00080 do ipg = 0, npg
00081 do itg = 0, ntg
00082 do irg = 0, 5
00083
00084 x = rg(irg)*sinthg(itg)*cosphig(ipg)
00085 y = rg(irg)*sinthg(itg)*sinphig(ipg)
00086 z = rg(irg)*costhg(itg)
00087
00088 call grgrad_4th_gridpoint_bhex(fnc,dfdx,dfdy,dfdz,irg,itg,ipg)
00089
00090 write(6,*) irg,itg,ipg
00091 write(6,20) x,y,z
00092 write(6,20) dfdx,dfdy,dfdz
00093
00094 df1 = 4*x**3*y*Cos(z)
00095 df2 = (1 + x**4 + 3*y**2)*Cos(z)
00096 df3 = -(((1 + x**4)*y + y**3)*Sin(z))
00097
00098 e1 = 0.0d0
00099 e2 = 0.0d0
00100 e3 = 0.0d0
00101 error = 0.0d0
00102 if (dabs(df1)+dabs(dfdx).ge.small) &
00103 e1 = dabs(df1-dfdx)/(dabs(df1)+dabs(dfdx))
00104 if (dabs(df2)+dabs(dfdy).ge.small) &
00105 e2 = dabs(df2-dfdy)/(dabs(df2)+dabs(dfdz))
00106 if (dabs(df3)+dabs(dfdz).ge.small) &
00107 e3 = dabs(df3-dfdz)/(dabs(df3)+dabs(dfdy))
00108 error = dmax1(e1,e2,e3)
00109 write(6,20) df1,df2,df3
00110 write(6,20) e1,e2,e3,error
00111 end do
00112 end do
00113 end do
00114
00115 20 format(1p,6e12.4)
00116 end subroutine fncdiff
