include Math

# functions for Montgomery method

def cm0(b_a)
# (cM)_0
	q=exp(-PI* b_a)
	du=1
	dd=1
	i=1
	begin
		qq=2*q**(i**2)
		du+=qq
		dd+=(-1)**i*qq
		i+=1
	end while qq>0.00000001
	return log(du/dd)
end

def h(b_a)
# H(b/a)
	return PI/4/cm0(b_a)
end

def ami(b_a)
# (aM)_\infty
	r=0
	l=0
	for m in 0...100
		r+=1/(sqrt(l**2+((2*m+1)*b_a)**2))*(-1)**l
	end
	begin
		l+=1
		for m in 0...100
			r+=2/(sqrt(l**2+((2*m+1)*b_a)**2))*(-1)**l
		end
	end while l<1000
	return 2*r
end

def am(b_a,c_a)
# (aM)
	r=0.0
	n=0.0
	l=0.0
	begin
		as=sqrt((2*l+1)**2+(n/c_a)**2)
		bs=PI*as*b_a
		r+=1/(as*(exp(bs)-exp(-bs))/2)
		l+=1
	end while as<1000 or bs<100
	begin
		n+=1
		l=0
		begin
			as=sqrt((2*l+1)**2+(n/c_a)**2)
			bs=PI*as*b_a
			r+=2/(as*(exp(bs)-exp(-bs))/2)
			l+=1
		end while as<1000 or bs<100
	end while l<2
	return 2*r/c_a
end

def e(b_a,c_a)
# E(b_a,c_a)/sqrt(ab)
	cm0(b_a)/am(b_a,c_a)/sqrt(b_a)
end

def rl(r)
# function R2/R1 to b/a
	x=1.0
	dlt=0.5
	dlt*=-1 if h(x)/h(1/x)>r
	for i in 0..30
		dlt/=2
		x+=dlt while h(x)/h(1/x)<r
		x-=dlt while h(x)/h(1/x)>=r
	end
	return x
end

def rho(r1,r2,a,b,c,cond=0)
# main function
# r1,r1: resistance
# a,b,c: size in cm
# cond: 0=>thin film, 1=>rho3=rho1
	l2_l1=rl(r2/r1) # l2/l1
	r2_r1=(l2_l1*a/b)**2 # rho2/rho1
	t=0
	if cond==0
		t=c
	else
		r214=sqrt(sqrt(r2_r1))
		ab2=sqrt(a*b)
		l3_l1=1/r214*c/ab2*sqrt(l2_l1)
		t=r214*e(l2_l1,l3_l1)*ab2
	end
	r1r2=(h(l2_l1)*t*r1)**2 # rho1*rho2
	rho1=sqrt(r1r2/r2_r1)
	rho2=sqrt(r1r2*r2_r1)
print "b/a=#{l2_l1},c/a=#{l3_l1}\n"
	return [rho1,rho2]
end