[VB.NET 4.0] Can you speed up this function? :)

[DragonQ]

I have some code for calculating a "polytrope" (a good model for some types of star). I was reading some optimisation threads on this forum today and it made me curious enough to try to make it faster. It was already probably the most efficient program I'd made since it was my first truly multithreaded application (that actually benefited hugely from it).

Moving "Math.PI" into a constant so it's only called once didn't make any measurable difference in the time taken. However, replacing two occurrences of "r ^ 2" with a new variable initialised to "r ^ 2", which was then used twice, improved the speed of the code by 35%!

I don't actually need this code as it's from an old program but I'm just curious to see if there are any other efficiency improvements one could make?

P.S. This is running in 64-bit mode so I don't think changing any of the Doubles into Singles will improve speed.

Public Const G As Double = 0.00000000006674     'gravitational constant
    Public Const Pi As Double = Math.PI             'Pi
 
    Public Function CalculatePolytrope(ByVal n As Double, ByVal StepSize As Double, ByVal Rho As Double, ByVal P As Double) As Double()
 
        Dim K As Double = P / (Rho ^ ((n + 1) / n)) 'constant of polytrope model
        Dim m As Double = 0                         'mass
        Dim r As Double = 0 - StepSize              'radius
        Dim dP_dr As Double                         'dP/dr
        Dim dm_dr As Double                         'dm/dr
        Dim rSquared As Double                      'r^2
 
        'Process loop until Rho or P are equal to (or less than) zero:
        Do
            'Increment r
            r += StepSize
 
            'This prevents dP/dr and dm/dr becoming NaN ("not a number") when r = 0
            If r = 0 Then
                dP_dr = 0
                dm_dr = 0
            Else
                rSquared = r ^ 2
                dP_dr = -(G * m * Rho) / rSquared
                dm_dr = 4 * Pi * Rho * rSquared
            End If
 
            'Calculate new mass and pressure
            m += StepSize * dm_dr
            P += StepSize * dP_dr
 
            'This prevents square rooting a minus number, which would make Rho NaN ("not a number")
            If P > 0 Then
                Rho = (P / K) ^ (n / (n + 1))
            Else
                Rho = 0
            End If
        Loop Until Rho <= 0 OrElse P <= 0
 
        'Create an array of doubles to hold our final M and R then return them to the calling subroutine.
        Dim Results() As Double = {m, r}
        Return Results
 
    End Function