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