c------------------------------------------------------------------------------
c     split(m2,mmin,sigma,iseed,dwmax,dw,nprog,mprog) 
c     
c     subroutine to do single binary split
c     
c     input:
c     ------
c     m2:    input mass
c     mmin:  minimum mass to resolve
c     iseed: random number seed
c     dwmax: max timstep dw (w=deltc(a))
c     sigma(m,dlsigdlm): FUNCTION to compute sigma(m) & dln(sigma)/dln(m)

c     output:
c     -------
c     dw:   actual time step in units of change in delta_c  
c     nprog: (=0,1,2) number of progenitors with m>mmin after split
c     mprog(2): array containing masses of progenitors with m>mmin
c     
c     ----------------------------
c     THEORY:
c     
c     Defining q=m/m2, qmin=mmin/m2
c     (1) generate 0 or 1 progenitors with qmin<q<1/2 according to the 
c     probability distribn dn/dq given by PS
c     (2) generate the complement q'=1-q
c     (3) compute the mean mass fraction f in progenitors with q<qmin according
c     to PS, and subtract this from the larger progenitor, so final q'=1-q-f
c     (4) discard any progenitors with q<qmin
c     
c     dn/dq  = (2/pi)^(1/2) alpha(q*m2)/q^2  dw
c     * sigma(q*m2)^2/(sigma(q*m2)^2-sigma(m2)^2)^(3/2)
c     
c     where alpha(m) = -dln(sigma)/dln(m)
c     
c     dn/dq = P(q) * R(q)  where
c     
c     P(q) = (2/pi)^(1/2)*alpha(m2/2)*B*dw /  q^(2-beta)                 
c     
c     and 
c     
c     R(q) = 1/(B*q^beta) *sigma(q*m2)^2/(sigma(q*m2)^2-sigma(m2)^2)^(3/2)
c     * alpha(q*m2)/alpha(m2/2)
c     
c     dn/dq is written this way so that P(q) is a simple power-law
c     while R(q)<1 for 0<q<1/2 . 
c     
c     If alpha(m)>0 and d(alpha)/dm>0, then the factor
c     v(q) = sigma(q*m2)^2/(sigma(q*m2)^2-sigma(m2)^2)^(3/2)
c     is monotonically increasing and ln(v) vs ln(q) is concave upwards 
c     for 0<q<1/2, so v(q) is bounded from above by a power law B*q^beta
c     chosen to equal v(q) at q=qmin and q=1/2
c     
c     alpha(q m2)/alpha(m2/2)  is also less than unity for q<1/2 .
c     
c     --------------------------------------
c     
c     Given a halo mass m2 and minimum mass mmin:
c     
c     1) Compute some values that are used repeatedly
c     
c     2) Find the power-law B*q^beta used in R(q)
c     
c     3) Compute a "timestep" dw  (dw = deltac1-deltac2 )
c     subject to two constraints specified by parameters eps1 and eps2
c     such that in this interval there is small but finite probablity
c     that the halo will have formed from the merger of two halos both 
c     of mass greater than mmin, and also the constraint dw<dwmax.
c     specifically:
c     (i) dw < eps1 *sqrt(2*(sigma(m2/2)^2-sigma(m2)^2))
c     to ensure dn/dq propto dw
c     (ii) nprog_av(qmin<q<1/2)(before rejection) < eps2 
c     (iii) dw <= dwmax to allow output at a specified w
c     
c     4) Compute f, the mean mass fraction in progenitors m<mmin
c     
c     5) Decide whether to create a fragment.
c     If yes  
c     
c     6) Generate a q from the power-law P(q)
c     
c     7) Keep or reject this q stochastically with probability R(q)
c     
c     8) If have generated a fragment q, then other fragment is q' = 1-q-f
c     If not, q' = 1-f.
c     
c     9) Discard q' if q'<qmin
c     
*****************************************************************************
c     
      subroutine split(m2,mmin,sigma,iseed,dwmax,dw,nprog,mprog) 
c     
      IMPLICIT none
      INTEGER iseed,nprog
      REAL m2,mmin,sigma,dwmax,dw,mprog(2),eps1,eps2
      EXTERNAL sigma
      COMMON /splitpars/ eps1,eps2
