real function contrainte(x)
      vector eta
      vector rho
      vector eta_up
      vector babar
      vector radplus
      vector radminus
      vector cenplus
      vector cenminus
      common/fit_param/cenvbc,sigvbc,cenvubcb,sigvubcb,cendelmd
     &,sigdelmd,cenepsK,sigepsK,cenmt,sigmt,cenbk,sigbk,
     &cenfbb,sigfbb,cenxi_s,sigxi_s,
     &cens2a,sigs2a,cens2b,sigs2b

      common/cte_val/pi,h,gf,xmw,xmbd,xmbs,xmk,delmk,xmc,xlamb,
     &xnu1,xnu2,xnu3,xnub,fk

      logical babar_flag,dms_flag
      common /fit_drive/babar_flag,dms_flag

      etamax = eta_up(1)

*- Initialisation :
      if (x.eq.0.) then 
        pi       = acos(-1.) 
        h        = 6.582122e-13
        gf       = 1.16639e-5   !  Fermi Constant
        xmw      = 80.33        !  W Mass ( error neglected )
        xmbd     = 5.279        !  B0d Mass ( error neglected )
        xmbs     = 5.369        !  B0s Mass ( error neglected ) 
        xmk      = .4977        !  K0 Mass ( error neglected )
        delmk    = 3.510e-15    !  DmK (error neglected) 
        xmc      = 1.3          !  c quark Mass ( error neglected )
        xlamb    = .2205        !  Lambda ( error = .0018 neglected ) 
        xnu1     = 1.38         !  eta parameters for EpsK (Buras)
        xnu2     = .57
        xnu3     = .47
        xnub     = .55
        fk       = .161         !  K decay constant
*--------------------------------------------------------------------------

        open (UNIT=10,FILE='parameters.input',STATUS='OLD')
c
        read(10,1010)flg_name,dms_flag
        read(10,1010)flg_name,babar_flag
 1010   format(a14,1x,l1)
        read(10,*)
c
c Input values 
c
        read(10,1020)var_name,cenvbc,sigvbc
        read(10,1020)var_name,cenvubcb,sigvubcb
        read(10,1020)var_name,cendelmd,sigdelmd
        read(10,1020)var_name,cenepsK,sigepsK
        cenepsK = cenepsK*1.e-03
        sigepsK = sigepsK*1.e-03
        read(10,1020)var_name,cenmt,sigmt
        read(10,1020)var_name,cenbk,sigbk
        read(10,1020)var_name,cenfbb,sigfbb    
        read(10,1020)var_name,cenxi_s,sigxi_s    
        read(10,1020)var_name,cens2a,sigs2a
        read(10,1020)var_name,cens2b,sigs2b
 1020   format(1x,a9,2(2x,f8.4))
        close(10)

        if (babar_flag) babar(1) = 1. 
      endif
c
c Dmd constraint 
c
      if(x.eq.1.) then 
        xt = (cenmt/xmw)**2
        C1 = gf**2/(6*pi**2)*xlamb**2*xmw**2*(xmbd/h)*xnub*si(xt)

        dmd = cendelmd +sigdelmd
        vcb = cenvbc - sigvbc
        fbb = cenfbb - sigfbb
        C2max =  dmd/(C1*vcb**2*fbb**2)

        dmd = cendelmd - sigdelmd
        vcb = cenvbc + sigvbc
        fbb = cenfbb + sigfbb
        C2min =  dmd/(C1*vcb**2*fbb**2)
        
        xmin = 1.-sqrt(C2min)
        xmax = 1.
        step = abs(xmax-xmin)/50.
        do i = 1,50
          x = xmin + float(i-1)*step 
          rho(i) = x
          eta(i) = dmd_fc(x,c2min,etamax) 
        enddo
        rho(51) = 1.
        eta(51) =  dmd_fc(1.,c2min,etamax) 
        ilast = 51
        xmin = 1
        xmax = 1. - sqrt(c2max)
        step = abs(xmax-xmin)/48.
        do i = 1,48
          x = xmin - float(i-1)*step 
          rho(i+ilast) = x
          eta(i+ilast) = dmd_fc(x,c2max,etamax) 
        enddo
        x =  1. - sqrt(c2max)
        rho(100) =  x
        eta(100) = dmd_fc(x,c2max,etamax) 
      endif 
