#
# theoretical analysis program for full-WPT
# Please refer my paper: Hiroshi Saito, Theoretical Analysis of Fully Wireless-Power-Transfer Node Networks, IEICE Trans. Commun., to appear.
#




# nonlinear

func_A <- function(r,k,beta,g,P_t,alpha)
{	work <- k*beta*g*P_t*(r^(-alpha))
	return((work*r)/(1+work))
}



quad_f <- function(q,lambda,P_t,a_v,a_w,s_v,s_w,g,beta,alpha,gamma_s,gamma_0,int_A0,theta_L,theta_H,time_slot)
{
	HtoL<- calculate_HtoL(P_t,a_v,a_w,theta_L,theta_H,time_slot)

	A1 <- exp(-2*pi*lambda*q*int_A0) 
	eta_0 <- gamma_s - gamma_0*A1 #expected amount of energy transferred from other nodes when a node is status 0, r_b=0
	A2 <- calculate_A2(q,beta,g,P_t,alpha,lambda,0) # A2(r_b=0)
	var_u <- calculate_var_u(q,A1,A2,gamma_s,gamma_0)
	var_uu <- var_u 

	LtoH <- calculate_LtoH(eta_0,var_uu,a_v,s_v,a_w,s_w,theta_L,theta_H,time_slot)

	error <- (q*(LtoH + HtoL)-HtoL)^2

	return(error)
}


quad_f0 <- function(q_0,lambda,P_t,a_v,a_w,s_v,s_w,g,beta,alpha,gamma_s,gamma_0,int_Ar_b,theta_L,theta_H,time_slot,q,r_b)
{
	HtoL<- calculate_HtoL(P_t,a_v,a_w,theta_L,theta_H,time_slot)

	A1 <- exp(-2*pi*lambda*q*int_Ar_b) #A_1(r_b)
	eta_r_b <- gamma_s - gamma_0*A1 #expected amount of energy transferred from other nodes when a node is status 0
	A2 <- calculate_A2(q,beta,g,P_t,alpha,lambda,r_b) #A_2(r_b)
	var_u0 <- calculate_var_u(q_0,A1,A2,gamma_s,gamma_0)
	var_uu0 <- var_u0 

	LtoH <- calculate_LtoH(eta_r_b,var_uu0,a_v,s_v,a_w,s_w,theta_L,theta_H,time_slot)

	error <- (q_0*(LtoH + HtoL)-HtoL)^2

	return(error)
}


calculate_A2 <- function(q,beta,g,P_t,alpha,lambda,r_b)
{	work <- integrate(func_A,r_b,Inf,k=2,beta=beta,g=g,P_t=P_t,alpha=alpha)
	int_A02 <- work$value
	A2 <- exp(-2*pi*lambda*q*int_A02) #A_2(r_b)
	return(A2)
}

calculate_var_u <- function(q,A1,A2,gamma_s,gamma_0)
{	var_u <- gamma_0^2*(A2-A1^2)
	return(var_u)
}



calculate_LtoH <- function(eta,var_uu,a_v,s_v,a_w,s_w,theta_L,theta_H,time_slot)
{	ave_total<- eta+a_v-a_w
	var_total <- var_uu+s_v^2+s_w^2
	work1 <- exp(-2*ave_total*theta_L/(time_slot*var_total))
	work2 <- exp(-2*ave_total*theta_H/(time_slot*var_total))
	fst_theory <- (theta_H-theta_L)/(time_slot*ave_total)
	fst_theory <- fst_theory-(var_total/(2*ave_total^2))*(work1-work2)
	return(fst_theory)
}

calculate_HtoL <- function(P_t,a_v,a_w,theta_L,theta_H,time_slot)
{	result <- (theta_H-theta_L)/(time_slot*(P_t-a_v+a_w))
	return(result)
}


# parameters

time_slot <- 0.01 #[s]
lambda <- 0.1
alpha <- 4
g <- 1
P_t <- 50 #[mW]
a_v <- 10 #[mW]
a_w <- 10 #[mW]
s_v <- 100 #[mW]
s_w <- 100 #[mW]
av_aw <- a_v-a_w #[mW]
gamma_s <- 11.05
gamma_0 <- 11.05
beta <- 0.1
theta_H <- 100 #[mW*slot]
theta_L <- 50 # [mW*slot]




#
# main part
#

