rho_wd Function

  DOUBLE PRECISION FUNCTION rho_wd(pr, rho_left_bracket, rho_right_bracket)

    !********************************************
    !
    !# Degenerate density as a function of pressure
    !
    !  Use bisection to find the value of the density,
    !  given the pressure
    !
    !  FT 19.12.2022
    !
    !********************************************

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN):: pr
    DOUBLE PRECISION, INTENT(IN):: rho_left_bracket
    DOUBLE PRECISION, INTENT(IN):: rho_right_bracket

    INTEGER,          PARAMETER:: tolerance_magnitude= 3
    DOUBLE PRECISION, PARAMETER:: tolerance_pr       = 1.D-5
    DOUBLE PRECISION, PARAMETER:: tolerance_rho      = 1.D-10
    DOUBLE PRECISION, PARAMETER:: pr_min             = 1.D-30

    INTEGER:: cnt
    DOUBLE PRECISION:: rho_left, rho_right, rho_mean, &
                       pr_left, pr_right, pr_mean, pr_cgs

    LOGICAL, PARAMETER:: debug= .FALSE.

    IF(pr < pr_min)THEN
      rho_wd= 0.D0
      RETURN
    ENDIF

    rho_left = FLOOR(LOG10(rho_left_bracket))
    rho_right= CEILING(LOG10(rho_right_bracket))
    ! Bisection in logarithmic scale, to find the approximate order of magnitude
    DO

      pr_left = pr_wd(1.D1**rho_left)
      pr_right= pr_wd(1.D1**rho_right)
      IF( pr_left <= pr .AND. pr_right > pr )THEN

        rho_mean= FLOOR((rho_left + rho_right)/2.D0)
        pr_mean = pr_wd(1.D1**rho_mean)

        IF(debug) PRINT *, FLOOR(ABS(LOG10(pr_right) - LOG10(pr_left)))
        IF(debug) PRINT *, LOG10(pr_right)
        IF(debug) PRINT *, LOG10(pr_left)

        IF( FLOOR(ABS(LOG10(pr_right) - LOG10(pr_left))) &
            <= tolerance_magnitude )THEN
          EXIT
        ENDIF

        IF(pr_mean < pr)THEN
          rho_left = rho_mean
        ELSE
          rho_right= rho_mean
        ENDIF

      ELSEIF( pr_left < pr .AND. pr_right < pr )THEN

        rho_right= rho_right + 2.D0

      ELSEIF( pr_left > pr .AND. pr_right > pr )THEN

        rho_left= rho_left - 2.D0

      ELSE

        PRINT *
        PRINT *, "** ERROR in SUBROUTINE rho_wd in MODULE wd_eos!"
        PRINT *, "rho_left=", rho_left
        PRINT *, "rho_right=", rho_right
        PRINT *, "pr_left=", pr_left
        PRINT *, "pr=", pr
        PRINT *, "pr_right=", pr_right
        PRINT *
        STOP

      ENDIF

    ENDDO

    ! Bisection in linear scale, to find the precise value
    pr_cgs   = pr
    rho_left = (1.D1**rho_left )
    rho_right= (1.D1**rho_right)
    IF(debug) PRINT *, "rho_left =", rho_left
    IF(debug) PRINT *, "rho_right=", rho_right
    IF(debug) PRINT *, "pr_cgs   =", pr_cgs
    cnt= 0
    DO

      pr_left = pr_wd(rho_left)
      pr_right= pr_wd(rho_right)
      IF( pr_left <= pr_cgs .AND. pr_right > pr_cgs )THEN

        rho_mean= (rho_left + rho_right)/2.D0
        pr_mean = pr_wd(rho_mean)
        IF(debug) PRINT *, "pr_left  =", pr_left
        IF(debug) PRINT *, "pr_right =", pr_right
        IF(debug) PRINT *, "rho_mean =", rho_mean
        IF(debug) PRINT *, "pr_mean  =", pr_mean
        IF(debug .AND. cnt > 100) &
          PRINT *,"ABS((pr_left - pr_right)/pr_right)=", &
          ABS((pr_left - pr_right)/pr_right)
        IF(debug .AND. cnt > 100) PRINT *, "pr_left  =", pr_left
        IF(debug .AND. cnt > 100) PRINT *, "pr_right =", pr_right
        IF(debug .AND. cnt > 100) PRINT *, "pr_mean  =", pr_mean
        IF(debug .AND. cnt > 100) PRINT *, "rho_left =", rho_left
        IF(debug .AND. cnt > 100) PRINT *, "rho_right=", rho_right
        IF(debug .AND. cnt > 100) PRINT *, "rho_mean =", rho_mean
        IF(debug .AND. cnt > 100) PRINT *, "pr_cgs   =", pr_cgs
        IF(debug) STOP
        IF( ABS((pr_left - pr_right)/pr_right) < tolerance_pr &
            .OR. &
            ABS((rho_left - rho_right)/rho_right) < tolerance_rho )THEN
          EXIT
        ENDIF

        IF(pr_mean < pr_cgs)THEN
          rho_left = rho_mean
        ELSE
          rho_right= rho_mean
        ENDIF

      ELSEIF( pr_left < pr_cgs .AND. pr_right < pr_cgs )THEN

        rho_right= 2.D0*rho_right

      ELSEIF( pr_left > pr_cgs .AND. pr_right > pr_cgs )THEN

        rho_left= rho_left/2.D0

      ELSE

        PRINT *
        PRINT *, "** ERROR in SUBROUTINE rho_wd in MODULE wd_eos!"
        PRINT *, "rho_left=", rho_left
        PRINT *, "rho_right=", rho_right
        PRINT *, "pr_left=", pr_left
        PRINT *, "pr_cgs=", pr_cgs
        PRINT *, "pr_right=", pr_right
        PRINT *
        STOP

      ENDIF

      cnt= cnt + 1
      IF(cnt > 150)THEN
        STOP
      ENDIF

    ENDDO

    IF( ABS(pr - pr_wd(rho_mean))/pr > tolerance_pr )THEN

      PRINT *, "** ERROR! The value of rho_mean found by FUNCTION rho_wd ", &
               "in MODULE wd_eos, does not give a value of pr_wd(rho_mean) ", &
               "compatible with the input pressure pr, within the required ", &
               "tolerance."
      PRINT *
      PRINT *, "rho_mean=", rho_mean
      PRINT *, "pr=", pr
      PRINT *, "pr_wd(rho_mean)=", pr_wd(rho_mean)
      PRINT *, "ABS(pr - pr_wd(rho_mean))/pr", ABS(pr - pr_wd(rho_mean))/pr
      PRINT *, "tolerance_pr=", tolerance_pr
      PRINT *
      STOP

    ENDIF

    rho_wd= rho_mean

  END FUNCTION rho_wd