In a recent article, I have explained the Euler’s method for solving ordinary differential equationsusing as a motivation the fictionalized version in the film Hidden Figures of the scientific contribution of Katherine Goble and her two colleagues to the NASA space program. I have also described as an example of application a program in awk programming language for calculating the orbits of planets of the solar system. My interest in astrodynamics come back to my juvenile age, when still going to high school, my parents decided to gift me a more sophisticated microcomputer than my previous one (the celebrated Commodore VIC 20). So I became a programmer of a Philips MSX VG 8010 that I still jelously own in its original box. So, powered by the versatile Federico Faggin’s Zilog Z80 processor with a clock 3.58 MHz, with an impressive (for a previous owner of a VIC20 with a mere 3.583 kB!) memory of 32 kB RAM , 16kB of video RAM and a dedicated tape-record device as storage system, I started to write more sophisticated in MSX Basic. At that time, I was eagerly following the department “Ricreazioni al Computer” by the famous computer scientist A. K. Dewdney on the magazine “Le Scienze”, the Italian edition of Scientific American. The new microcomputer allowed me to experiment with the fascinating computational topics that Dewdney was offering every month. One of these topics was dedicated to the simulation of stars using the algorithm based on the Euler integration of the Newton equation. Following the instruction of Dewney, I managed to write a small program in MSX basic and this was the starting of my interest in computational astronomy.

