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") # #