subroutine binormal(x,n,nit1,nit2,c,M,MM,k1,k2,tol,pi1,pi2,mu1,mu2,v1,v2,pi1pen,pi2pen,mu1pen,mu2pen,v1pen,v2pen,mu0,v0,lrt,mlrt)
	implicit integer*4(i-k,n)
	implicit real*8(a-h,l-m,o-z)
	dimension x(1:n) 
	dimension tau1(1:n),tau2(1:n),f1(1:n),f2(1:n),f0(1:n)
	pi = 3.14159265
	i = 0
	tmp1 = 0
	tmp2 = 0
	tmp3 = 0 
  101   continue
	i = i + 1
	oldpi1 = pi1
	oldpi2 = pi2
	oldmu1 = mu1
	oldmu2 = mu2
	oldv1 = v1
	oldv2 = v2
	do 102  k = 1, n 
	       tmp1 = x(k) - mu1
		 tmp2 = x(k) - mu2
		 tmp3 = pi1*exp(-tmp1*tmp1/(2*v1))/sqrt(v1) + pi2*exp(-tmp2*tmp2/(2*v2))/sqrt(v2)
		 tau1(k) = pi1*exp(-tmp1*tmp1/(2*v1)) / (tmp3*sqrt(v1))
		 tau2(k) = 1 - tau1(k)
  102   continue
	tmp4 = 0
	tmp5 = 0
	tmp6 = 0
	do 103  j = 1, n  
		 tmp4 = tmp4 + tau1(j)
		 tmp5 = tmp5 + tau1(j)*x(j)
		 tmp6 = tmp6 + tau2(j)*x(j)
  103   continue	  
	pi1 = tmp4 / n
	pi2 = 1 - pi1
	mu1 = tmp5 / tmp4
	mu2 = tmp6 / (n - tmp4)
	if (mu1.ge.M) mu1 = M
	if (mu1.le.-M) mu1 = -M
	if (mu2.ge.M) mu2 = M
	if (mu2.le.-M) mu2 = -M
	tmp5 = 0
	tmp6 = 0
	do 104  j = 1, n
		tmp5 = tmp5 + tau1(j)*(x(j)-mu1)*(x(j)-mu1)
		tmp6 = tmp6 + tau2(j)*(x(j)-mu2)*(x(j)-mu2)
  104 continue
	v1 = tmp5/tmp4
	v2 = tmp6/(n-tmp4)
	if (v1.ge.MM) v1=MM
	if (v1.le.(1/MM)) v1=(1/MM)
	if (v2.ge.MM) v2=MM
	if (v2.le.(1/MM)) v2=(1/MM)
	dtmp1 = sqrt((mu1-oldmu1)*(mu1-oldmu1))/sqrt(mu1*mu1)
	dtmp2 = sqrt((mu2-oldmu2)*(mu2-oldmu2))/sqrt(mu2*mu2)
	dtmp3 = sqrt((pi1-oldpi1)*(pi1-oldpi1))/sqrt(pi1*pi1)
	dtmp4 = sqrt((pi2-oldpi2)*(pi2-oldpi2))/sqrt(pi2*pi2) 
	dtmp5 = sqrt((v1-oldv1)*(v1-oldv1))/sqrt(v1*v1)
	dtmp6 = sqrt((v2-oldv2)*(v2-oldv2))/sqrt(v2*v2) 
	nit1 = i
	if (i.ge.k1) dtmp4=-(dtmp1+dtmp2+dtmp3+dtmp5+dtmp6)
	if ((dtmp1+dtmp2+dtmp3+dtmp4+dtmp5+dtmp6).ge.tol) goto 101
	lrt = 0
        do 105  j = 1, n
        	tmp1 = x(j) - mu1
        	tmp2 = x(j) - mu2
        	tmp3 = x(j) - mu0
        	f1(j) = pi1*exp(-tmp1*tmp1/(2*v1))/sqrt(2*pi*v1)
        	f2(j) = pi2*exp(-tmp2*tmp2/(2*v2))/sqrt(2*pi*v2)
        	f0(j) = exp(-tmp3*tmp3/(2*v0))/sqrt(2*pi*v0)
        	lrt = lrt + log(f1(j)+f2(j))  
        	lrt = lrt - log(f0(j))
  105   continue
  	lrt = 2*lrt
	pi1pen = pi1
	pi2pen = pi2
	mu1pen = mu1
	mu2pen = mu2
	v1pen = v1
	v2pen = v2
	i = 0
	tmp1 = 0
	tmp2 = 0
	tmp3 = 0 
  106   continue
	i = i + 1
	oldpi1 = pi1pen
	oldpi2 = pi2pen
	oldmu1 = mu1pen
	oldmu2 = mu2pen
	oldv1 = v1pen
	oldv2 = v2pen
	do 107  k = 1, n 
	       tmp1 = x(k) - mu1pen
		 tmp2 = x(k) - mu2pen
		 tmp3 = pi1pen*exp(-tmp1*tmp1/(2*v1pen))/sqrt(v1pen) + pi2pen*exp(-tmp2*tmp2/(2*v2pen))/sqrt(v2pen)
		 tau1(k) = pi1pen*exp(-tmp1*tmp1/(2*v1pen)) / (tmp3*sqrt(v1pen))
		 tau2(k) = 1 - tau1(k)
  107   continue
	tmp4 = 0
	tmp5 = 0
	tmp6 = 0
	do 108  j = 1, n  
		 tmp4 = tmp4 + tau2(j)
		 tmp5 = tmp5 + tau1(j)*x(j)
		 tmp6 = tmp6 + tau2(j)*x(j)
  108   continue	
	pi1pen = 1 - ((tmp4+c) / (n+2*c))
	pi2pen = 1 - pi1pen
	mu1pen = tmp5 / (n - tmp4)
	mu2pen = tmp6 / (tmp4)
	if (mu1pen.ge.M) mu1pen = M
	if (mu1pen.le.-M) mu1pen = -M
	if (mu2pen.ge.M) mu2pen = M
	if (mu2pen.le.-M) mu2pen = -M
	tmp5 = 0
	tmp6 = 0
	do 109  j = 1, n
		tmp5 = tmp5 + tau1(j)*(x(j)-mu1pen)*(x(j)-mu1pen)
		tmp6 = tmp6 + tau2(j)*(x(j)-mu2pen)*(x(j)-mu2pen)
  109 continue
	v1pen = tmp5/(n-tmp4)
	v2pen = tmp6/tmp4
	if (v1pen.ge.MM) v1pen=MM
	if (v1pen.le.(1/MM)) v1pen=(1/MM)
	if (v2pen.ge.MM) v2pen=MM
	if (v2pen.le.(1/MM)) v2pen=(1/MM)

	dtmp1 = sqrt((mu1pen-oldmu1)*(mu1pen-oldmu1))/sqrt(mu1pen*mu1pen)
	dtmp2 = sqrt((mu2pen-oldmu2)*(mu2pen-oldmu2))/sqrt(mu2pen*mu2pen)
	dtmp3 = sqrt((pi1pen-oldpi1)*(pi1pen-oldpi1))/sqrt(pi1pen*pi1pen)
	dtmp4 = sqrt((pi2pen-oldpi2)*(pi2pen-oldpi2))/sqrt(pi2pen*pi2pen) 
	dtmp5 = sqrt((v1pen-oldv1)*(v1pen-oldv1))/sqrt(v1pen*v1pen)
	dtmp6 = sqrt((v2pen-oldv2)*(v2pen-oldv2))/sqrt(v2pen*v2pen) 
	nit2=i  
	if (i.ge.k2) dtmp4=-(dtmp1+dtmp2+dtmp3+dtmp5+dtmp6)
	if ((dtmp1+dtmp2+dtmp3+dtmp4+dtmp5+dtmp6).ge.tol) goto 106
	mlrt = 0
        do 110  j = 1, n
        	tmp1 = x(j) - mu1pen
        	tmp2 = x(j) - mu2pen
        	tmp3 = x(j) - mu0
        	f1(j) = pi1pen*exp(-tmp1*tmp1/(2*v1pen))/sqrt(2*pi*v1pen)
        	f2(j) = pi2pen*exp(-tmp2*tmp2/(2*v2pen))/sqrt(2*pi*v2pen)
        	f0(j) = exp(-tmp3*tmp3/(2*v0))/sqrt(2*pi*v0)
        	mlrt = mlrt + log(f1(j)+f2(j))  
        	mlrt = mlrt - log(f0(j))
  110 continue
        mlrt = mlrt + c*log(4*pi1pen*pi2pen)
  	mlrt = 2*mlrt
	return
	end