c
c Epsilon K  constraint 
c
      if(x.eq.2.) then 
        epsK = cenepsK
        xt    = (cenmt/xmw)**2 
        xc    = (xmc/xmw)**2 

        C1 = gf**2*fk**2*xmk*xmw**2/(6.*sqrt(2.)*pi**2*delmk)
        B1_f0 = -xnu1*si(xc) + xnu3*sij(xc,xt)
        B2_f0 = vcb**2*xnu2*si(xt)

        vcb = cenvbc + sigvbc
        bk = cenbk + sigbk 
        C2max =  C1*BK*vcb**2*xlamb**2
        B2_f0max = vcb**2*xnu2*si(xt)

        vcb = cenvbc - sigvbc
        bk = cenbk - sigbk 
        C2min =  C1*BK*vcb**2*xlamb**2
        B2_f0min = vcb**2*xnu2*si(xt)

        xmin = -1.
        xmax = 1. 
        step = abs(xmax-xmin)/50.
        do i = 1,50
          x = xmin + float(i-1)*step 
          y = fc_eps(epsk,C2max,B1_f0,b2_f0max,x,etamax)
          rho(i) = x 
          eta(i) = y 
        enddo  
        x = 1.
        y =  fc_eps(epsk,C2max,B1_f0,b2_f0max,x,etamax)
        rho(51) = x 
        eta(51) = y 
        y =  fc_eps(epsk,C2min,B1_f0,b2_f0min,x,etamax)
        rho(52) = x 
        eta(52) = y 
        ilast = 52
        xmin = 1
        xmax = -1.
        step = abs(xmax-xmin)/48.
        do i = 1,48
          x = xmin - float(i)*step 
          y = fc_eps(epsk,C2min,B1_f0,b2_f0min,x,etamax)
          rho(ilast+i) = x 
          eta(ilast+i) = y 
        enddo  
      endif  

c
c Vub/Vcb  constraint 
c
      if(x.eq.3.) then 
        valm = cenvubcb - sigvubcb
        valp = cenvubcb + sigvubcb
        xmin = -valp/xlamb
        xmax = -valm/xlamb
        step = abs(valp/xlamb-valm/xlamb)/10.
        do i = 1,10
          x = xmin + float(i)*step 
          rho(i) = x
          eta(i) = 0. 
          ilast = i 
        enddo
        xmin = -valm/xlamb
        xmax = valm/xlamb
        step=abs(xmax-xmin)/30.
        do i = 1,30 
          x = xmin + float(i)*step 
          rho(ilast+i)=x
          y = (valm/xlamb)**2-x**2
          if (y.le.0) then
            y = 0.
          else
            y = sqrt(y)
          endif
          eta(ilast+i)=y 
        enddo
        ilast = 40
        xmin = valm/xlamb
        xmax = valp/xlamb
        step=abs(xmax-xmin)/10.
        do i=1,10
          x = xmin + float(i)*step 
          rho(ilast+i) = x 
          eta(ilast+ i) = 0.
        enddo
        ilast = 50 
        xmin = -valp/xlamb
        xmax = valp/xlamb
        step = abs(xmax-xmin)/50.
        do i = 1,50
          x = xmax- float(i)*step 
          rho(ilast+i)=x
          y = (valp/xlamb)**2-x**2
          if (y.le.0) then
            y = 0.
          else
            y = sqrt(y)
          endif
          eta(ilast+i)=y 
        enddo  
      endif

c
c Sin(2B) constraint 
c
      if(x.eq.4.) then 
        b = cens2b - sigs2b
        if (b.lt.-0.99) b = -.99  
        if (b.gt.+0.99) b = +.99  
        xmin = -1.
        xmax = 1.
        step = abs(xmax-xmin)/50.
        do i = 1,50
          x = xmin + float(i-1)*step
          y = (1.-x)/rplus(b)
          rho(i) = x 
          eta(i) = y 
        enddo   
        ilast = 50
        b = cens2b + sigs2b
        if (b.lt.-0.99) b = -.99  
        if (b.gt.+0.99) b = +.99  
        xmin = -1.
        xmax = 1.
        step = abs(xmax-xmin)/50.
        do i = 1,50
          x = xmax - float(i-1)*step
          y = (1.-x)/rplus(b)
          rho(ilast+i) = x 
          eta(ilast+i) = y 
        enddo   
      endif  