c     
c     calls SIGMA, RAN3 (Num.Rec.)
c     
      REAL ROOT2,RT2_PI,EPSBETA,EPSQ,dlsigdlm,sigsq_qmin,sigsq_m2,sigsq_hf,
     &     qmin,diff_qmin,diff12_qmin,diff32_qmin,v_qmin,diff_hf,diff12_hf,
     &     diff32_hf,v_hf,sfac,beta,two_pow_beta,betam1,betam1_inv,
     &     two_pow_1mbeta,qminexp,ffac,dn_dw,n_av,f,random,q_pow_betam1,
     &     q,alpha_q,sigsq_q,diff32_q,p_accept,sig_qmin,sig_m2,sig_hf,sig_q,
     &     diff12_q,b,ran3,mminlast,n_eps2,alpha_hf
      PARAMETER(ROOT2=1.4142136,RT2_PI=0.7978845608,EPSBETA=1e-6,EPSQ=1.5e-6)
c     
d     INTEGER ncall,nc1,nc2
      SAVE mminlast,sigsq_qmin
d     COMMON /monitor_split/ ncall,nc1,nc2
      DATA mminlast/0.0/
d     DATA ncall/0/
d     DATA nc1/0/
d     DATA nc2/0/
c     


c     
*******************************************************************************
cd      write(0,*) 'start split'
c     0)Since mmin will typically remain fixed we need only compute
c     sigma(mmin) and slope alpha on the first call of this routine
c     
      if (mmin.ne.mminlast) then  
         dlsigdlm=0.0           ! set to avoid compiler warning
         sig_qmin = sigma(mmin,dlsigdlm)
         sigsq_qmin = sig_qmin**2 
         mminlast=mmin
      end if
d     ncall=ncall+1       !count number of calls to this routine
c     -------------------------- 
c     1)Compute some useful numbers
c     
      sig_m2 = sigma(m2,dlsigdlm) !sigma(m2)**2    and slope alpha
      sigsq_m2 = sig_m2**2 
      sig_hf = sigma(0.5*m2,dlsigdlm) !sigma(m2/2) and slope alpha_2
      sigsq_hf = sig_hf**2
      alpha_hf = -dlsigdlm
      qmin = mmin/m2            !minimum mass ratio qmin
c     
      if(qmin.lt.(0.5-EPSQ)) then      !generate fragment with qmin<q<1/2
c     
c     2)Find B and beta such that B*q^beta = v(q) at q=qmin and q=1/2
c     where v(q) = sigma(q*m2)^2/(sigma(q*m2)^2-sigma(m2)^2)^(3/2)
c     
         diff_qmin = sigsq_qmin - sigsq_m2
         diff12_qmin = sqrt(diff_qmin)
         diff32_qmin = diff_qmin*diff12_qmin
cd        if (diff32_qmin.le.0) stop 'split() diff32_qmin.le.0 !'
         v_qmin = sigsq_qmin/diff32_qmin
c     
         diff_hf = sigsq_hf - sigsq_m2
         diff12_hf = sqrt(diff_hf)
         diff32_hf = diff_hf*diff12_hf
cd        if (diff32_hf.le.0) stop 'split() diff32_hf.le.0 !'
         v_hf = sigsq_hf/diff32_hf

c     sfac = sqrt(2*(sigma(m2/2)^2 - sigma(m2)^2)), used in (3) below
         sfac = ROOT2*diff12_hf   
c     
cd        if (v_hf.le.0) stop 'split() v_hf.le.0 !'
cd        if (qmin.le.0) stop 'split() qmin.le.0 !'
         beta = log(v_qmin/v_hf)/log(2*qmin)
         if (beta.le.0.0) then
            write(0,*) 'qmin=',qmin, ' log(2*qmin)=',log(2*qmin)
            write(0,*) 'v_qmin=',v_qmin
            write(0,*) 'v_hf=',v_hf
            write(0,*) 'log(v_qmin/v_hf)=',log(v_qmin/v_hf)
            write(0,*) 'beta=',beta
            write(0,*) 'This can probably be avoided by increasing EPSQ'
            stop ' split() in split.f: beta<0'
         end if
         two_pow_beta = 2**beta
         b = v_hf*two_pow_beta
         betam1 = beta-1.0
