Tony YaBB Administrator Posts: 1058 Gender: Re: Runge Kutta 4th order integrator Reply #2 - 10/30/11 at 18:42:52   Sorry, Wil.  I didn't mean to ignore your post.  When I get back into programming mode, I'll try your suggestion.   The lastest beta version contains 4 integrators:  Euler-Cromer, Runge Kutta 4, Leap frog, and Verlet.  Here's the source code for Verlet and Leap frog:   Verlet: Code:```    If Not firstTimeThrough Then        For k = 1 To objects            a1x(k) = 0            a1y(k) = 0            a1z(k) = 0        Next k        For k = 1 To objects            Mu = ObjMass(k) * G            For j = k + 1 To objects                Mu2 = ObjMass(j) * G                dx = (objx(j) - objx(k))                dy = (objy(j) - objy(k))                dz = (objz(j) - objz(k))                D = Sqr(dx * dx + dy * dy + dz * dz)                If Mu > 0 Then 'ignore massless particles                    a1 = -(1 / D) * (1 / D) * Mu                    a1x(j) = a1x(j) + (dx / D) * a1                    a1y(j) = a1y(j) + (dy / D) * a1                    a1z(j) = a1z(j) + (dz / D) * a1                End If                If Mu2 > 0 Then 'ignore massless particles                    a1 = -(1 / D) * (1 / D) * Mu2                    a1x(k) = a1x(k) - (dx / D) * a1                    a1y(k) = a1y(k) - (dy / D) * a1                    a1z(k) = a1z(k) - (dz / D) * a1                End If            Next j        Next k        StepEven = False        firstTimeThrough = True    End If    If StepEven Then        For k = 1 To objects            a1x(k) = 0            a1y(k) = 0            a1z(k) = 0        Next k        For k = 1 To objects            objx(k) = objx(k) + objvx(k) * timescale + 0.5 * a2x(k) * timescale * timescale            objy(k) = objy(k) + objvy(k) * timescale + 0.5 * a2y(k) * timescale * timescale            objz(k) = objz(k) + objvz(k) * timescale + 0.5 * a2z(k) * timescale * timescale        Next k        For k = 1 To objects            Mu = ObjMass(k) * G            For j = k + 1 To objects                Mu2 = ObjMass(j) * G                dx = (objx(j) - objx(k))                dy = (objy(j) - objy(k))                dz = (objz(j) - objz(k))                D = Sqr(dx * dx + dy * dy + dz * dz)                If Mu > 0 Then 'ignore massless particles                    a1 = -(1 / D) * (1 / D) * Mu                    a1x(j) = a1x(j) + (dx / D) * a1                    a1y(j) = a1y(j) + (dy / D) * a1                    a1z(j) = a1z(j) + (dz / D) * a1                End If                If Mu2 > 0 Then 'ignore massless particles                    a1 = -(1 / D) * (1 / D) * Mu2                    a1x(k) = a1x(k) - (dx / D) * a1                    a1y(k) = a1y(k) - (dy / D) * a1                    a1z(k) = a1z(k) - (dz / D) * a1                End If            Next j        Next k        For k = 1 To objects            objvx(k) = objvx(k) + 0.5 * (a1x(k) + a2x(k)) * timescale            objvy(k) = objvy(k) + 0.5 * (a1y(k) + a2y(k)) * timescale            objvz(k) = objvz(k) + 0.5 * (a1z(k) + a2z(k)) * timescale        Next k    Else        For k = 1 To objects            a2x(k) = 0            a2y(k) = 0            a2z(k) = 0        Next k        For k = 1 To objects            objx(k) = objx(k) + objvx(k) * timescale + 0.5 * a1x(k) * timescale * timescale            objy(k) = objy(k) + objvy(k) * timescale + 0.5 * a1y(k) * timescale * timescale            objz(k) = objz(k) + objvz(k) * timescale + 0.5 * a1z(k) * timescale * timescale        Next k        For k = 1 To objects            Mu = ObjMass(k) * G            For j = k + 1 To objects                Mu2 = ObjMass(j) * G                dx = (objx(j) - objx(k))                dy = (objy(j) - objy(k))                dz = (objz(j) - objz(k))                D = Sqr(dx * dx + dy * dy + dz * dz)                If Mu > 0 Then 'ignore massless particles                    a2 = -(1 / D) * (1 / D) * Mu                    a2x(j) = a2x(j) + (dx / D) * a2                    a2y(j) = a2y(j) + (dy / D) * a2                    a2z(j) = a2z(j) + (dz / D) * a2                End If                If Mu2 > 0 Then 'ignore massless particles                    a2 = -(1 / D) * (1 / D) * Mu2                    a2x(k) = a2x(k) - (dx / D) * a2                    a2y(k) = a2y(k) - (dy / D) * a2                    a2z(k) = a2z(k) - (dz / D) * a2                End If            Next j        Next k        For k = 1 To objects            objvx(k) = objvx(k) + 0.5 * (a1x(k) + a2x(k)) * timescale            objvy(k) = objvy(k) + 0.5 * (a1y(k) + a2y(k)) * timescale            objvz(k) = objvz(k) + 0.5 * (a1z(k) + a2z(k)) * timescale        Next k            End If 'if stepeven    StepEven = Not StepEven ```   Leap frog: Code:```    If Not firstTimeThrough Then        For k = 1 To objects            a1x(k) = 0            a1y(k) = 0            a1z(k) = 0        Next k        For k = 1 To objects           Mu = ObjMass(k) * G           For j = k + 1 To objects               Mu2 = ObjMass(j) * G               dx = (objx(j) - objx(k))               dy = (objy(j) - objy(k))               dz = (objz(j) - objz(k))               D = Sqr(dx * dx + dy * dy + dz * dz)               If Mu > 0 Then 'ignore massless particles                   a1 = -(1 / D) * (1 / D) * Mu                   a1x(j) = a1x(j) + (dx / D) * a1                   a1y(j) = a1y(j) + (dy / D) * a1                   a1z(j) = a1z(j) + (dz / D) * a1               End If               If Mu2 > 0 Then 'ignore massless particles                   a1 = -(1 / D) * (1 / D) * Mu2                   a1x(k) = a1x(k) - (dx / D) * a1                   a1y(k) = a1y(k) - (dy / D) * a1                   a1z(k) = a1z(k) - (dz / D) * a1               End If           Next j       Next k       StepEven = True       firstTimeThrough = True    End If    If StepEven Then        For k = 1 To objects            a2x(k) = 0            a2y(k) = 0            a2z(k) = 0        Next k        For k = 1 To objects            objx(k) = objx(k) + objvx(k) * timescale + timescale * timescale * a1x(k) * 0.5            objy(k) = objy(k) + objvy(k) * timescale + timescale * timescale * a1y(k) * 0.5            objz(k) = objz(k) + objvz(k) * timescale + timescale * timescale * a1z(k) * 0.5        Next k        For k = 1 To objects            Mu = ObjMass(k) * G            For j = k + 1 To objects                Mu2 = ObjMass(j) * G                dx = (objx(j) - objx(k))                dy = (objy(j) - objy(k))                dz = (objz(j) - objz(k))                D = Sqr(dx * dx + dy * dy + dz * dz)                If Mu > 0 Then 'ignore massless particles                    a2 = -(1 / D) * (1 / D) * Mu                    a2x(j) = a2x(j) + (dx / D) * a2                    a2y(j) = a2y(j) + (dy / D) * a2                    a2z(j) = a2z(j) + (dz / D) * a2                End If                If Mu2 > 0 Then 'ignore massless particles                    a2 = -(1 / D) * (1 / D) * Mu2                    a2x(k) = a2x(k) - (dx / D) * a2                    a2y(k) = a2y(k) - (dy / D) * a2                    a2z(k) = a2z(k) - (dz / D) * a2                End If            Next j        Next k    Else 'stepeven        For k = 1 To objects            a1x(k) = 0            a1y(k) = 0            a1z(k) = 0        Next k        For k = 1 To objects            objx(k) = objx(k) + objvx(k) * timescale + timescale * timescale * a2x(k) * 0.5            objy(k) = objy(k) + objvy(k) * timescale + timescale * timescale * a2y(k) * 0.5            objz(k) = objz(k) + objvz(k) * timescale + timescale * timescale * a2z(k) * 0.5        Next k        For k = 1 To objects            Mu = ObjMass(k) * G            For j = k + 1 To objects                Mu2 = ObjMass(j) * G                dx = (objx(j) - objx(k))                dy = (objy(j) - objy(k))                dz = (objz(j) - objz(k))                D = Sqr(dx * dx + dy * dy + dz * dz)                If Mu > 0 Then 'ignore massless particles                    a1 = -(1 / D) * (1 / D) * Mu                    a1x(j) = a1x(j) + (dx / D) * a1                    a1y(j) = a1y(j) + (dy / D) * a1                    a1z(j) = a1z(j) + (dz / D) * a1                End If                If Mu2 > 0 Then 'ignore massless particles                    a1 = -(1 / D) * (1 / D) * Mu2                    a1x(k) = a1x(k) - (dx / D) * a1                    a1y(k) = a1y(k) - (dy / D) * a1                    a1z(k) = a1z(k) - (dz / D) * a1                End If            Next j        Next k    End If ' stepeven        For k = 1 To objects            objvx(k) = objvx(k) + timescale * (a1x(k) + a2x(k)) * 0.5            objvy(k) = objvy(k) + timescale * (a1y(k) + a2y(k)) * 0.5            objvz(k) = objvz(k) + timescale * (a1z(k) + a2z(k)) * 0.5        Next k    StepEven = Not StepEven ``` Back to top IP Logged