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 jealously 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 Dewdney, 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. I have run the cassette on a Walkman tape recorder connecting the audio output to the mic input of my laptop. The audio signal was recorded and stored in wav format using the freeware software audacity. The MSX as most of the microcomputers from the ’80, stored the information on audio tapes as square waves at 1200 or 2400 baud. The encoding used was the frequency shifting key, or FSK. This is a frequency modulation scheme in which digital information is transmitted through discrete frequency changes of a carrier signal. Nowadays, this method is is commonly used over telephone lines for displaying the caller number of for other digital information.
To retrieve these programs, I have used the open-source MSX emulator OpenMSX (https://openmsx.org). However, using this emulator (not to mention using the old hardware) to improve the program was 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
daniloroccatano.blog 