###  Input data  ###

Spin <- 3           # 3,4,5,6,7; Need to choose a relatively large value compared to the range of integration
steps <- 30         # number of steps
range <- 0.5        # range of integration
divisions <- 360    # for consideration of infinitesimal rotations

###  Construction of fuzzy sphere  ###
## spin-n/2 representation (bosonic cases) ##
s <- Spin
n <- 2*s
Lz <- matrix(0,n+1,n+1)
diag(Lz) <- s:(-s)

m <- (s-1):(-s)
lp <- sqrt((s-m)*(s+m+1))
Ls <- matrix(0,n,n)
diag(Ls) <- lp
Lp <- rbind(Ls, c(rep(0,n)))
Lp <- cbind(c(rep(0,n+1)),Lp)
Lm <- rbind(c(rep(0,n)), Ls)
Lm <- cbind(Lm, c(rep(0,n+1)))

Lx <- 0.5*(Lp + Lm)
Ly <- -0.5i*(Lp - Lm)

# (Lx%*%Lx + Ly%*%Ly + Lz%*%Lz) leads to s(s+1) times identity matrix.
# Now we consider an axial rotation which we can regard as a time-evolution.
# Since the background coordinates should be fixed, 
# the operation is purely on fluctuations A.

I <- matrix(0, n+1, n+1)
diag(I) <- rep(1, n+1)
R <- function(k) {I + (2*k * pi / divisions) * 1i* Lz} 
# representing infinitesimal rotation with k being from 0 to steps

############ need to modify to obtain full S^2 symmetry
## construction of off-diagonal matrices of arbitrary dimensions ##

fx.Ap <- function(a){  
nn <- length(a)
ret <- matrix(0,nn,nn)
diag(ret) <- a
ret <- rbind(ret, rep(0,nn))
ret <- cbind(rep(0,nn+1), ret)
ret
}

fx.Am <- function(a){ -t(fx.Ap(a)) }
fx.Ax <- function(a){ 0.5*(fx.Ap(a) + fx.Am(a)) }
fx.Ay <- function(a){ -0.5i*(fx.Ap(a) - fx.Am(a)) }

##### Construction of action   with length(a) = n = 2s+1 #####

Dx <- function(a,k){ Lx + R(k) %*% fx.Ax(a) }
Dy <- function(a,k){ Ly + R(k) %*% fx.Ay(a) }
Dz <- Lz

### Lagrangian and action to be used
L <- function(a,k){ Dz %*% Dx(a,k) %*% Dy(a,k) - Dz %*% Dy(a,k) %*% Dx(a,k) }  
S <- function(a,k){ -2i* sum(diag(L(a,k))) }
expS <- function(a,k){ exp(1i * S(a,k)) }

fx.AA <- function(a,k){ 0.5*( R(k) %*% fx.Ap(a) %*% R(k) %*% fx.Am(a) + R(k) %*% fx.Am(a) %*% R(k) %*% fx.Ap(a)) }
AA0 <- function(a,k){ diag(fx.AA(a,k))[s+1] * expS(a,k) } 

### evaluated at the equator which is nonzero in general ###
######  to execute multi-dimensional integral  ######

library(adapt)

for( l in 0: (steps) ){
    if(l==0){ expVal <- 0; nor <-0; ave <- 0 }
    ret_norm <- adapt(n, lo=(rep(-range, n)), up=(rep(range, n)), functn =expS,  minpts = 100, maxpts = NULL, eps = 0.15, k=l)
    norm <- ret_norm$value  # to be converged
    ret_AA0 <-  adapt(n, lo=(rep(-range, n)), up=(rep(range, n)), functn =AA0, minpts=100, maxpts = NULL, eps = 0.15, k=l)
    ave_AA0 <- ret_AA0$value
    
    aveAA <- ave_AA0 / norm # this gives trace of aveAA
    expVal <- cbind(expVal, aveAA)
    nor <- cbind(nor, norm)
    ave <- cbind(ave, ave_AA0)
}

expVal <- replace(expVal, c(1), NA)
expVal <- c(expVal[!is.na(expVal)],NA)
nor <- replace(nor, c(1), NA)
nor <- c(nor[!is.na(nor)],NA)
ave <- replace(ave, c(1), NA)
ave <- c(ave[!is.na(ave)],NA)

#expVal #nor #ave

# z <- file("/file/to/be/saved/expVal40_{Spin}.txt", "w")
# cat(expVal, file = z)
# close(z)

plot(expVal)

expVal