From 46c22f8ac30a4a4edb98374f5ae9cfe4da85088e Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Tue, 30 Jan 2018 11:06:49 -0500
Subject: [PATCH] Ink bottle seems right

---
 example/InkBottle/InkBottle.R | 48 ++++++++++++++++++++++++++---------
 1 file changed, 36 insertions(+), 12 deletions(-)

diff --git a/example/InkBottle/InkBottle.R b/example/InkBottle/InkBottle.R
index fd95faf9..d1b1c9e7 100644
--- a/example/InkBottle/InkBottle.R
+++ b/example/InkBottle/InkBottle.R
@@ -1,4 +1,6 @@
-require(ggplot2)
+require("ggplot2")
+
+GGTHEME<-theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
 
 PI=3.14159
 
@@ -37,30 +39,52 @@ Rb=1.2
 # Length of the bulb
 Hb=sqrt(Rb^2-r1^2)+sqrt(Rb^2-r2^2)
 
-time=seq(1,100)
+time=seq(0,100)
 time_inlet<-time
 
-M0_inlet<-CapVolume(r1,r1)+CylVolume(r1,time*0.03)
-M1_inlet<-CapArea(r1,r1)+CylArea(r1,time*0.03)
+M0_inlet<-CapVolume(r1,r1)+CylVolume(r1,time*0.01*L1)
+M1_inlet<-CapArea(r1,r1)+CylArea(r1,time*0.01*L1)
 M2_inlet<-2*CapArea(r1,r1)/r1+CylArea(r1,time*0.03)/r1
 Jwn_inlet<-2/r1
 
 # radius of the fluid interface while contact line is pinned
 rt_pore<-r1+0.01*(Rb-r1)*time
 # height of the spherical cap at the inlet side 
-ht_pore<-(-rt_pore+sqrt(rt_pore^2+r1^2))
+ht_pore<-(rt_pore-sqrt(rt_pore^2-r1^2))
 # depth of penetration into the pore
 xt_pore<-2*rt_pore-ht_pore
 
-time_pore<-time+100;
-M0_pore<-CylVolume(r1,L1)+CapVolume(r1,xt_pore)
-M1_pore<-CylArea(r1,L1)+CapArea(rt,xt_pore)
-M2_pore<-2*CapArea(r1,r1)/r1+CylArea(r1,L1)/r1+2*CapVolume(r1,xt_pore)/rt_pore
+time_pore<-time+100
+M0_pore<-CylVolume(r1,L1)+pi*xt_pore*(3*r1^2+xt_pore^2)/6
+M1_pore<-CylArea(r1,L1)+2*pi*rt_pore*xt_pore
+M2_pore<-2*CapArea(r1,r1)/r1+CylArea(r1,L1)/r1+2*pi*xt_pore
 Jwn_pore<-2/rt_pore
 
+rt_pin<-Rb+(r2-Rb)*time*0.01
+ht_pin<-rt_pin-sqrt(rt_pin^2-r2^2)
+M0_pin<-CylVolume(r1,L1)+BulbVolume(r1,r2,Hb)+pi*ht_pin*(3*r2^2+ht_pin^2)/6
+M1_pin<-CylArea(r1,L1)+2*pi*Rb*Hb+2*pi*rt_pin*ht_pin
+M2_pin<-2*CapArea(r1,r1)/r1+CylArea(r1,L1)/r1+2*pi*Hb+2*pi*ht_pin
+Jwn_pin<-2/rt_pin
+
 time_outlet<-time_pore+100;
 Jwn_outlet<-2/r2
-M0_outlet<-CylVolume(r1,L1)+BubVolume(r1,r2,Hb)+CapVolume(r2,r2)+CylVolume(r2,time*0.03)
-M1_outlet<-CylArea(r1,L1)+CapArea(Rb,Hb)+CapArea(r2,r2)+CylArea(r2,time*0.03)
-M2_outlet<-2*CapArea(r2,r2)/r2+CylArea(r2,time*0.03)/r2
+M0_outlet<-CylVolume(r1,L1)+BulbVolume(r1,r2,Hb)+CapVolume(r2,r2)+CylVolume(r2,time*0.01*L2)
+M1_outlet<-CylArea(r1,L1)+CapArea(Rb,Hb)+CapArea(r2,r2)+CylArea(r2,time*0.01*L2)
+M2_outlet<-2*CapArea(r2,r2)/r2+CylArea(r2,time*0.01*L2)/r2
+
+M0<-as.vector(rbind(M0_inlet,M0_pore,M0_pin,M0_outlet))
+Jwn<-as.vector(rbind(Jwn_inlet,Jwn_pore,Jwn_pin,Jwn_outlet))
+
+Vol=max(M0)
+sw<-M0/Vol
+
+p1<-ggplot()+geom_line(aes(time,M0_inlet,color="inlet"))+
+	geom_line(aes(time+100,M0_pore,color="pore"))+
+	geom_line(aes(time+200,M0_pin,color="pin"))+
+	geom_line(aes(time+300,M0_outlet,color="outlet"))
+
+p2<-ggplot()+geom_line(aes(sw,Jwn))+GGTHEME+
+	xlab(expression(s^w))+ylab(expression(p^c))
+