# calculate Berry and Pranger fall speed

calculate.eta_function(temp=0){
	#dynamic viscosity is only a function of temp
	# page 102 Rogers and Yau
	# temp is in K
	#eta is in kg m-1 s-1
	temp[temp<233]_temp+273.15 # convert to Kelvin
	print(paste("Kelvin temp is",temp))
	eta_1.72*10^-5*(393.0/(temp+120.0))*((temp/273.0)^1.5)
	return(eta)
}

#calculate.eta(seq(-40,30,10))

calculate.rho_function(temp=0,press=1000){
	#rho=press/(R*T)
	#rho units are kg m-3
	if (temp<233) temp_temp+273.15 # convert to Kelvin
		print(paste("Kelvin temp is",temp))
	# convert to hPa or mbar to Pascals
	press_press*100
	rho_press/(287*temp)
	return(rho)
}

#calculate.rho(0)

# from Berry and Pranger (JAM 1974, pp108-113)

calculateX_function(diameter,temp,press){
	# from eqn [3]
	#make sure diameter is in meters
	g_9.81 # m/s
	rhow_1000 # density of water
	X_(4/3)*g*rhow*calculate.rho(temp,press)*diameter^3*(calculate.eta(temp)^-2)
	return(X)
}

calculate.Re_function(diameter,temp,press){
	# calculate the Reynolds number
	#eqn [13]
	X_calculateX(diameter,temp,press)
	print(paste("X=",X))
	a20_-0.312611*10^1
	a21_0.101338*10^1
	a22_-0.191182*10^(-1)
	a23_0
	#debug from Table 1, X_1090
	lnRe_a20+a21*log(X)+a22*(log(X))^2+a23*(log(X))^3
	return(exp(lnRe))
}

calculate.Re(10^-3,0,1000)

berrypranger.fallspeed_function(diameter,temp=0,press=0){
	# eqn [1]
	# diameter must be in meters!!
	# V is in m/sec
	V_calculate.Re(diameter,temp,press)*calculate.eta(temp)/
			(diameter*calculate.rho(temp,press))
	return(V)
}

a_berrypranger.fallspeed(10^-2,20,1000) # as in Figure 3
berrypranger.fallspeed(2*10^-3,20,1013) 

# From Hess page 7
# 1 mm Hg = 1.3333mb
# 1 atmos = 1013.25 = 760 mm
testbpfallspeed_function(diameters=seq(2*10^-4,2*10^-2,10^-4),temp=20,press=1013.25){
	num_length(diameters)
	Ves_NULL
	for (i in 1:num){
		Ves_append(Ves,berrypranger.fallspeed(diameters[i],temp,press))
	}
	return(Vt=Ves,D=diameters)
}

b_testbpfallspeed()

plot(b$D/2*100,b$V*100,log="xy",xlab="Drop Radius (cm)", ylab="Vt (cm/sec)",type="l",
xlim=c(0.01,1),col=1)
abline(h=seq(100,1000,100),v=seq(0.1,5,.1),col=1)
title("Berry and Pranger Fall Speeds",col=1)
mtext(date(),side=1,outer=T,adj=1,line=-2,col=1)