c
c Sin(2a) constraint 
c
      if(x.eq.5.) then 
c Minus 1 sigma 
        a = cens2a - sigs2a
        if (a.lt.-0.99) a = -.99  
        if (a.gt.+0.99) a = +.99  
        cenminus(1) = rminus(a)/2.
        radminus(1) = .5*sqrt(1.+rminus(a)**2)
        xmin = .5 - sqrt((1.+rminus(a)**2)/4.)
        xmax = .5 + sqrt((1.+rminus(a)**2)/4.)
        step = abs(xmax-xmin)/50.
        do i = 1,50
          x = xmin + float(i-1)*step
          y1 = (1+rminus(a)**2)/4. - (.5-x)**2
          if (y1.lt.0.) y1 = 0.
          y = rminus(a)/2. + sqrt(y1)
          if (y.lt.0.) y=0.
          rho(i) = x 
          eta(i) = y 
        enddo   
        ilast = 50
        do i = 1,50
          x = xmax - float(i-1)*step
          y1 = (1+rminus(a)**2)/4. - (.5-x)**2
          if (y1.lt.0.) y1 = 0.
          y = rminus(a)/2. - sqrt(y1)
          if (y.lt.0.) y=0.
          rho(i+ilast) = x 
          eta(i+ilast) = y
        enddo   
        rho(100) = rho(1)
        eta(100) = eta(1)
c Now Plus 1 sigma 
        ilast = 100
        a = cens2a + sigs2a
        if (a.lt.-0.99) a = -.99  
        if (a.gt.+0.99) a = +.99  
        cenplus(1) = rminus(a)/2.
        radplus(1) = .5*sqrt(1.+rminus(a)**2)
        xmin = .5 - sqrt((1.+rminus(a)**2)/4.)
        xmax = .5 + sqrt((1.+rminus(a)**2)/4.)
        step = abs(xmax-xmin)/50.
        do i = 1,50
          x = xmin + float(i-1)*step
          y1 = (1+rminus(a)**2)/4. - (.5-x)**2
          if (y1.lt.0.) y1= 0. 
          y = rminus(a)/2. + sqrt(y1)
          if (y.lt.0.) y=0.
          rho(i+ilast) = x 
          eta(i+ilast) = y
        enddo   
        ilast = 150
        do i = 1,50
          x = xmax - float(i-1)*step
          y1 = (1+rminus(a)**2)/4. - (.5-x)**2
          if (y1.lt.0.) y1 = 0.
          y = rminus(a)/2. - sqrt(y1)
          if (y.lt.0.) y=0.
          rho(i+ilast) = x 
          eta(i+ilast) = y
        enddo   
        rho(200) = rho(101)
        eta(200) = eta(101)
      endif  


 999  continue
      return 
      end
*************************************************************
      real function rplus(z) 
      rplus = (1./z)*(1.+sqrt(1-z**2))
      return
      end
*************************************************************
      real function rminus(z) 
      rminus = (1./z)*(1.-sqrt(1-z**2))
      return
      end
*************************************************************
      function fc_eps(epsk,C2,B1_f0,b2_f0,x,etamax)
      y = epsk/(C2*(B1_f0+B2_f0*(1.-x)))
      if (y.gt.etamax) y = etamax
      if( y .lt.0.) y = 0.
      fc_eps = y 
      return
      end
*************************************************************
      function dmd_fc(x,cte,etamax)
      y = cte - (1.-x)**2
      if (y.le.0) then
        y = 0.
      else
        y = sqrt(y)
      endif
      if(y.gt.etamax) y = etamax
      dmd_fc = y
      return
      end
*************************************************************
      FUNCTION si(x)

      si = 1./4.+9./4.*1./(1.-x)-3./2.*1./(1.-x)**2
      si = si-3./2.*x**2/(1.-x)**3*alog(x)
      si = x*si

      return
      end


*************************************************************
      FUNCTION sij(xi,xj)

      sij1 = (1./4.+3./2.*1./(1.-xi)-3./4.*1./(1.-xi)**2)*
     &     alog(xi)/(xi-xj)
      sij2 = (1./4.+3./2.*1./(1.-xj)-3./4.*1./(1.-xj)**2)*
     &     alog(xj)/(xj-xi)
      sij3 = -3./4.*1./((1.-xi)*(1.-xj))

      sij = xi*xj*(sij1+sij2+sij3)

      return
      end