C#######calculate the probabilities of pij depending on psi1 and psi2    

      subroutine mcpij(th,th1,th2,psi1,psi2,p)
      implicit double precision(a-h,o-z)
      dimension p(4)

      ps1=psi1-1
      ps2=psi2-1

C----- p=(p00,p10,p01,p11)

      if (dabs(ps1).gt.1.0e-10.and.dabs(ps2).gt.1.0e-10) then

      d0=dsqrt(1+ps1*(psi1*(th-th1)**2-(th+th1)**2+2*(th+th1)))

      d=dsqrt(1+ps1*(psi1*(th1-th2)**2-(th1+th2)**2+2*(th1+th2)))

      d1=dsqrt(th1**2+(ps2/(4*ps1**2))*(ps2*(d0-d+ps1*(th2-th))**2-
     *  4*(d-1)*(d0-1)+4*ps1*((d-1)*th+(d0-1)*th2)+
     *  4*ps1**2*(th1**2-th*th2)))

      d2=dsqrt((1-th1)**2+(ps2/(4*ps1**2))*(ps2*(d-d0+ps1*(th2-th))**2-
     *   4*(1-d)*(1-d0)+4*ps1*((d-1)*(1-th)+(d0-1)*(1-th2))+
     *   4*ps1**2*((1-th1)**2+(1-th)*(th2-1)) ) )

      p(1)=(ps2*((d0-d)+ps1*(th-th2))+2*ps1*(th1-1+d2))/
     *   (2*ps2*(1-d+ps1*(2-th1-th2)))

      p(2)=(ps2*((d+d0-2)+ps1*(th-2*th1+th2))+2*ps1*(1-th1-d2))/ 
     *   (2*ps2*(d-1+ps1*(th2-th1)))

      p(3)=(ps2*(d-d0+ps1*(th-th2))+2*ps1*(d1-th1))/ 
     *   (2*ps2*(d-1+ps1*(th1-th2)))

      p(4)=(ps2*(2-d0-d+ps1*(th+2*th1+th2))-2*ps1*(d1-th1))/
     *   (2*ps2*(1-d+ps1*(th1+th2)))

      else if (dabs(ps1).gt.1.0e-10.and.dabs(ps2).lt.1.0e-10) then
C      else if (dabs(ps1).gt.1.0e-10.and.ps2.eq.0.0) then

      d0=dsqrt(1+ps1*(psi1*(th-th1)**2-(th+th1)**2+2*(th+th1)))

      p(1)=((d0-1)+ps1*(th-th1))/(2*ps1*(1-th1))

      p(2)=((d0-1)+ps1*(th-th1))/(2*ps1*(1-th1))
	
      p(3)=((1-d0)+ps1*(th+th1))/(2*ps1*th1)

      p(4)=((1-d0)+ps1*(th+th1))/(2*ps1*th1)


      else if (dabs(ps1).lt.1.0e-10.and.dabs(ps2).gt.1.0e-10) then
C      else if (ps1.eq.0.0.and.dabs(ps2).gt.1.0e-10) then

      d3=dsqrt(1+ps2**2*(th2-th)**2+ps2*(2*th+2*th2-4*th*th2))

      p(1)=((d3-1)+ps2*(th-th2))/(2*ps2*(1-th2))

      p(2)=((1-d3)+ps2*(th+th2))/(2*ps2*th2)

      p(3)=((d3-1)+ps2*(th-th2))/(2*ps2*(1-th2))

      p(4)=((1-d3)+ps2*(th+th2))/(2*ps2*th2)
       else

      p(1)=th
      p(2)=th
      p(3)=th
      p(4)=th

      end if
      return
      end