It is part of a 3-body solar system simulation. The original code is:

! Sum accelerations and integrate DO AXIS=1,3 ! Total accelerations of each body A(AXIS,1)=VF(AXIS,1,2)+VF(AXIS,1,3) A(AXIS,2)=VF(AXIS,2,1)+VF(AXIS,2,3) A(AXIS,3)=VF(AXIS,3,1)+VF(AXIS,3,2) ! ! Integrate for each body DO BODY1=1,3 ! Keep previous velocities for Tustin's method PV(AXIS,BODY1)=V(AXIS,BODY1) ! Trapezoidal integration for velocities V(AXIS,BODY1)=V(AXIS,BODY1)+(3*A(AXIS,BODY1)-PA(AXIS,BODY1)) 1 *DT*0.5D0 ! Tustin's method for positions S(AXIS,BODY1)=S(AXIS,BODY1)+0.5D0*(V(AXIS,BODY1)+PV(AXIS, 1 BODY1))*DT ! Keep previous accelerations for next frame PA(AXIS,BODY1)=A(AXIS,BODY1) ENDDO ENDDO