From 65bd35b00f8c828653085b71a556d6377dc4dcf8 Mon Sep 17 00:00:00 2001
From: "Dr. Markus Schulze" <markus.schulze@hu-berlin.de>
Date: Fri, 24 Nov 2017 18:25:49 +0100
Subject: [PATCH] Added lightcone checks
---
mod_LoopSpace.f90 | 59 +++++++++++++++++++++++++++++++++++++++--------
mod_Misc.f90 | 15 +++++++++++-
2 files changed, 64 insertions(+), 10 deletions(-)
diff --git a/mod_LoopSpace.f90 b/mod_LoopSpace.f90
index 7d17646..c8db663 100755
--- a/mod_LoopSpace.f90
+++ b/mod_LoopSpace.f90
@@ -35,9 +35,10 @@ real(8) :: sinZ,cosZ,sinT,cosT,sinP,cosP
! urnd(1:4)=(/0.9999999999d0, 0.146446609406726d0, 0.37581005831787d0, 0.422403883500291d0 /)
-
-
! urnd(2) = 0.5d0; print *, "fixing zeta to zero"
+! urnd(3) = 0.0000001d0; print *, "fixing theta to zero"
+! urnd(3) = 0.9999999d0; print *, "fixing theta to one"
+
kE = dsqrt( mu_Loop1_sq*dtan( 0.5d0*Pi*uRnd(1) ) )
kE_sq=kE**2
@@ -52,7 +53,7 @@ real(8) :: sinZ,cosZ,sinT,cosT,sinP,cosP
kLoop1(2) = kE * sinZ * sinT * sinP
kLoop1(3) = kE * sinZ * sinT * cosP
kLoop1(4) = kE * sinZ * cosT
-
+
Jacobian = 2d0*Pi**2*kE_sq/mu_Loop1_sq * ( kE_sq**2+mu_Loop1_sq**2 )*sinZ
!note: for zeta=1/4,3/4Pi we have sinZ^2=cosZ^2 such that kLoop1.kLoop1 = kE^2 ( sinZ^2 - cosZ^2 ) = 0, without constraint on kE
@@ -561,7 +562,7 @@ integer :: j
kappa_Ext(1) = c_pl * k_pl(1) + c_mi * k_mi(1)
kappa_Ext(2:4) =-c_pl * k_pl(2:4) - c_mi * k_mi(2:4)
-
+
! print *, "kTilde",kTilde
! print *, "absvec kTilde",dsqrt(kTilde(2)**2+kTilde(3)**2+kTilde(4)**2),1000d0*0.01d0
! print *, "kappa_Ext",kappa_Ext
@@ -569,8 +570,8 @@ integer :: j
! print *, "ki(1:4,2).dot.kappa_Ext",ki(1:4,2).dot.kappa_Ext(1:4),(ki(1:4,2).dot.ki(1:4,2))-msqi(2)
- kappa_Int(:)=0d0 ! removing kappa_int for now...
-! compute kappa_Int
+! kappa_Int(:)=0d0 ! removing kappa_int for now...
+! !compute kappa_Int
! vij(:,ij12) = 0.5d0*( qi(:,1)+qi(:,2) - (dsqrt(msqi(1))-dsqrt(msqi(2)))*(qi(:,1)-qi(:,2))/Get_MInv(qi(1:4,1)-qi(1:4,2)) )
! vij(:,ij13) = 0.5d0*( qi(:,1)+qi(:,3) - (dsqrt(msqi(1))-dsqrt(msqi(3)))*(qi(:,1)-qi(:,3))/Get_MInv(qi(1:4,1)-qi(1:4,3)) )
! vij(:,ij23) = 0.5d0*( qi(:,2)+qi(:,3) - (dsqrt(msqi(2))-dsqrt(msqi(3)))*(qi(:,2)-qi(:,3))/Get_MInv(qi(1:4,2)-qi(1:4,3)) )
@@ -632,6 +633,46 @@ integer :: j
! pause
+
+
+! if( ((kTilde).dot.(Get_kappa) ).lt.0d0 .and. CloseToLightcone(qi(1:4,1),kTilde,msqi(1)) ) then ! if wrong deformation and close to lightcone of q1
+!
+! kLoop1Full(1:4) = kLoop1(1:4) + (0d0,1d0)*Get_kappa(1:4)
+! print *, "-----"
+! print *, "close to lightcone q1"
+! print *, "D1tilde",Get_MInv2(kTilde-qi(1:4,1))-msqi(1)
+! print *, ""
+! print *, "inside lightcone q1",InsideLightcone(qi(1:4,1),kTilde,msqi(1))
+! print *, "kappa.ktilde",(kTilde).dot.(Get_kappa)
+! print *, ""
+! print *, "-----"
+! pause
+!
+! endif
+
+
+
+! if( IsInInterior(kTilde) ) then
+!
+! print *, "in interior"
+! call printMom( kTilde(1:4) )
+! pause
+! else
+! print *, "in exterior"
+!
+! call printMom( kTilde(1:4) )
+!
+! endif
+
+
+! if( kTilde(4).lt.-5d0 ) then
+! print *, "IsInInterior",IsInInterior(kTilde)
+! call printMom( kTilde(1:4) )
+! pause
+! endif
+
+
+
RETURN
END FUNCTION
@@ -640,10 +681,10 @@ END FUNCTION
-FUNCTION IsInInterior(kTilde)
+FUNCTION IsInInterior(k)
use ModMisc
implicit none
-real(8) :: kTilde(1:4),P_pl(1:4),P_mi(1:4)
+real(8) :: k(1:4),P_pl(1:4),P_mi(1:4)
logical :: IsInInterior
IsInInterior=.false.
@@ -651,7 +692,7 @@ logical :: IsInInterior
P_pl = get_Zpl(qi(:,1),+qi(:,1))
P_mi = get_Zmi(qi(:,2),-qi(:,2))
- if( InsideFWDLightcone(P_pl,kTilde) .and. InsideBWDLightcone(P_mi,kTilde) ) IsInInterior=.true.
+ if( InsideFWDLightcone(P_pl,k) .and. InsideBWDLightcone(P_mi,k) ) IsInInterior=.true.
RETURN
END FUNCTION
diff --git a/mod_Misc.f90 b/mod_Misc.f90
index 40701ee..5a38f3b 100755
--- a/mod_Misc.f90
+++ b/mod_Misc.f90
@@ -235,7 +235,7 @@ real(8) :: RefMom(1:4),TestMom(1:4)
real(8) :: Qsq,RefScale
real(8),optional :: Msq
! *******************************************
-real(8),parameter :: RelSpaceing = 1d-3 !*
+real(8),parameter :: RelSpaceing = 1d-2 !*
! *******************************************
@@ -1419,6 +1419,19 @@ END SUBROUTINE
+
+SUBROUTINE printMom(p)
+implicit none
+real(8) :: p(1:4)
+
+ write(*,"(A,1F22.16,A,1F22.16,A,1F22.16,A,1F22.16)") "E=",p(1)," pX=",p(2)," pY=",p(3)," pZ=",p(4)
+
+
+
+END SUBROUTINE
+
+
+
SUBROUTINE SwitchEnergyComponent(p)
implicit none
real(8) :: p(:,:),px,py,pz,E
--
GitLab