c---------------------------------------------------------------
	program mass_conversion_driver
	implicit none

	real*8 Mconvert,Mold,DeltaOld,DeltaNew 
 	real*8 ConcentOld,ConcentNew,concent
	real*8 Cconvert,MconvertCNew

	write(6,*) 'Original Halo Mass?'
	read*, Mold
	write(6,*) 'Original Concentration (<0 for default)?'
	read*, ConcentOld
	if (ConcentOld.lt.0) then
	  ConcentOld=concent(Mold,1d0)
	  write(6,*) ' Replacing with c=',ConcentOld
	end if

	write(6,*) 'Original Spherical Overdensity?'
	read*, DeltaOld

	write(6,*) 'Target Spherical Overdensity?'
	read*, DeltaNew

	write(6,*) 'M_Target='
	write(6,'(1E15.7)') Mconvert(Mold,DeltaOld,DeltaNew,ConcentOld)

	write(6,*) 'c_target='
	ConcentNew= Cconvert(Mold,DeltaOld,DeltaNew,ConcentOld)
	write(6,'(1E15.7)') ConcentNew

c	Use this routine if the new rather than old concentration is known 
c	write(6,*) 'M_Target='
c	write(6,'(1E15.7)') MconvertCNew(Mold,DeltaOld,DeltaNew,ConcentNew)
c


	end
c---------------------------------------------------------------


c---------------------------------------------------------------
        real*8 function Mconvert(Mold,DeltaOld,DeltaNew,ConcentOld) 
	implicit none
	real*8 Mold,DeltaOld,DeltaNew,ConcentOld,fval,rat,convinv

	rat = DeltaOld/DeltaNew
        fval = log(1.d0+ConcentOld)-ConcentOld/(1.d0+ConcentOld)
        fval = fval/rat/ConcentOld**3

        Mconvert = Mold * ((ConcentOld*convinv(fval))**(-3)/rat)

	return
	end 
c---------------------------------------------------------------


c---------------------------------------------------------------
        real*8 function Cconvert(Mold,DeltaOld,DeltaNew,ConcentOld) 
	implicit none
	real*8 Mold,DeltaOld,DeltaNew,ConcentOld,fval,rat,convinv

	rat = DeltaOld/DeltaNew
        fval = log(1.d0+ConcentOld)-ConcentOld/(1.d0+ConcentOld)
        fval = fval/rat/ConcentOld**3

        Cconvert = (convinv(fval))**(-1)

	return
	end 
c---------------------------------------------------------------



c---------------------------------------------------------------
        real*8 function MconvertCNew(Mold,DeltaOld,DeltaNew,ConcentNew)
        implicit none
        real*8 Mold,DeltaOld,DeltaNew,ConcentNew,fval,rat,convinv

        rat = DeltaOld/DeltaNew
        fval = log(1.d0+ConcentNew)-ConcentNew/(1.d0+ConcentNew)
        fval = fval*rat/ConcentNew**3

        MconvertCNew = Mold * ((ConcentNew*convinv(fval))**3/rat)

        return
        end
c---------------------------------------------------------------


c---------------------------------------------------------------
        real*8 function convinv(f)
c       fitting function from Hu&Kravtsov 
c---------------------------------------------------------------
         implicit none
         real*8 f
         real*8 a2,a3,a4,a5,p

         a2 = 0.5116
         a3 = -1.285/3.
         a4 = -3.13E-3
         a5 = -3.52E-5
         p = a3+a4*log(f)+a5*(log(f))**2
         convinv = a2*f**(2*p)+(3./4.)**2
         convinv = 1./sqrt(convinv) + 2*f
        return
        end
c---------------------------------------------------------------



c---------------------------------------------------------------
        real*8 function concent(M,a)
c	default to Bullock et al in an LCDM sigma8=0.8 cosmology
c	assumes M is in h^-1 Ms
c---------------------------------------------------------------

        implicit none
        real*8 Mstar,a,M


	 Mstar = 2.77E12  ! M_* in h^-1 Ms
	
         concent = 9.*a*(M/Mstar)**(-0.13)
         return
        end
c---------------------------------------------------------------