require(graphics)
require(ggplot2)
require(rootSolve)

### Testvalues ###
v <- seq(0.34,6,0.01)

### theta-function ###
theta <- function(x) as.numeric(x > 0)

### trapez integration ###
trap <- function(x, y){
    idx = 2:length(x);
    return (as.double( (x[idx] - x[idx - 1]) %*% (y[idx] + y[idx - 1])) / 2)
}

### vWG ###

pvw_tv   <- function(t,v) theta(3*v - 1)*(8*t/(3*v - 1) - 3/v^2)
pvw_sg_v <- function(v) theta(3*v - 2)*((3*v - 2)/v^3)
pvw_mm   <- function(t,v) 4*t*v^3 - 9*v^2 + 6*v - 1

fv1v2    <- function(pm,v1,v2) p = integrate(pvw_tv,v1,v2,t=0.95)$value - pm*(v2 - v1)
ftv1v2   <- function(t,pm,v1,v2) p = integrate(pvw_tv,v1,v2,t=t)$value - pm*(v2 - v1)

### determining upper and lower boundary and value (v_1, v_2, p_m) ###

vw_msl <- function(t) { 

  dv <- 0.001
  vmin <- uniroot(function(v) pvw_mm(t,v),c(0.34,1))$root
  vmax <- uniroot(function(v) pvw_mm(t,v),c(1,10))$root
  pma  <- (pvw_tv(t,vmax) + pvw_tv(t,vmin))/2 
  res  <- 1

  while(res > 0){
  
    pma <- pma + dv
    v1a <- uniroot(function(v) pvw_tv(t,v) - pma, c(0.34,vmin))$root
    v2a <- uniroot(function(v) pvw_tv(t,v) - pma, c(vmax,100))$root
    res <- ftv1v2(t,pma,v1a,v2a)

  } 

 return(c(v1a,v2a,pma))

}


### dataframes of all curves ###
pvw_sg <- pvw_sg_v(v); 
pvw_t085 <- pvw_tv(0.85,v);
pvw_t090 <- pvw_tv(0.90,v);
pvw_t095 <- pvw_tv(0.95,v);
pvw_t1 <- pvw_tv(1,v);
pvw_t105 <- pvw_tv(1.05,v);
pvw_t110 <- pvw_tv(1.10,v);
pvw_t115 <- pvw_tv(1.15,v);
vw_df <- data.frame(v,pvw_sg, pvw_t085, pvw_t090, pvw_t095, pvw_t1, pvw_t105, pvw_t110, pvw_t115)

vw_msl_t095 = vw_msl(0.95)
vw_msl_t090 = vw_msl(0.90)
vw_msl_t085 = vw_msl(0.85)


### plotting all data ###
ggplot(data = vw_df) +
  geom_line(aes(v, pvw_sg, color = "SG"), size = 1.2) + 
  geom_line(aes(v, pvw_t085, color = "t=0.85"), size = 0.7) + 
  geom_line(aes(v, pvw_t090, color = "t=0.90"), size = 0.7) + 
  geom_line(aes(v, pvw_t095, color = "t=0.95"), size = 0.7) + 
  geom_line(aes(v, pvw_t1, color = "t=1"), size = 0.7) + 
  geom_line(aes(v, pvw_t105, color = "t=1.05"), size = 0.7) + 
  geom_line(aes(v, pvw_t110, color = "t=1.10"), size = 0.7) + 
  geom_line(aes(v, pvw_t115, color = "t=1.15"), size = 0.7) + 
  geom_segment(aes(x=vw_msl_t085[1],y=vw_msl_t085[3],xend=vw_msl_t085[2],yend=vw_msl_t085[3]))  + 
  geom_segment(aes(x=vw_msl_t090[1],y=vw_msl_t090[3],xend=vw_msl_t090[2],yend=vw_msl_t090[3]))  + 
  geom_segment(aes(x=vw_msl_t095[1],y=vw_msl_t095[3],xend=vw_msl_t095[2],yend=vw_msl_t095[3]))  + 
  xlim(0,4) + ylim(0,2) +
  labs(x = "v", y = "p", colour = " ") +
  theme(legend.position = c(0.9, 0.75))

ggsave("Phasendiagramm-vWG.pdf", width = 15, height = 15, units = "cm")

#  
#