00001 subroutine calc_redblue_shift_WL
00002 use phys_constant, only : long, pi
00003 use coordinate_grav_r, only : rg
00004 use grid_parameter, only : nrf, ntf, npf, &
00005 & ntfpolp, ntfeq, ntfxy, npfxzp, npfyzp
00006 use trigonometry_grav_theta, only : sinthg, costhg
00007 use trigonometry_grav_phi, only : sinphig, cosphig
00008 use def_metric, only : psi, alph, bvxd, bvyd, bvzd, bvxu, bvyu, bvzu
00009 use def_metric_hij, only : hxxd, hxyd, hxzd, hyyd, hyzd, hzzd
00010 use def_matter, only : rs
00011 use def_matter_parameter, only : ome, ber, radi
00012 use def_quantities, only : zrb_xp_plus, zrb_xp_minus, &
00013 & zrb_yp_plus, zrb_yp_minus, &
00014 & zrb_zp_plus, zrb_zp_minus
00015 use interface_interpo_radial1p_grav
00016 implicit none
00017 integer :: ir, it, ip, ii
00018 real(long) :: rv, xx, yy, zz
00019 real(long) :: vphixd, vphiyd, vphizd, vphixu, vphiyu, vphizu
00020 real(long) :: psiw, alphw
00021 real(long) :: bvxdw, bvydw, bvzdw, bvxuw, bvyuw, bvzuw
00022 real(long) :: hxxdw, hxydw, hxzdw, hyydw, hyzdw, hzzdw
00023 real(long) :: gtt, gtp, gpp, hhgc, utgc, bplus, bminus
00024 real(long) :: kdott_plus, kdotu_plus, kdott_minus, kdotu_minus
00025
00026
00027
00028
00029 do ii = 1, 2
00030 if (ii.eq.1) then; it = ntfeq; ip = npfxzp; end if
00031 if (ii.eq.2) then; it = ntfeq; ip = npfyzp; end if
00032
00033 rv = rg(nrf)*rs(it,ip)
00034 xx = rv*sinthg(it)*cosphig(ip)
00035 yy = rv*sinthg(it)*sinphig(ip)
00036 zz = rv*costhg(it)
00037 vphixu = - yy
00038 vphiyu = xx
00039 vphizu = 0.0d0
00040 hhgc = 1.0d0
00041 utgc = hhgc/ber
00042 call interpo_radial1p_grav(psi,psiw,rv,it,ip)
00043 call interpo_radial1p_grav(alph,alphw,rv,it,ip)
00044
00045 call interpo_radial1p_grav(bvxd,bvxdw,rv,it,ip)
00046 call interpo_radial1p_grav(bvyd,bvydw,rv,it,ip)
00047 call interpo_radial1p_grav(bvzd,bvzdw,rv,it,ip)
00048 call interpo_radial1p_grav(bvxu,bvxuw,rv,it,ip)
00049 call interpo_radial1p_grav(bvyu,bvyuw,rv,it,ip)
00050 call interpo_radial1p_grav(bvzu,bvzuw,rv,it,ip)
00051
00052 call interpo_radial1p_grav(hxxd,hxxdw,rv,it,ip)
00053 call interpo_radial1p_grav(hxyd,hxydw,rv,it,ip)
00054 call interpo_radial1p_grav(hxzd,hxzdw,rv,it,ip)
00055 call interpo_radial1p_grav(hyyd,hyydw,rv,it,ip)
00056 call interpo_radial1p_grav(hyzd,hyzdw,rv,it,ip)
00057 call interpo_radial1p_grav(hzzd,hzzdw,rv,it,ip)
00058
00059 vphixd = (1.0d0+hxxdw)*vphixu + hxydw*vphiyu + hxzdw*vphizu
00060 vphiyd = hxydw*vphixu + (1.0d0+hyydw)*vphiyu + hyzdw*vphizu
00061 vphizd = hxzdw*vphixu + hyzdw*vphiyu + (1.0d0+hzzdw)*vphizu
00062
00063 gtt = - alphw**2 + psiw**4*(bvxdw*bvxuw + bvydw*bvyuw + bvzdw*bvzuw)
00064 gtp = psiw**4*(bvxdw*vphixu + bvydw*vphiyu + bvzdw*vphizu)
00065 gpp = psiw**4*(vphixd*vphixu + vphiyd*vphiyu + vphizd*vphizu)
00066 bplus = (- gtp + sqrt(gtp**2 - gtt*gpp))/gpp
00067 bminus = (- gtp - sqrt(gtp**2 - gtt*gpp))/gpp
00068
00069 kdott_plus = gtt + bplus*gtp
00070 kdotu_plus = utgc*(gtt + (bplus + ome)*gtp + bplus*ome*gpp)
00071 kdott_minus = gtt + bminus*gtp
00072 kdotu_minus = utgc*(gtt + (bminus + ome)*gtp + bminus*ome*gpp)
00073
00074 if (ii.eq.1) then
00075 zrb_xp_plus = kdotu_plus/kdott_plus - 1.0d0
00076 zrb_xp_minus = kdotu_minus/kdott_minus - 1.0d0
00077 end if
00078 if (ii.eq.2) then
00079 zrb_yp_plus = kdotu_plus/kdott_plus - 1.0d0
00080 zrb_yp_minus = kdotu_minus/kdott_minus - 1.0d0
00081 end if
00082
00083 end do
00084
00085 hhgc = 1.0d0
00086 utgc = hhgc/ber
00087 zrb_zp_plus = utgc - 1.0d0
00088 zrb_zp_minus = utgc - 1.0d0
00089
00090 end subroutine calc_redblue_shift_WL