In a nostalgic retro-programming archaeology survey of my ancient software creations, some years ago I managed to successfully recover also this program from 35 years old magnetic tapes. To retrieve these programs, I have used the opensource MSX emulator OpenMSX (https://openmsx.org) but editing the application using this emulator (not to mention using the old hardware) is not in my retro programming intentions! Fortunately, the MSX Basic is a Microsoft product, and as such, it is indeed very similar to the GW basic. Therefore GW-Basic interpreter or compilers can also run my ancient program codes. Looking on the Web, I could recently find a great freeware program called QB64 by Rob Galleon and QBTeam that run on my MacBook (and also on Linux and MS Windows) that does a good job. Being a QuickBasic compiler, it first transforms the program code from Basic to C++ language, and then compile it using the gnu C++ compiler. Therefore, it can be used for developing programs in BASIC on modern computers.

With some very little modification, I could compile with QB64 the old MSX code, but in doing that I got hooked in the pleasure of programming again in Basic language, and so I decided to update a bit the computer code to use the higher graphics screen resolution of my MacBook and to improve the algorithm and the graphic output. The final BASIC program (see Appendix) uses an algorithm very similar to the one in the P*lanetario.awk program* of my previous article but with some little variation.

The starting positions of the planets, stored in the DATA section of the program, correspond to those of the previous article (taken from the NASA ephemerides of the 1/7/1969). The program starts with a menu that allows you to change the time step (in seconds) and the length of the simulation (in days).

Once the simulation is running, it produces the following output with the trajectory of the planets projected on the planes x-y and x-z. The program allows a limited zoom of the trajectory using the keys indicated on the screen.

**HAPPY SIMULATIONS!**

### APPENDIX

The listing of the Program PLANETARIO.

' PLANETARIO ' Questo programma esegue una simulazione dei pianeti ' del sistema solare ' (c) 1986-2019 D. Roccatano '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ' DIM SHARED X(100), Y(100), Z(100), X1(100), Y1(100), Z1(100) AS DOUBLE DIM SHARED VX(100), VY(100), VZ(100), VX1(100), VY1(100), VZ1(100) AS DOUBLE DIM SHARED AX(100), AY(100), AZ(100) AS DOUBLE DIM SHARED MASS(100) AS DOUBLE DIM SHARED NAM$(100) ' PARAMETRI E COSTANTI AU = 1.496E11 IAU = 1.0 / AU LM = 10 * AU LZ = 0.25 * AU FX = 150. / (LM) FY = 150. / (LM) FZ = 150. / (LZ) G = 6.67430E-11 MU = 1.0E24 K2M = 1000 D2S = 24 * 3600 S2D = 1.0 / D2S D2Y = 1.0 / 365 VCV = AU / D2S DT = 600 SIML = 3650 ' Carica le cordinate e le velocita' delle stelle. READDATA ' Inizializza lo schermo e disegna il sistema di coordinate. R = 0 SCREEN _NEWIMAGE(800, 600, 256) MENU NSTEPS = SIML * D2S / DT DISEGNA FOR K = 1 TO N MASS(K) = MASS(K) * MU X(K) = X1(K) * AU Y(K) = Y1(K) * AU Z(K) = Z1(K) * AU VX(K) = VX1(K) * VCV VY(K) = VY1(K) * VCV VZ(K) = VZ1(K) * VCV NEXT K ' Inizia la simulazione FOR T = 1 TO NSTEPS FOR I = 1 TO N X(I) = X(I) + VX(I) * DT Y(I) = Y(I) + VY(I) * DT Z(I) = Z(I) + VZ(I) * DT AX(I) = 0 AY(I) = 0 AZ(I) = 0 FOR J = 1 TO N IF J <> I THEN D2 = (X(I) - X(J)) ^ 2 + (Y(I) - Y(J)) ^ 2 + (Z(I) - Z(J)) ^ 2 ID = 1. / SQR(D2) F = G * MASS(J) / D2 AX(I) = AX(I) + F * (X(J) - X(I)) * ID AY(I) = AY(I) + F * (Y(J) - Y(I)) * ID AZ(I) = AZ(I) + F * (Z(J) - Z(I)) * ID END IF NEXT J VX(I) = VX(I) + AX(I) * DT VY(I) = VY(I) + AY(I) * DT VZ(I) = VZ(I) + AZ(I) * DT IF (X(I) < LM AND X(I) > -LM) AND (Y(I) > -LM AND Y(I) < LM) AND (Z(I) < LZ AND Z(I) > -LZ) THEN PSET (225 + FX * X(I), 200 + FY * Y(I)), I PSET (575 + FX * X(I), 200 + FZ * Z(I)), I END IF NEXT I R = R + DT DAYS = R * S2D LOCATE 32, 5: PRINT USING "TIME LAPSE: ######.## Days -> ####.## Years"; DAYS, DAYS * D2Y K$ = INKEY$ IF K$ = "E" OR K$ = "e" THEN SIGN = 1: RESCALE IF K$ = "R" OR K$ = "r" THEN SIGN = -1: RESCALE NEXT T IF INKEY$ <> "" THEN RUN END SUB MENU SHARED DT, SIML '***** SUBROUTINE I ***** ' WHILE (HH <> 3) CLS LOCATE 5, 3: PRINT "SETUP SIMULATION PARAMETERS" LOCATE 8, 10: PRINT USING "1] TIME STEP. Current value: ####.#### (Seconds)"; DT LOCATE 9, 10: PRINT USING "2] SIMULATION LENGHT. Current value: ###### (Days)"; SIML LOCATE 10, 10: PRINT "3] RUN IT" LOCATE 14, 12: INPUT "Scegli la funzione {1} {2} {3}"; HH ' IF (HH = 1) THEN LOCATE 16, 12: INPUT "Inserisci il valore del time step in secondi: "; DT IF (HH = 2) THEN LOCATE 18, 12: INPUT "Inserisci la lunghezza della simulazione in giorni: "; SIML WEND END SUB SUB READDATA SHARED N, MASS, X1, Y1, Z1, VX1, VY1, VZ1, NAM$ ' ******* SUBROUTINE I ******* ' Coord. e vel. dei pianeti del sistema solare il 01/07/1969 alle 00:00 ' Numero di corpi : 10 N = 10 FOR I = 1 TO N READ MASS(I), X1(I), Y1(I), Z1(I), VX1(I), VY1(I), VZ1(I), NAM$(I) NEXT DATA 1.989E6,4.501422142050814E-03,8.265096431911833E-04,-6.106498754821928E-05,-3.725998016167921E-07,5.629018018191787E-06,-1.420343330057636E-08,"SOL" DATA 0.330,3.626352905971110E-01,-3.112613637670733E-02,-3.556949387212367E-02,-2.936459346589251E-03,2.929519404193396E-02,2.661199336367969E-03,"MERCURY" DATA 4.872,6.434726285492286E-01,-3.454394151201869E-01,-4.166240919470099E-02,9.516484358480345E-03,1.770408177673933E-02,-3.091459518364233E-04,"VENUS" DATA 6.04566673,1.707828894316037E-01,-1.002169564426745E+00,-1.278736821922684E-04,1.669197477793333E-02,2.756121398804654E-03,2.571053551433845E-07,"EARTH+MOON" DATA 0.642,-6.669858220909308E-02,-1.458105888577555E+00,-2.885595380269896E-02,1.450798593015611E-02,5.244375487322264E-04,-3.467167007807005E-04,"MARS" DATA 1898,-5.376385639107029E+00,-8.726835067297002E-01,1.240832480992071E-01,1.122255348949764E-03,-7.096896640081010E-03,4.095137444408094E-06,"JUPITER" DATA 568,7.884725882151399E+00,4.851428514519507E+00,-3.981641564307521E-01,-3.226967396201162E-03,4.740175258901204E-03,4.511263076949473E-05,"SATURN" DATA 86.8,-1.826473469434352E+01,-1.177885710585230E+00,2.326769071350472E-01,2.237727077038867E-04,-4.108566058385769E-03,-1.826199809165856E-05,"URANUS" DATA 102.0,-1.604710777852071E+01,-2.571063138942750E+01,8.990992317848128E-01,2.642815824514614E-03,-1.643946715455203E-03,-2.678926484808933E-05,"NEPTUNE" DATA 0.01464,-3.048234462543693E+01,2.734464418107592E+00,8.524018333082706E+00,3.230523482089385E-04,-3.313034717024442E-03,2.612068684771845E-04,"PLUTO+MOONS" END SUB SUB DISEGNA SHARED IAU, DT, SIML, NSTEP, LM, LZ 7000 ' ******* SUBROUTINE II ******* ' DISEGNA IL SISTEMA DI RIFERIMENTO CLS LOCATE 2, 45: PRINT "PLANETARIO" 'BOX intorno al contatore LINE (25, 490)-(500, 515), 15, B ' DISEGNA I DUE GRAFICI LINE (75, 50)-(375, 350), 15, B LINE (425, 50)-(725, 350), 15, B FOR I = 25 TO 325 STEP 10 LINE (50 + I, 350)-(50 + I, 350 - 3), 15 LINE (50 + I, 50)-(50 + I, 50 + 3) LINE (400 + I, 350)-(400 + I, 350 - 3), 15 LINE (400 + I, 50)-(400 + I, 50 + 3) NEXT FOR I = 300 TO 0 STEP -10 LINE (75, 50 + I)-(75 + 3, 50 + I), 15 LINE (375, 50 + I)-(375 - 3, 50 + I), 15 LINE (425, 50 + I)-(425 + 3, 50 + I), 15 LINE (725, 50 + I)-(725 - 3, 50 + I), 15 NEXT LINE (75, 200)-(375, 200), 12, , 63 LINE (225, 50)-(225, 350), 12, , 63 LINE (425, 200)-(725, 200), 12, , 63 LINE (575, 50)-(575, 350), 12, , 63 LOCATE 24, 29: PRINT "X": LOCATE 5, 8: PRINT "Y" LOCATE 24, 74: PRINT "X": LOCATE 5, 52: PRINT "Z" ' DISEGNA LEGENDA DEI COLORI DEI PIANETI LINE (500, 390)-(740, 515), 15, B LOCATE 26, 67: PRINT "COLORS OF THE PLANETS" FOR K = 1 TO 5 COLOR K: LOCATE 27 + K, 70: PRINT NAM$(K) NEXT K FOR K = 6 TO 11 COLOR K: LOCATE 22 + K, 80: PRINT NAM$(K) NEXT K LOCATE 26, 5: COLOR 15: PRINT "Coordinate limits in AU" LOCATE 28, 5: PRINT USING " X: +/- ######.## "; LM * IAU LOCATE 29, 5: PRINT USING " Y: +/- ######.## "; LM * IAU LOCATE 30, 5: PRINT USING " Z: +/- ######.## "; LZ * IAU LINE (25, 390)-(240, 490), 15, B LOCATE 34, 5: COLOR 15: PRINT "KEY COMMANDS:" LOCATE 34, 20: PRINT "E/e: Expand X/Y limits of 10 AU" LOCATE 35, 20: PRINT "R/r: Shrink X/Y limits of 10 AU" LOCATE 34, 60: PRINT "Z/z: Expand Z limits of 10 AU" LOCATE 35, 60: PRINT "C/c: Shrink Z limits of 10 AU" LINE (25, 525)-(740, 560), 15, B LOCATE 26, 32: COLOR 15: PRINT "Parametri di Simulazione" LOCATE 28, 32: PRINT USING " TIME STEP (seconds): #####.##"; DT LOCATE 29, 32: PRINT USING " SIM. LENGHT (days) : ########"; SIML LOCATE 30, 32: PRINT USING " NUMBER OF STEPS : ########"; NSTEPS LINE (240, 390)-(500, 490), 15, B 'GOSUB 9000 REPLOT END SUB SUB RESCALE SHARED LM, AU, FX, FY, SIGN, AU 8000 ' ****** SUBROUTINE III ******* ' ****** RESCALE ****** LM = LM + SIGN * AU FX = 150 / (LM) FY = 150 / (LM) DISEGNA END SUB SUB REPLOT SHARED FX, FY, FZ, TY, N ' ****** SUBROUTINE IV ******** ' ***** REPLOT ***** FOR K = 1 TO N PSET (225 + FX * X1(K), 200 + FY * Y1(K)), K DRAW "D2U1R3L6" PSET (575 + FX * X1(K), 200 + FZ * Z1(K)), K DRAW "D2U1R3L6" NEXT K END SUB