Retro programming nostalgia III: the MSX Microcomputer and the Orbit of the Planets in the Solar System

In a recent article, I have explained the Euler’s method for solving ordinary differential equations using 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. As an example of application, I have also shown a program written in the awk programming language for calculating the orbits of planets of the solar system. However, 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 Planetario.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

About Danilo Roccatano

I have a Doctorate in chemistry at the University of Roma “La Sapienza”. I led educational and research activities at different universities in Italy, The Netherlands, Germany and now in the UK. I am fascinated by the study of nature with theoretical models and computational. For years, my scientific research is focused on the study of molecular systems of biological interest using the technique of Molecular Dynamics simulation. I have developed a server (the link is in one of my post) for statistical analysis at the amino acid level of the effect of random mutations induced by random mutagenesis methods. I am also very active in the didactic activity in physical chemistry, computational chemistry, and molecular modeling. I have several other interests and hobbies as video/photography, robotics, computer vision, electronics, programming, microscopy, entomology, recreational mathematics and computational linguistics.
This entry was posted in Programming, What is new. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.