Gravity Balloon pressure, radius, and mass relations

Function Gval()
  Gval = 6.67384E-11
End Function
Function Pi()
  Pi = 3.14159265358979
End Function
' Functions based on (M,R,t) triple
Function M_Rt(R As Double, t As Double, rho As Double) As Double
  M_Rt = (4# / 3#) * Pi() * rho * ((R + t) ^ 3 - R ^ 3)
End Function
Function R_Mt(M As Double, t As Double, rho As Double) As Double
  R_Mt = (-t + (t ^ 2 - 4 * ((t ^ 2) / 3 - M / (4 * Pi() * rho * t))) ^ (1 / 2)) / (2)
End Function
Function t_RM(R As Double, M As Double, rho As Double) As Double
  t_RM = (M / ((4 / 3) * Pi() * rho) + R ^ 3) ^ (1 / 3) - R
End Function
Function dtdR_RM(R As Double, M As Double, rho As Double) As Double
  dtdR_RM = R ^ 2 / (3 * M / (4 * Pi() * rho) + R ^ 3) ^ (2 / 3) - 1
End Function
' Composite function
Function dPdR_RM(R As Double, M As Double, rho As Double) As Double
  dPdR_RM = dPdR_Rt(R, t_RM(R, M, rho), rho) + dPdt_Rt(R, t_RM(R, M, rho), rho) * dtdR_RM(R, M, rho)
End Function
' Functions based on (P,R,t) triple
Function P_Rt(R As Double, t As Double, rho As Double) As Double
  P_Rt = (2# / 3) * Pi() * Gval() * (t * rho) ^ 2 * (3 * R + t) / (R + t)
End Function
Function dPdt_Rt(R As Double, t As Double, rho As Double) As Double
  dPdt_Rt = (4# / 3) * rho ^ 2 * Pi() * Gval() * t * (3 * R ^ 2 + 3 * R * t + t ^ 2) / (R + t) ^ 2
End Function
Function dPdR_Rt(R As Double, t As Double, rho As Double) As Double
  dPdR_Rt = 4 * t ^ 3 * Gval() * Pi() * rho ^ 2 / (3 * (R + t) ^ 2)
End Function
Function R_Pt(P As Double, t As Double, rho As Double) As Double
  R_Pt = (1# / 3) * t * (3 * P - 2 * rho ^ 2 * Pi() * Gval() * t ^ 2) / (-P + 2 * rho ^ 2 * Pi() * Gval() * t ^ 2)
End Function
Function t_RP(R As Double, P As Double, rho As Double) As Double
  Dim n As Integer, delta As Double
  t_RP = 1.366 * (P / (2 * Gval() * Pi())) ^ (1# / 2) / rho
  n = 1
  Do
    delta = -(P_Rt(R, t_RP, rho) - P) / dPdt_Rt(R, t_RP, rho)
    t_RP = t_RP + delta
    n = n + 1
    If (Math.Abs(delta / t_RP) < 0.000001) Then
      Exit Do
    ElseIf (n > 20) Then
      t_RP = -2
      Exit Do
    End If
  Loop
End Function
' Functions based on the (P,R,M) triple
Function M_RP(R As Double, P As Double, rho As Double) As Double
  M_RP = M_Rt(R, t_RP(R, P, rho), rho)
End Function
Function P_RM(R As Double, M As Double, rho As Double) As Double
  P_RM = P_Rt(R, t_RM(R, M, rho), rho)
End Function
Function R_MP(M As Double, P As Double, rho As Double) As Double
  Dim n As Integer, delta As Double
  R_MP = R_Mt(M, (P / (2 * Pi() * Gval())) ^ (1 / 2) / rho, rho) - 0.5 * (P / (2 * Pi() * Gval())) ^ (1 / 2) / rho
  n = 1
  Do
    delta = -(P_RM(R_MP, M, rho) - P) / dPdR_RM(R_MP, M, rho)
    R_MP = R_MP + delta
    n = n + 1
    If (Math.Abs(delta * rho / (P / (2 * Gval() * Pi()))) ^ (1 / 2) < 0.000001) Then ' divide by est. value of t
      Exit Do
    ElseIf (n > 20) Then
      R_MP = -2
      Exit Do
    End If
  Loop
End Function