!program for predicting splitting of air, water system in a t using 
!shoham model
program shoham
	!implicit none
	real::a,r0,al,hd,ag,abg,si,d,rhoair,rhowat,b,sb,ab1,fc,k,t,viscl,h,tna
	real::deltar,vgs,vls,fbg,vl,vg,visg
	real::an,ad,at,adn,sad,tatley
	real::sa,fbl,abn,abl1,solve
	COMMON d,sb,abg,b,ab1,ad,rhoair,rhowat,viscl,visg,vgs,vls 
	print*,"enter the lighter fluid flow rate "
	read*,vgs
	print*,"enter the flow rate of heavier fluid"  
	read*,vls
	print*,"enter tube dia"
	read*,d
	print*,"enter T dia"
	read*,dt


	OPEN(UNIT=5,FILE="fbgl.xls")

	rhoair=1.18

	rhowat=998
	viscl=84e-5
	visg=1.9e-5
	hd=tatley(vls,vgs,rhoair,rhowat,viscl,visg,d)
	h=hd*d
	at=(d**2)*3.14/4.
	print*,"assume value for fbg(fraction of gas flow rate through branch)"
	!read*,fbg
	
do	fbg=0.01, 0.9,0.01
	al=(d**2)*(3.14 - ACOS(2*hd-1)+(2*hd-1)*(1-(2*hd-1)**2)**0.5)/4
	ag=at-al
	abg=fbg*ag
	si=d*sqrt(1-(2*hd-1)**2)
	b=0.5*(d-si)
	sb=d*sqrt(1-(2*b/d-1)**2)
	ab1=0.25*(d**2)*(3.14-acos(2*b/d-1)+(2*b/d-1)*sqrt(1-(2*b/d-1)**2))
	ad=solve()	
	a=ad*d
	r0=((dt**2)+(a**2))/(2*a)
	vl = vls*at/al
	vg = vgs*at/ag
	!fc=((rhowat*vl**2)-(rhoair*vg**2))/(rhowat*r0)
	!k=3*viscl/(rhowat*(h**2))
	!t=(r0*asin(dt/r0))/vl
	!deltar=fc*(t-(1-exp(-k*t))/k)/k
	!an=a-deltar
	!if(an>d)then
	!an=d
	!end if
	!adn=an/d
	tna=2.*ad-1.
	sad=sqrt(1-tna**2)
	sa=sad*d
	abn=.25*(d**2)*(3.14-acos(tna)+tna*sad)
	abl1=.5*(abn-(an-b)*sb-ab1)
	if(b.ge.an)abn=0.
	if(an.ge.(b+si))abn=al
	fbg1=abg/ag
	fbl=abl1/al
	print*,"fbl=",fbl
	print*,"fbg=",fbg1
	write(5,*)fbl,fbg1
	
	
	end do
end program shoham 


real function solve
	implicit none
	common ad
	REAL::a1=0.0,a2=1.0,AD,EPS=.00001 
	REAL::x1,x2,f0,f1,f2,x0
	REAL,EXTERNAL::F
	x1=a1
	x2=a2
	f1=f(x1)
	f2=f(x2)
	DO
		x0=(x1+x2)/2.
		f0=f(x0)
		IF(f0.eq.0.0)THEN
			solve=x0
			EXIT
		END IF
		IF(f1*f0.LT.0.0)THEN
			x2=x0
			f2=f0
		ELSE
			x1=x0
			f1=f0
		END IF
		IF(ABS((x2-x1))<EPS)THEN
			solve=(x1+x2)/2.
			EXIT
			ELSE
			CONTINUE
		END IF
	END DO
end function solve
 
 
real FUNCTION  f(ad)
	REAL::ad,d,sb,abg,b,ab1,ab
	COMMON d,sb,abg,b,ab1
	ab=0.25*d**2*(3.14-ACOS(2*ad-1)+(2*ad-1)*SQRT(1-(2*ad-1)**2))
	F =(ad*d*sb+ab-(2*abg+b*sb-ab1))
	RETURN
END	FUNCTION F


real function tatley(vls,vgs,rhoair,rhowat,viscl,visg,d )
	 implicit none
	 real::hd1,hd,TOL,rhowat,rhoair,viscl,PI=3.14,visg,towi,towl,towg,al,at,ee
	 real::fi,fg,fl,eng,cg,enl,cl,vg,vl,reg,rel,dl,dg,sl,sg,si,a,b,h,d,ag,vls,vgs
	 a=0.
	 b=d
	 at=.25*PI*d**2
	 DO 
		h = (a+b)/2  
		hd = h/d
	      ! AREA OF CROSS SECTION OF THE UPPER LAYER	 IN M2
		ag = (d**2)*(ACOS(2*hd-1)-(2*hd-1)*(1-(2*hd-1)**2)**0.5)/4
		  ! AREA OF CROSS SECTION OF THE LOWER LAYER	 IN M2
		al = (d**2)*(3.14 - ACOS(2*hd-1)+(2*hd-1)*(1-(2*hd-1)**2)**0.5)/4
		al = at - ag
			  ! PERIMETER OF UPPER LAYER IN M
		sg = d*ACOS(2*hd-1)
				! PERIMETER OF LOWER LAYER IN M
		sl = d*(3.14-ACOS(2*hd-1))
		Sl = 3.14*d - sg
			 ! PERIMETER OF INTERFACE LAYER IN M
		si = d*((1-(2*hd-1)**2)**0.5)
			! PHASE VELOCITIES, EQUALENT DIAMETER AND REYNOLDS NUMBERS
		vl = vls*at/al
		vg = vgs*at/ag
		dl = 4*al/sl
		dg = 4*ag/(sg+si)
		rel = dl*vl*rhowat/viscl
		reg = dg*vg*rhoair/visg
			! FRICTION FACTORS AND SHEAR STRESSES
		IF(rel.LT. 2100) THEN
			 cl= 16
			 enl = 1
		ELSE
			 cl= 0.046
			 enl = 0.2
		ENDIF

		IF(reg .LT. 2100) THEN
			 cg = 16
			 eng = 1
		ELSE
			 cg = 0.046
			 eng = 0.2
		ENDIF
		fl = cl/(rel**enl)
		fg = cg/(reg**eng)
		IF(fg > 0.014) THEN
		  fi = fg
		ELSE
		  fi = 0.014
		ENDIF
		ee = al/at
		towg = (fg*rhoair*(vg**2))/2
		towl = (fl*rhowat*(vl**2))/2
		towi = (fi*rhoair*(vg-vl)*ABS(vg-vl))/2
		TOL = towg*sg/ag - towl*sl/al + towi*si*((1/al)+(1/ag))			  
		IF( ABS(TOL) < 0.01)  THEN
		GOTO 10
		 ELSEIF(TOL < 0 ) THEN
		 a= h								  
		 ELSE
		 b= h
		 ENDIF
	ENDDO
		10 tatley=hd
		hd1=hd
END function tatley
