﻿
```

Sub Methode_2_PP(pp)

Dim j, istart1 As Integer
Dim dx, dy, r2, a, r, a2
Dim loopc As Integer

Inkey\$ = ""
istart1 = istart + pp - 1
amax = 0
loopc = 0
sync0(pp) = 0 : sync1(pp) = 0
For j = istart1 To totstar Step np
kax0(j) = 0 : kay0(j) = 0
kax1(j) = 0 : kay1(j) = 0
kax2(j) = 0 : kay2(j) = 0
kax3(j) = 0 : kay3(j) = 0
kbx0(j) = 0 : kby0(j) = 0
kbx1(j) = 0 : kby1(j) = 0
kbx2(j) = 0 : kby2(j) = 0
kbx3(j) = 0 : kby3(j) = 0
loopc = loopc + 1
Next j

nnout(pp) = nnout(pp) + loopc

Synchronisation(pp, 1)
' calculate ka0

For j = istart To totstar
For k = istart1 To totstar Step np
If j <> k Then
dx = x0(j) - x0(k)
dy = y0(j) - y0(k)
r2 = dx * dx + dy * dy
r = Sqrt(r2)
a = G * m(j) / r2
If mond = 1 Then
a2 = a * monda0
If a2 > 0 Then
a = Sqrt(a2)
Else
a = -Sqrt(-a2)
End If
End If
kax0(k) = kax0(k) + h * a * dx / r
kay0(k) = kay0(k) + h * a * dy / r
End If
Next k
Next j

' calculate kb0
For j = istart1 To totstar Step np
kbx0(j) = h * vx0(j)
kby0(j) = h * vy0(j)
Next j

Synchronisation(pp, 2)
' calculate ka1
For j = istart To totstar
For k = istart1 To totstar Step np
If j <> k Then
dx = x0(j) + kbx0(j) / 2 - x0(k) - kbx0(k) / 2
dy = y0(j) + kby0(j) / 2 - y0(k) - kby0(k) / 2
r2 = dx * dx + dy * dy
r = Sqrt(r2)
a = G * m(j) / r2
If a > amax Then amax = a
If mond = 1 Then
a2 = a * monda0
If a2 > 0 Then
a = Sqrt(a2)
Else
a = -Sqrt(-a2)
End If
End If
kax1(k) = kax1(k) + h * a * dx / r
kay1(k) = kay1(k) + h * a * dy / r
End If
Next k
Next j

' calculate kb1
For j = istart1 To totstar Step np
kbx1(j) = h * (vx0(j) + kax0(j) / 2)
kby1(j) = h * (vy0(j) + kay0(j) / 2)
Next j

Synchronisation(pp, 3)
' calculate ka2
For j = istart To totstar
For k = istart1 To totstar Step np
If j <> k Then
dx = x0(j) + kbx1(j) / 2 - x0(k) - kbx1(k) / 2
dy = y0(j) + kby1(j) / 2 - y0(k) - kby1(k) / 2
r2 = dx * dx + dy * dy
r = Sqrt(r2)
a = G * m(j) / r2
If mond = 1 Then
a2 = a * monda0
If a2 > 0 Then
a = Sqrt(a2)
Else
a = -Sqrt(-a2)
End If
End If
kax2(k) = kax2(k) + h * a * dx / r
kay2(k) = kay2(k) + h * a * dy / r
End If
Next k
Next j

' calculate kb2
For j = istart1 To totstar Step np
kbx2(j) = h * (vx0(j) + kax1(j) / 2)
kby2(j) = h * (vy0(j) + kay1(j) / 2)
Next j

Synchronisation(pp, 4)
' calculate ka3
For j = istart To totstar
For k = istart1 To totstar Step np
If j <> k Then
dx = x0(j) + kbx2(j) - x0(k) - kbx2(k)
dy = y0(j) + kby2(j) - y0(k) - kby2(k)
r2 = dx * dx + dy * dy
r = Sqrt(r2)
a = G * m(j) / r2
If a > amax Then amax = a
If mond = 1 Then
a2 = a * monda0
If a2 > 0 Then
a = Sqrt(a2)
Else
a = -Sqrt(-a2)
End If
End If
kax3(k) = kax3(k) + h * a * dx / r
kay3(k) = kay3(k) + h * a * dy / r
End If
Next k
Next j

' calculate kb3
For j = istart1 To totstar Step np
kbx3(j) = h * (vx0(j) + kax2(j))
kby3(j) = h * (vy0(j) + kay2(j))
Next j

Synchronisation(pp, 5)
'    y i+1  = y i   + (ka0     + 2 * ka1     + 2 * ka2     + ka3) / 6
For j = istart1 To totstar Step np
vx1(j) = vx0(j) + (kax0(j) + 2 * kax1(j) + 2 * kax2(j) + kax3(j)) / 6
vy1(j) = vy0(j) + (kay0(j) + 2 * kay1(j) + 2 * kay2(j) + kay3(j)) / 6
Next j

'   z i+1 = z i   + (kb0     + 2 * kb1     + 2 * kb2     + kb3) / 6
For j = istart1 To totstar Step np
x1(j) = x0(j) + (kbx0(j) + 2 * kbx1(j) + 2 * kbx2(j) + kbx3(j)) / 6
y1(j) = y0(j) + (kby0(j) + 2 * kby1(j) + 2 * kby2(j) + kby3(j)) / 6
Next j

Synchronisation(pp, 6)
' Update t = 0 values with t = 1 values

For j = istart1 To totstar Step np
x0(j) = x1(j) : y0(j) = y1(j)
vx0(j) = vx1(j) : vy0(j) = vy1(j)
Next j

End Sub

```