cd        if (two_pow_beta.le.0) stop 'split() two_pow_beta.le.0 !'
         two_pow_1mbeta = 2/two_pow_beta ! 2^(1-beta)
c     
c     3)Set time step dw and compute mean number of progeny before rejection
c     NB: special case beta=1
         if(abs(betam1).gt.EPSBETA) then ! beta != 1
            betam1_inv = 1/betam1
            qminexp = qmin**betam1 !used again in section 4)
            ffac  = two_pow_1mbeta - qminexp !used again in section 4)
            dn_dw =  RT2_PI*alpha_hf*b*betam1_inv*ffac
         else                   ! beta = 1
            dn_dw =  -RT2_PI*alpha_hf*b*log(2*qmin)
         end if
c     
         if (dn_dw.gt.0.0) then
            n_eps2=eps2/dn_dw
         else
            n_eps2=dwmax
         end if
         dw = min(eps1*sfac,n_eps2,dwmax)
         n_av = dn_dw*dw
c     
c     4)Compute the mean fraction of mass in objects of mass less than mmin
cd        if (diff12_qmin.le.0) stop 'split() diff12_qmin.le.0 !'
         f = RT2_PI*dw/diff12_qmin
c     
c     5)Randomly choose one or zero progeny depending on n_av
         random = ran3(iseed)
         if (random.le.n_av) then
d           nc1 = nc1 + 1 
c     
c     6)Select a value of q with a power-law q**(beta-2) distribution
c     NB: beta=1 is special case
            random = ran3(iseed)
c     
            if(abs(betam1).gt.EPSBETA) then ! beta != 1
               q_pow_betam1 = qminexp + random*ffac ! q^(beta-1)
               q = q_pow_betam1**betam1_inv
            else                ! beta = 1
               q = qmin*(2*qmin)**(-random)
            end if
c     
c     7)Compute rejection probability
            sig_q = sigma(q*m2,dlsigdlm) !sigma(q m2), alpha_q
            sigsq_q = sig_q**2
            alpha_q = -dlsigdlm
            diff12_q = sqrt(sigsq_q - sigsq_m2)
            diff32_q = diff12_q**3
c     acceptance probability = p_accept<1
            p_accept = sigsq_q*alpha_q/(b*q*q_pow_betam1*diff32_q*alpha_hf)
c     
            random=ran3(iseed)
            if (random.ge.p_accept) q=0.0 !reject fragment
d        if (random.lt.p_accept) nc2=nc2+1 !count number of acceptances
         else
            q=0.0               !no fragment
         end if

c     
c     calc mass of other progenitor in pair, and count how many with m>mmin
         mprog(2) = q*m2
         if(mprog(2).le.mmin) then
            nprog = 1
         else
            nprog = 2
         end if
c     
         mprog(1) = (1-q-f)*m2  !adjust for mass frac in progenitors m<mmin
         if(mprog(1).le.mmin) then
            nprog = nprog - 1
            mprog(1) = mprog(2)
            mprog(2) = 0
         end if
c     
      else                      ! qmin>1/2
         diff_hf = sigsq_hf - sigsq_m2
         diff12_hf = sqrt(diff_hf)
         sfac = ROOT2*diff12_hf   
         dw = min(eps1*sfac,dwmax) ! improve this?
         diff_qmin = sigsq_qmin - sigsq_m2
        if (diff_qmin.lt.0.0) then
c            diff_qmin<0 which should be impossible but may have
c            occured as within rounding error m2=mmin
             diff_qmin=0.0
         end if
         diff12_qmin = sqrt(diff_qmin)
         if(diff12_qmin.gt.RT2_PI*dw)then
cd        if (diff12_qmin.le.0) stop 'split() diff12_qmin.le.0 !' 
            f = RT2_PI*dw/diff12_qmin
         else
            f = 1
         end if
         mprog(1) = (1-f)*m2    !adjust for mass frac in progenitors m<mmin
         if(mprog(1).gt.mmin) then
            nprog = 1
         else
            nprog = 0
            mprog(1) = 0
         end if
         mprog(2) = 0
      end if
c     
cd      write(0,*) 'done split'
      return
      end