work <- integrate(func_A,0,Inf,k=1,beta=beta,g=g,P_t=P_t,alpha=alpha)
int_A0 <- work$value # A_k, k=1, r_b=0

opt <- optimize(quad_f,interval = c(0,0.99999),lambda,P_t,a_v,a_w,s_v,s_w,g,beta,alpha,gamma_s,gamma_0,int_A0,theta_L,theta_H,time_slot) # solve nonlinear equation and obtain q
q <-opt$minimum # prob(status = 1)

HtoL <- calculate_HtoL(P_t,a_v,a_w,theta_L,theta_H,time_slot) # FST(LtoH) of non-typical node
A1 <- exp(-2*pi*lambda*q*int_A0)
eta_0 <- gamma_s - gamma_0*A1 #expected amount of energy transferred from other nodes, r_b=0
ave_u <- eta_0
A2 <- calculate_A2(q,beta,g,P_t,alpha,lambda,0) #A2 with b_r=0
var_u <- calculate_var_u(q,A1,A2,gamma_s,gamma_0) #r_b=0
var_uu <- var_u 
LtoH <- calculate_LtoH(eta_0,var_uu,a_v,s_v,a_w,s_w,theta_L,theta_H,time_slot) # FST(HtoL) of non-typical node



q_0_v <- NULL; eta_r_b_v <- NULL; LJR <- NULL; lost <- NULL; consume <- NULL
r_b_v <- c(0,1,2,3,4,5,10) # first element must be 0.

for(r_b in r_b_v)
{	work <- integrate(func_A,r_b,Inf,k=1,beta=beta,g=g,P_t=P_t,alpha=alpha)
	int_Ar_b <- work$value #lower-bound of intg is r_b
	A1r_b <- exp(-2*pi*lambda*q*int_Ar_b) #A_1(r_b)
	eta_r_b<- gamma_s - gamma_0*A1r_b #expected amount of energy transferred from other nodes to a typical node when the node is status 0 and no other nodes within r_b
	eta_r_b_v<- c(eta_r_b_v,eta_r_b)
	opt <- optimize(quad_f0,interval = c(0,0.99999),lambda,P_t,a_v,a_w,s_v,s_w,g,beta,alpha,gamma_s,gamma_0,int_Ar_b,theta_L,theta_H,time_slot,q,r_b) # solve nonlinear equation and obtain q0
	q_0 <- opt$minimum # prob(status of typical node = 1)
	q_0_v <- c(q_0_v,q_0)
	ave_u0 <- eta_r_b*(1-q_0)
	LJR <- c(LJR, 1-(eta_r_b*(1-q_0)+a_v)/(P_t*q_0+a_w)) # lost job ratio
	lost <- c(lost, P_t*q_0+a_w-a_v-ave_u0)
	consume <- c(consume,(P_t*q_0+a_w))
}	

#
# output
#


path <- "C:\\Users\\Hiroshi Saito\\Documents\\papers\\fullWPT\\"
file_name <- "theory_public_out.txt"
out_file <- paste(path,file_name,sep="")
parameters <- c("node density","alpha","gamma_sat","gamma_0","beta","gain", "Tx-power","ave_v","ave_w","sd_v","sd_w","theta_H","theta_L","time-slot")
values <- c(lambda,alpha,gamma_s,gamma_0,beta,g,P_t,a_v,a_w,s_v,s_w,theta_H,theta_L,time_slot)
write.table(data.frame(parameters,values),out_file,append=T,quote=F,col.names=F,row.names=F)
write("statistics of non-typical nodes",out_file,append=T)
description <- c("prob(status=1):","first passge L to H:","first passge H to L:", "mean energy at status 0:","var energy at status 0:")
values <- c(q,LtoH,HtoL,ave_u,var_u)
write.table(data.frame(description,values),out_file,append=T,quote=F,col.names=F,row.names=F)
write("statistics of typical node",out_file,append=T)
description <- c("r_b:","lost job ratio:", "prob(status of node_0=1):","energy to node_0:", "energy from 0:", "lost:", "consume")
df <- data.frame(r_b_v,LJR,q_0_v,eta_r_b_v*(1-q_0_v),P_t*q_0_v,lost,consume)
colnames(df) <- description
write.table(df,out_file,append=T,quote=F,col.names=T,row.names=F)
write("\n\n",out_file,append=T)