## omega_lm
## this needs qtilde as input!
omegalm <- function(l=0, m=0, q, gamma=1, dvec=c(0,0,0)) {
return( LuescherZeta(qsq = q^2, l=l, m=m, gamma=gamma, dvec=dvec)/(pi^(3/2)*sqrt(2*l+1)*q^(l+1)*gamma) )
}
## center of mass energy from q^2 and Mpi
Ecmofqsq <- function(q, Mpi) {
return(2*acosh(cosh(Mpi) + 2*sin(q/2)^2 ))
}
Ecmofqsq.cont <- function(q, Mpi) {
return(2*sqrt(Mpi^2 + q^2))
}
## energy level in MF from Ecm
EofEcm <- function(Ecm, dvec, L) {
return(acosh( cosh(Ecm) + 2*sum(sin(pi*dvec/L)^2) ))
}
EofEcm.cont <- function(Ecm, dvec, L) {
return(sqrt(Ecm^2 + sum((2*pi*dvec/L)^2)))
}
## computes the boost factor as a function of q^2
gammaofq <- function(q, dvec, Mpi, L) {
if(lat.disp) {
Ecm <- Ecmofqsq(q, Mpi)
return(EofEcm(Ecm, dvec, L)/Ecm)
}
Ecm <- Ecmofqsq.cont(q, Mpi)
## E/Ecm
return(EofEcm.cont(Ecm, dvec, L)/Ecm)
}
## computes the energy level given q^2
Eofqsq <- function(q, dvec, Mpi, L) {
if(lat.disp) {
Ecm <- Ecmofqsq(q, Mpi)
return(EofEcm(Ecm, dvec, L))
}
Ecm <- Ecmofqsq.cont(q, Mpi)
return(EofEcm.cont(Ecm, dvec, L))
}
## center of mass frame formulae
## this function is indendet for a rough scan in q
## and it is independent of the irrep
prepdetEqCMscan <- function(q, Mpi, L) {
W <- list()
W$dvec <- c(0,0,0)
W$gamma <- gammaofq(q, dvec=W$dvec, Mpi=Mpi, L=L)
W$qt <- L*q/2/pi
W$L <- L
W$Mpi <- Mpi
W$q <- q
W$w00 <- Re(omegalm(l=0, m=0, q=W$qt, gamma=W$gamma))
W$w40 <- Re(omegalm(l=4, m=0, q=W$qt, gamma=W$gamma))
return(W)
}
## A1 irrep
detEqCMA1scan <- function(par, W) {
cd0 <- par[1] + 0.5*par[2]*W$q^2
return(cd0 - W$q*W$w00)
}
detEqCMA1 <- function(q, par, Mpi, L) {
gamma <- rep(1., times=(length(q)))
qt <- L*q/2/pi
cd0 <- par[1] + 0.5*par[2]*q^2
return(cd0 - q*Re(omegalm(l=0, m=0, q=qt, gamma=gamma)))
}
## E irrep
detEqCMEscan <- function(par, W) {
cd2 <- par[3]
return(cd2 - W$q^5*(W$w00 + 18/7*W$w40))
}
detEqCME <- function(q, par, Mpi, L) {
gamma <- rep(1., times=(length(q)))
qt <- L*q/2/pi
w00 <- Re(omegalm(l=0, m=0, q=qt, gamma=gamma))
w40 <- Re(omegalm(l=4, m=0, q=qt, gamma=gamma))
cd2 <- par[3]
return(cd2 - q^5*(w00 + 18/7*w40))
}
## T2 irrep
detEqCMT2scan <- function(par, W) {
cd2 <- par[3]
return(cd2 - W$q^5*(W$w00 - 12/7*W$w40))
}
detEqCMT2 <- function(q, par, Mpi, L) {
gamma <- rep(1., times=(length(q)))
qt <- L*q/2/pi
w00 <- Re(omegalm(l=0, m=0, q=qt, gamma=gamma))
w40 <- Re(omegalm(l=4, m=0, q=qt, gamma=gamma))
cd2 <- par[3]
return(cd2 - q^5*(w00 - 12/7*w40))
}
## moving frame with d=(0,0,1)
prepdetEqMF1scan <- function(q, Mpi, L) {
W <- list()
W$dvec <- c(0,0,1)
W$gamma <- gammaofq(q, dvec=W$dvec, Mpi=Mpi, L=L)
W$qt <- L*q/2/pi
W$L <- L
W$Mpi <- Mpi
W$q <- q
W$w00 <- Re(omegalm(l=0, m=0, q=W$qt, gamma=W$gamma, dvec=W$dvec))
W$w20 <- Re(omegalm(l=2, m=0, q=W$qt, gamma=W$gamma, dvec=W$dvec))
W$w40 <- Re(omegalm(l=4, m=0, q=W$qt, gamma=W$gamma, dvec=W$dvec))
return(W)
}
## A1 irrep
detEqMF1A1scan <- function(par, W) {
cd0 <- par[1] + 0.5*par[2]*W$q^2
if(length(par) < 3) {
return(cd0 - W$q*W$w00)
}
cd2 <- par[3]
return( (cd0 - W$q*W$w00)*
(cd2 - W$q^5*(W$w00 + 10/7*W$w20 + 18/7*W$w40))
-W$q^6*5*W$w20^2
)
}
detEqMF1A1 <- function(q, par, Mpi, L) {
dvec <- c(0,0,1)
gamma <- gammaofq(q, dvec=dvec, Mpi=Mpi, L=L)
qt <- L*q/2/pi
w00 <- Re(omegalm(l=0, m=0, q=qt, gamma=gamma, dvec=dvec))
cd0 <- par[1] + 0.5*par[2]*q^2
if(length(par)<3) {
return(cd0 - q*w00)
}
w20 <- Re(omegalm(l=2, m=0, q=qt, gamma=gamma, dvec=dvec))
w40 <- Re(omegalm(l=4, m=0, q=qt, gamma=gamma, dvec=dvec))
cd2 <- par[3]
return( (cd0 - q*w00)*
(cd2 - q^5*(w00 + 10/7*w20 + 18/7*w40))
-q^6*5*w20^2
)
}
## moving frame with d=(1,1,0)
prepdetEqMF2scan <- function(q, Mpi, L) {
W <- list()
W$dvec <- c(1,1,0)
W$gamma <- gammaofq(q, dvec=W$dvec, Mpi=Mpi, L=L)
W$qt <- L*q/2/pi
W$L <- L
W$Mpi <- Mpi
W$q <- q
W$w00 <- Re(omegalm(l=0, m=0, q=W$qt, gamma=W$gamma, dvec=W$dvec))
W$w20 <- Re(omegalm(l=2, m=0, q=W$qt, gamma=W$gamma, dvec=W$dvec))
W$w40 <- Re(omegalm(l=4, m=0, q=W$qt, gamma=W$gamma, dvec=W$dvec))
W$w22 <- omegalm(l=2, m=2, q=W$qt, gamma=W$gamma, dvec=W$dvec)
W$w44 <- Re(omegalm(l=4, m=4, q=W$qt, gamma=W$gamma, dvec=W$dvec))
W$w42 <- omegalm(l=4, m=2, q=W$qt, gamma=W$gamma, dvec=W$dvec)
return(W)
}
## A1 irrep
detEqMF2A1scan <- function(par, W) {
cd0 <- (par[1]/W$q + 0.5*par[2]*W$q)
if(length(par) < 3 ) {
return(cd0 - W$q*W$w00)
}
cd2 <- par[3]/W$q^5
return(Re(W$q^11*(10/49*(10*(-2*cd0 + 2*W$w00 + 7*W$w20)*W$w22^2 +
3*sqrt(15)*(4*cd0 - 4*W$w00 - 7*W$w20)*W$w22*W$w42 + 27*(W$w00 - cd0)*W$w42^2)
+(-cd2 + W$w00 + 2/7*(5*W$w20 + 9*W$w40))*(10*W$w22^2 + 1/7*(cd0 - W$w00)*(7*cd2 - 7*W$w00 + 10*W$w20 - 3*W$w40 + 3*sqrt(70)*W$w44)) +
5/7*W$w20*(20*W$w22^2 - 6*sqrt(15)*W$w22*W$w42 + W$w20*(7*cd2 - 7*W$w00 + 10*W$w20 - 3*W$w40 + 3*sqrt(70)*W$w44))
)
))
}
detEqMF2A1 <- function(q, par, Mpi, L) {
dvec=c(1,1,0)
gamma <- gammaofq(q, dvec=dvec, Mpi=Mpi, L=L)
qt <- L*q/2/pi
w00 <- Re(omegalm(l=0, m=0, q=qt, gamma=gamma, dvec=dvec))
cd0 <- (par[1]/q + 0.5*par[2]*q)
if(length(par) < 3) {
return(cd0 - q*w00)
}
w20 <- Re(omegalm(l=2, m=0, q=qt, gamma=gamma, dvec=dvec))
w22 <- omegalm(l=2, m=2, q=qt, gamma=gamma, dvec=dvec)
w40 <- Re(omegalm(l=4, m=0, q=qt, gamma=gamma, dvec=dvec))
w44 <- Re(omegalm(l=4, m=4, q=qt, gamma=gamma, dvec=dvec))
w42 <- omegalm(l=4, m=2, q=qt, gamma=gamma, dvec=dvec)
cd2 <- par[3]/q^5
return(Re(q^11*(10/49*(10*(-2*cd0 + 2*w00 + 7*w20)*w22^2 +
3*sqrt(15)*(4*cd0 - 4*w00 - 7*w20)*w22*w42 + 27*(w00 - cd0)*w42^2)
+(-cd2 + w00 + 2/7*(5*w20 + 9*w40))*(10*w22^2 + 1/7*(cd0 - w00)*(7*cd2 - 7*w00 + 10*w20 - 3*w40 + 3*sqrt(70)*w44)) +
5/7*w20*(20*w22^2 - 6*sqrt(15)*w22*w42 + w20*(7*cd2 - 7*w00 + 10*w20 - 3*w40 + 3*sqrt(70)*w44))
)
))
}
## moving frame with d=(1,1,1)
prepdetEqMF3scan <- function(q, Mpi, L) {
W <- list()
W$dvec <- c(1,1,1)
W$gamma <- gammaofq(q, dvec=W$dvec, Mpi=Mpi, L=L)
W$qt <- L*q/2/pi
W$L <- L
W$Mpi <- Mpi
W$q <- q
W$w00 <- Re(omegalm(l=0, m=0, q=W$qt, gamma=W$gamma, dvec=W$dvec))
W$w40 <- Re(omegalm(l=4, m=0, q=W$qt, gamma=W$gamma, dvec=W$dvec))
W$w22 <- omegalm(l=2, m=2, q=W$qt, gamma=W$gamma, dvec=W$dvec)
W$w42 <- omegalm(l=4, m=2, q=W$qt, gamma=W$gamma, dvec=W$dvec)
return(W)
}
## A1 irrep
detEqMF3A1scan <- function(par, W) {
cd0 <- (par[1] + 0.5*par[2]*W$q^2)
if(length(par) < 3) {
return(cd0 - W$q*W$w00)
}
cd2 <- par[3]
return(Re(( cd0 - W$q*W$w00 )*
( cd2 - W$q^5*(W$w00 - 12/7*W$w40 - 12*sqrt(10)/7*1i*W$w42 - 10*sqrt(6)/7*1i*W$w22)) +
W$q^6*30*W$w22^2
)
)
}
detEqMF3A1 <- function(q, par, Mpi, L) {
dvec=c(1,1,1)
gamma <- gammaofq(q, dvec=dvec, Mpi=Mpi, L=L)
qt <- L*q/2/pi
w00 <- Re(omegalm(l=0, m=0, q=qt, gamma=gamma, dvec=dvec))
cd0 <- (par[1] + 0.5*par[2]*q^2)
if(length(par) < 3) {
return(cd0 - q*w00)
}
w40 <- Re(omegalm(l=4, m=0, q=qt, gamma=gamma, dvec=dvec))
w22 <- omegalm(l=2, m=2, q=qt, gamma=gamma, dvec=dvec)
w42 <- omegalm(l=4, m=2, q=qt, gamma=gamma, dvec=dvec)
cd2 <- par[3]
return(Re(( cd0 - q*w00 )*
( cd2 - q^5*(w00 - 12/7*w40 - 12*sqrt(10)/7*1i*w42 - 10*sqrt(6)/7*1i*w22)) +
q^6*30*w22^2
)
)
}
findSignChanges <- function(fn, makeplot=FALSE, par, W, no=3, threshold=100) {
res <- fn(par=par, W=W)
if(makeplot) {
if(interactive()) X11()
plot(W$q^2, res, type="l", ylim=c(-threshold, threshold))
abline(a=0, b=0, col="red")
}
ii <- integer(0)
j <- 1
for(i in c(2:(length(W$q)-1))) {
## find all sign-changes
##cat(res[i], res[i+1], "\n")
if(!is.na(res[i]) && !is.na(res[i+1])) {
if(res[i]/abs(res[i]) != res[i+1]/abs(res[i+1])) {
## remove the poles
if((res[i] < 0 && res[i-1] < res[i]) ||
(res[i] > 0 && res[i-1] > res[i])) {
ii[j] <- i
if(j == no) break
j <- j+1
}
}
}
}
return(ii)
}
findZeros <- function(fn, q, ii, makeplot=FALSE, tol=1.e-14, ...) {
## now determine the root more precisely
zeros <- numeric(0.)
for(j in c(1:length(ii))) {
z <- uniroot(fn, interval=c(q[ii[j]], q[ii[j]+1]), tol = tol, ...)
if(z$f.root < 1) {
zeros[j] <- z$root
##cat("root", j, ":", z$root, z$f.root, z$estim.prec, "\n")
}
}
if(makeplot) {
abline(v=zeros^2, col="blue")
}
return(zeros)
}
getEvalues <- function(par, Wlist, Lvalues, nolist, framelist, scanfnlist, fnlist, irreplist) {
evalues <- c()
j <- 1
k <- 1
for(L in Lvalues) {
for(f in framelist) {
for(i in irreplist) {
if(nolist[[k]] > 0) {
ii <- findSignChanges(fn=scanfnlist[[i]][[f]], par=par, makeplot=FALSE, W=Wlist[[j]], no=nolist[[k]])
if(length(ii) != nolist[[k]]) return(NA)
zeros <- findZeros(fn=fnlist[[i]][[f]], q=Wlist[[j]]$q, ii=ii, L=L, Mpi=Wlist[[j]]$Mpi, par=par, makeplot=FALSE)
if(length(zeros) != nolist[[k]]) return(NA)
evalues <- c(evalues, Eofqsq(zeros, dvec=Wlist[[j]]$dvec, Mpi=Wlist[[j]]$Mpi, L=L))
}
k <- k+1
}
j <- j+1
}
}
return(evalues)
}
getEDiff <- function(x, par, index, level, scanfn, dfn, W, E) {
par[index] <- x
ii <- findSignChanges(fn=scanfn, par=par, makeplot=FALSE, W=W, no=level)
if(length(ii) != level) return(NA)
zero <- findZeros(fn=dfn, q=W$q, ii=ii[level], L=W$L, Mpi=W$Mpi, par=par, makeplot=FALSE)
if(length(zero) != 1) return(NA)
return(Eofqsq(zero, dvec=W$dvec, Mpi=W$Mpi, L=W$L)-E)
}
getq <- function(par, level, scanfn, dfn, W, E) {
ii <- findSignChanges(fn=scanfn, par=par, makeplot=FALSE, W=W, no=level)
if(length(ii) != level) return(NA)
zero <- findZeros(fn=dfn, q=W$q, ii=ii[level], L=W$L, Mpi=W$Mpi, par=par, makeplot=FALSE)
if(length(zero) != 1) return(NA)
return(zero)
}
solveforCotDelta <- function(par, index, level, scanfn, fn, W, E, interval) {
n <- 0
if(missing(interval)) {
interval <- c(-pi/2.,-0.1)
if(index == 3) interval <- c(-0.025,-0.001)
}
l <- getEDiff(x=interval[1], par, index=index, level=level, scanfn=scanfn, dfn=fn, W=W, E=E)
r <- getEDiff(x=interval[2], par, index=index, level=level, scanfn=scanfn, dfn=fn, W=W, E=E)
while(is.na(l) || is.na(r) || l/abs(l) == r/abs(r)) {
if(n > 50*index) {
cat("NA\n")
return(NA)
}
if(index == 1) interval[1] <- interval[1]*9/10.
if(index == 3) interval[1] <- interval[1]*1.1
l <- getEDiff(x=interval[1], par, index=index, level=level, scanfn=scanfn, dfn=fn, W=W, E=E)
r <- getEDiff(x=interval[2], par, index=index, level=level, scanfn=scanfn, dfn=fn, W=W, E=E)
n <- n+1
}
qcotdelta <- uniroot(f = getEDiff, interval=interval, tol = 1.e-14, par=par, index=index, level=level, scanfn=scanfn, dfn=fn, W=W, E=E)
cat("level", level, qcotdelta$root, qcotdelta$f.root, qcotdelta$estim.prec, "\n")
return(qcotdelta$root)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.