10 REM
20 REM          The Apsidal Motion Program that Roared
30 REM
40 REM   Setting: Dept. of Astronomy and Astrophysics
50 REM                    Villanova University
60 REM                        Summer 1986              We bad!
70 REM                   as told to Sean Carroll
80 REM
90 REM                Revised, Summers 1995 and 1996
100 REM                  as told to James Marshall
110 REM
120 SCREEN 0
130 CLS
140 ON ERROR GOTO 2460
150 DIM l$(20), k$(20)
160 VORR = 2
170 STAR$ = "Your Star"
180 wrwk#(1) = -1#: wrwk#(2) = -1#
190 kgr# = 9.287200000000001D-03
200 k2#(1) = .006#: k2#(2) = .006#
210 i = 90
220 pi# = 3.1415926535#: ir = i * pi# / 180
230 G# = 9.903D+24
240 FOR l% = 1 TO 13
250 READ k$(l%)
260 NEXT l%
270 REM Introdcution screen
280 CLS
290 PRINT "Theoretical Apsidal Motion Calculator"
300 PRINT "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
310 PRINT
320 PRINT "This program computes theoretical apsidal motions for 2- and 3-star systems."
330 PRINT "Classical, general relativistic, and NGT apsidal motions are auotmatically"
340 PRINT "calculated, and there are options to iterate a third-body solution and"
350 PRINT "to compute formal errors.  You may enter parameters by hand, or use a"
360 PRINT ".apm input file from a previous run."
370 PRINT
380 PRINT "Enter the file name and path (no extension) of the .apm input file."
390 INPUT "Hit ENTER to input parameters by hand:  ", IFILE$
400 IF IFILE$ = "" THEN GOTO 810
410 REM open data file and read parameters from it
420 OPEN IFILE$ + ".apm" FOR INPUT AS #1
430 CLS
440 LINE INPUT #1, l$(0): STAR$ = LEFT$(l$(0), 9)
450 LINE INPUT #1, k$(1)
460 LINE INPUT #1, l$(1): p# = VAL(l$(1))
470 LINE INPUT #1, k$(2)
480 LINE INPUT #1, l$(2): m#(1) = VAL(l$(2))
490 LINE INPUT #1, k$(3)
500 LINE INPUT #1, l$(3): m#(2) = VAL(l$(3))
510 LINE INPUT #1, k$(4)
520 LINE INPUT #1, l$(4): k2#(1) = VAL(l$(4))
530 LINE INPUT #1, k$(5)
540 LINE INPUT #1, l$(5): k2#(2) = VAL(l$(5))
550 LINE INPUT #1, k$(6)
560 LINE INPUT #1, l$(6): e# = VAL(l$(6))
570 LINE INPUT #1, k$(7)
580 LINE INPUT #1, l$(7): r#(1) = VAL(l$(7))
590 LINE INPUT #1, k$(8)
600 LINE INPUT #1, l$(8): r#(2) = VAL(l$(8))
610 LINE INPUT #1, k$(11)
620 LINE INPUT #1, l$(11): i = VAL(l$(11)): ir = i * pi# / 180
630 LINE INPUT #1, l$(19): VORR = VAL(l$(19))
640 IF VORR = 1 THEN GOTO 700
650 LINE INPUT #1, k$(9)
660 LINE INPUT #1, l$(9): wrwk#(1) = VAL(l$(9))
670 LINE INPUT #1, k$(10)
680 LINE INPUT #1, l$(10): wrwk#(2) = VAL(l$(10))
690 GOTO 750
700 REM V sin(i)
710 LINE INPUT #1, k$(12)
720 LINE INPUT #1, l$(12): v1sini = VAL(l$(12))
730 LINE INPUT #1, k$(13)
740 LINE INPUT #1, l$(13): v2sini = VAL(l$(13))
750 CLOSE #1
760 IF VORR = 1 THEN
770 GOTO 1300
780 ELSE GOTO 1450
790 END IF
800 REM section to enter parameters manually
810 CLS
820 PRINT "Data Entry Section"
830 PRINT "~~~~~~~~~~~~~~~~~~"
840 PRINT
850 PRINT "(type ENTER for default values [in brackets])"
860 PRINT "Enter name of star"
870 PRINT "             ["; STAR$; : INPUT "]  ", i$
880 IF i$ <> "" THEN STAR$ = i$
890 PRINT "Enter P, the period of revolution in days"
900 PRINT "             ["; p#; : INPUT "]  ", i$
910 IF i$ <> "" THEN p# = VAL(i$)
920 PRINT "Enter M, the mass in solar masses"
930 PRINT "     primary ["; m#(1); : INPUT "]  ", i$
940 IF i$ <> "" THEN m#(1) = VAL(i$)
950 PRINT "   secondary ["; m#(2); : INPUT "]  ", i$
960 IF i$ <> "" THEN m#(2) = VAL(i$)
970 PRINT "Enter k2, the internal concentration parameter"
980 PRINT "     primary ["; k2#(1); : INPUT "]  ", i$
990 IF i$ <> "" THEN k2#(1) = VAL(i$)
1000 PRINT "   secondary ["; k2#(2); : INPUT "]  ", i$
1010 IF i$ <> "" THEN k2#(2) = VAL(i$)
1020 PRINT "Enter e, the eccentricity of the orbit"
1030 PRINT "             ["; e#; : INPUT "]  ", i$
1040 IF i$ <> "" THEN e# = VAL(i$)
1050 IF wrwk#(1) < 0 THEN wrwk#(1) = 1# + .5 * (12# * e# * e#) ^ (3# / 8#)
1060 IF wrwk#(2) < 0 THEN wrwk#(2) = 1# + .5 * (12# * e# * e#) ^ (3# / 8#)
1070 PRINT "Enter r, the fractional radius (R/A)"
1080 PRINT "     primary ["; r#(1); : INPUT "]  ", i$
1090 IF i$ <> "" THEN r#(1) = VAL(i$)
1100 PRINT "   secondary ["; r#(2); : INPUT "]  ", i$
1110 IF i$ <> "" THEN r#(2) = VAL(i$)
1120 PRINT "Enter i, the inclination of the orbit"
1130 PRINT "             ["; i; : INPUT "]  ", i$
1140 IF i$ <> "" THEN i = VAL(i$): ir = i * pi# / 180
1150 PRINT
1160 PRINT "Do you want to use: 1) v"; CHR$(250); "sini or 2) the ratio of angular rotation speed"
1170 PRINT "to mean angular orbital velocity? (type 1 or 2)"
1180 PRINT "           ["; VORR; : INPUT "]  ", i$
1190 IF i$ <> "" THEN VORR = VAL(i$)
1200 PRINT
1210 IF VORR = 1 THEN GOTO 1240
1220 IF VORR = 2 THEN GOTO 1380
1230 GOTO 1150: REM only continue if 1 or 2 was entered
1240 REM V sin(i)
1250 PRINT "Enter v"; CHR$(250); "sini"
1260 PRINT "     primary ["; v1sini; : INPUT "]  ", i$
1270 IF i$ <> "" THEN v1sini = VAL(i$)
1280 PRINT "   secondary ["; v2sini; : INPUT "]  ", i$
1290 IF i$ <> "" THEN v2sini = VAL(i$)
1300 v1 = v1sini / SIN(ir): v2 = v2sini / SIN(ir)
1310 a# = 149600000 * ((m#(1) + m#(2)) * ((p# / 365.25) ^ 2)) ^ (1 / 3): REM semi-major axis
1330 br1 = r#(1) * a#: br2 = r#(2) * a#: REM true stellar radii
1340 C1 = 2 * pi# * a#
1350 t1 = C1 / (p# * 24 * 3600)
1360 wrwk#(1) = (v1 / br1) / (t1 / a#): wrwk#(2) = (v2 / br2) / (t1 / a#)
1370 GOTO 1450
1380 REM Wr/Wk
1390 PRINT "Enter Wr/Wk, the ratio of angular axial rotation speed of the star to the "
1400 PRINT "mean orbital angular Keplerian velocity"
1410 PRINT "     primary ["; wrwk#(1); : INPUT "]  ", i$
1420 IF i$ <> "" THEN wrwk#(1) = VAL(i$)
1430 PRINT "   secondary ["; wrwk#(2); : INPUT "]  ", i$
1440 IF i$ <> "" THEN wrwk#(2) = VAL(i$)
1450 REM Output screen
1460 CLS
1470 PRINT k$(1); " "; p#
1480 PRINT k$(2); " "; m#(1)
1490 PRINT k$(3); " "; m#(2)
1500 PRINT k$(4); " "; k2#(1)
1510 PRINT k$(5); " "; k2#(2)
1520 PRINT k$(6); " "; e#
1530 PRINT k$(7); " "; r#(1)
1540 PRINT k$(8); " "; r#(2)
1550 PRINT k$(11); " "; i
1560 IF VORR = 2 THEN GOTO 1600
1570 PRINT k$(12); " "; v1sini
1580 PRINT k$(13); " "; v2sini
1590 GOTO 1620
1600 PRINT k$(9); " "; wrwk#(1)
1610 PRINT k$(10); " "; wrwk#(2)
1620 REM Classical apsidal motion
1630 f2# = (1 + 1.5 * e# * e# + .125 * e# ^ 3) / ((1 - e# * e#) ^ 5)
1640 wcl# = k2#(1) * r#(1) ^ 5 * (15# * f2# * m#(2) / m#(1) + wrwk#(1) * wrwk#(1) * ((1 + m#(2) / m#(1)) / (1 - e# * e#) ^ 2))
1650 wcl# = wcl# + k2#(2) * r#(2) ^ 5 * (15 * f2# * m#(1) / m#(2) + wrwk#(2) * wrwk#(2) * ((1 + m#(1) / m#(2)) / (1 - e# * e#) ^ 2))
1660 wcl# = (365.25 * 360 / p#) * wcl#
1670 PRINT
1680 PRINT USING "* * * * * *  \        \  Wcl =     ####.### degrees/100 years     * * * * * *"; STAR$; wcl# * 100
1690 REM General Relativistic (Albert's) apsidal motion
1700 wgr# = kgr# * (m#(1) + m#(2)) ^ (2 / 3) / (p# / 2 / pi#) ^ (5 / 3) / (1 - e# * e#)
1710 PRINT USING "* * * * * *              Wgr =     ####.### degrees/100 years     * * * * * *"; wgr# * 100
1720 PRINT USING "* * * * * *              Wcl+gr =  ####.### degrees/100 years     * * * * * *"; (wgr# + wcl#) * 100
1730 REM NGT (Non-symmetric Gravitational Theory) apsidal motion
1740 l# = (.2 * (m#(1) + m#(2)) + .08) * 1E+09
1750 L1# = (.2 * m#(1) + .08) * 1E+09
1760 L2# = (.2 * m#(2) + .08) * 1E+09
1770 LP# = (L1# ^ 2 + L2# ^ 2) ^ .5
1780 LAMBDA# = 1 - (4.626E-35) * l# * l# * l# * l# * (1 + e# * e# / 4#) / ((p# / (2 * pi#)) ^ (4 / 3) * (1 - e# * e#) ^ 2 * (m#(1) + m#(2)) ^ (8 / 3))
1790 wngt# = wgr# * LAMBDA#
1800 PRINT USING "* * * * * *              Wngt =    ####.### degrees/100 years     * * * * * *"; wngt# * 100
1810 PRINT USING "* * * * * *              Wcl+ngt = ####.### degrees/100 years     * * * * * *"; (wcl# + wngt#) * 100
1820 REM Menu to choose what to do next
1830 PRINT
1840 PRINT "Choose: 1) change parameters               4) iterate third-body solutions"
1850 PRINT "        2) save parameters and results     5) compute formal errors"
1860 PRINT "        3) start a new star                6) quit"
1870 PRINT
1880 INPUT "Enter your choice: ", yn
1890 IF yn = 1 THEN GOTO 810: REM branch back to data entry section
1900 IF yn = 2 THEN GOTO 1970: REM branch to save file section
1910 IF yn = 3 THEN CLEAR : GOTO 10: REM branch all the way back to the top
1920 IF yn = 4 THEN GOSUB 2500: REM branch to 3rd-body section
1930 IF yn = 5 THEN GOTO 3350: REM banch to formal errors section
1940 IF yn = 6 THEN GOTO 4580: REM end it
1950 GOTO 1820: REM go back to menu if choice is not 1-6
1960 REM Create a data file for output of parameters and results
1970 INPUT "Enter the filename with path (.apm automatic):  ", OFILE$
1980 OPEN OFILE$ + ".apm" FOR OUTPUT AS #2
1990 REM parameters
2000 PRINT #2, STAR$
2010 PRINT #2, k$(1)
2020 PRINT #2, p#; "days"
2030 PRINT #2, k$(2)
2040 PRINT #2, m#(1); "solar masses"
2050 PRINT #2, k$(3)
2060 PRINT #2, m#(2); "solar masses"
2070 PRINT #2, k$(4)
2080 PRINT #2, k2#(1)
2090 PRINT #2, k$(5)
2100 PRINT #2, k2#(2)
2110 PRINT #2, k$(6)
2120 PRINT #2, e#
2130 PRINT #2, k$(7)
2140 PRINT #2, r#(1)
2150 PRINT #2, k$(8)
2160 PRINT #2, r#(2)
2170 PRINT #2, k$(11)
2180 PRINT #2, i
2190 IF VORR = 2 THEN GOTO 2260
2200 PRINT #2, "1 -- We will use Vsin(i), not Wr/Wk"
2210 PRINT #2, k$(12)
2220 PRINT #2, v1sini
2230 PRINT #2, k$(13)
2240 PRINT #2, v2sini
2250 GOTO 2320
2260 PRINT #2, "2 -- We will use Wr/Wk, not Vsin(i)"
2270 PRINT #2, k$(9)
2280 PRINT #2, wrwk#(1)
2290 PRINT #2, k$(10)
2300 REM results
2310 PRINT #2, wrwk#(2)
2320 PRINT #2, "Wcl"
2330 PRINT #2, (wcl# * 100)
2340 PRINT #2, "Wgr"
2350 PRINT #2, (wgr# * 100)
2360 PRINT #2, "Wcl+gr"
2370 PRINT #2, ((wcl# + wgr#) * 100)
2380 PRINT #2, "Wngt"
2390 PRINT #2, (wngt# * 100)
2400 PRINT #2, "Wcl+ngt"
2410 PRINT #2, ((wcl# + wngt#) * 100)
2420 CLOSE #2
2430 PRINT "File "; F$; ".apm completed"
2440 GOTO 1820: REM go back to menu
2450 END
2460 REM Error trapping routine
2470 IF ERR = 53 THEN PRINT "File not found -- try again": RESUME 580
2480 PRINT "Error number"; ERR; "on line"; ERL
2490 STOP
2500 REM Third-body routine
2510 IPMIN = 0: IPMAX = 90: DIP = 5
2520 M3MIN = .5: M3MAX = 3: DM3 = .5
2530 PPMIN = 1: PPMAX = 10: DPP = 1
2540 EPMIN = 0: EPMAX = .001: DEP = .1
2550 CLS
2560 PRINT "Third-Body Calculations"
2570 PRINT "~~~~~~~~~~~~~~~~~~~~~~~"
2580 PRINT
2590 PRINT "We are about to explore a parameter-space describing possible"
2600 PRINT "orbits and masses of a third body around "; STAR$; "."
2610 PRINT "Your job is to enter the upper and lower bounds of these values,"
2620 PRINT "as well as the value we will increment by in our exploration."
2630 PRINT
2640 REM get 3rd-body parameters
2650 PRINT "Enter i', the inclination of the third body orbit"
2660 PRINT "    minimum ["; IPMIN; : INPUT "]  ", i$
2670 IF i$ <> "" THEN IPMIN = VAL(i$)
2680 PRINT "    maximum ["; IPMAX; : INPUT "]  ", i$
2690 IF i$ <> "" THEN IPMAX = VAL(i$)
2700 PRINT "  increment ["; DIP; : INPUT "]  ", i$
2710 IF i$ <> "" THEN DIP = VAL(i$)
2720 PRINT "Enter M3, the mass of the third body"
2730 PRINT "    minimum ["; M3MIN; : INPUT "]  ", i$
2740 IF i$ <> "" THEN M3MIN = VAL(i$)
2750 PRINT "    maximum ["; M3MAX; : INPUT "]  ", i$
2760 IF i$ <> "" THEN M3MAX = VAL(i$)
2770 PRINT "  increment ["; DM3; : INPUT "]  ", i$
2780 IF i$ <> "" THEN DM3 = VAL(i$)
2790 PRINT "Enter P', the period of the 3rd body (yrs)"
2800 PRINT "    minimum ["; PPMIN; : INPUT "]  ", i$
2810 IF i$ <> "" THEN PPMIN = VAL(i$)
2820 PRINT "    maximum ["; PPMAX; : INPUT "]  ", i$
2830 IF i$ <> "" THEN PPMAX = VAL(i$)
2840 PRINT "  increment ["; DPP; : INPUT "]  ", i$
2850 IF i$ <> "" THEN DPP = VAL(i$)
2860 PRINT "Enter e', the eccentricity of the third body orbit"
2870 PRINT "    minimum ["; EPMIN; : INPUT "]  ", i$
2880 IF i$ <> "" THEN EPMIN = VAL(i$)
2890 PRINT "    maximum ["; EPMAX; : INPUT "]  ", i$
2900 IF i$ <> "" THEN EPMAX = VAL(i$)
2910 PRINT "  increment ["; DEP; : INPUT "]  ", i$
2920 IF i$ <> "" THEN DEP = VAL(i$)
2930 PRINT
2940 PRINT "The results will be entered into an output file with extension .apm"
2950 IF OFILE$ <> "" THEN F$ = OFILE$ ELSE F$ = IFILE$
2960 PRINT "Enter the name of your data file ["; F$; : INPUT "]  ", i$
2970 IF i$ <> "" THEN F$ = i$
2980 OPEN F$ + ".apm" FOR APPEND AS #3: REM append results to end of file
2990 REM calculate it
3000 PRINT : PRINT "Running . . ."
3010 PS# = p# / 365.25#: ns# = 2 * pi# / PS#
3020 FOR EP = EPMIN TO EPMAX STEP DEP
3030 FOR PP = PPMIN TO PPMAX STEP DPP
3040 FOR M3 = M3MIN TO M3MAX STEP DM3
3050 MSUM# = m#(1) + m#(2) + M3
3060 PRINT #3, ""
3070 PRINT #3, "  Eccentricity ="; EP; "       Period ="; PP; "       Mass ="; M3
3080 PRINT #3, "    Inclination    Omega(100yrs)   LTC(days)"
3090 FOR IP = IPMIN TO IPMAX STEP DIP
3100 GAMMA# = SIN(IP * pi# / 360)
3110 MS# = PS# / PP
3120 LLS# = M3 / MSUM#
3130 ZETA# = 1 - EP ^ 2
3140 C1# = (3 / 8) * MS# ^ 2 * ns# * LLS# / ZETA# ^ 1.5
3150 C2# = (9 / 64) * MS# ^ 3 * ns# * LLS# ^ 2 * (1 + 2 * EP * EP / 3) / ZETA# ^ 3
3160 W3B# = C1# * (2 - e# * e# - 16 * GAMMA# ^ 2 - .25 * e# ^ 4 + 23 * (e# * GAMMA#) ^ 2 - 30 * GAMMA# ^ 4 - .25 * e# ^ 6)
3170 W3B# = W3B# + C1# * (-37 / 4 * e# ^ 4 * GAMMA# ^ 2 + 345 / 2 * e# ^ 2 * GAMMA# ^ 4 - 150 * GAMMA# ^ 6)
3180 W3B# = W3B# + C2# * (50 - 75 * e# ^ 2 - 168 * GAMMA# ^ 2 + 18 * e# ^ 2 * GAMMA# ^ 2 + 492 * GAMMA# ^ 4)
3190 W3B# = W3B# + C2# * (1425 / 4 * e# ^ 4 * GAMMA# ^ 2 - 7179 / 2 * e# ^ 2 * GAMMA# ^ 4 + 2430 * GAMMA# ^ 6)
3200 W3B# = W3B# * 57.29577951# * 100
3210 I3# = (i - IP) * pi# / 180#
3220 ROOT# = (PP ^ 2 * MSUM#) ^ (1# / 3#)
3230 LTC# = (LLS# * 498.37071# * ROOT# * SIN(I3#)) / 86400#
3240 LOCATE 24, 1
3250 PRINT USING "e'= #.###   P'= ###.##   M3= ###.##   i'= ####.#   W3b= ####.###### "; EP; PP; M3; IP; W3B#;
3260 PRINT #3, USING "     ####.#        ####.#####      ##.#######"; IP; W3B#; LTC#
3270 NEXT IP
3280 NEXT M3
3290 NEXT PP
3300 NEXT EP
3310 CLOSE #3
3320 PRINT : PRINT "Done!"
3330 PRINT "Calculations complete and file written."
3340 GOTO 1820: REM go back to menu
3350 REM Compute formal errors
3360 dwgrde# = ABS(kgr# * (m#(1) + m#(2)) ^ (2 / 3) * (2 * pi# / p#) ^ (5 / 3) * (1 - 2 * e#) / (1 - e# * e#) ^ 2)
3370 dwgrdp# = ABS(kgr# * (m#(1) + m#(2)) ^ (2 / 3) / (1 - e# * e#) * (2 * pi#) ^ (5 / 3) / p# ^ (8 / 3))
3380 dwgrdm# = ABS(kgr# * (m#(1) + m#(2)) ^ (-1 / 3) * (2 * pi# / p#) ^ (5 / 3) / (1 - e# * e#))
3390 f2# = (1 + 1.5 * e# * e# + .125 * e# ^ 4) / ((1 - e# * e#) ^ 5)
3400 f3# = ((1 - e# * e#) ^ 5 * (1 + 3 * e# + e# ^ 3 / 2) - (1 + 1.5 * e# * e# + .125 * e# ^ 4) * 5 * (1 - e# * e#) ^ 4 * (-2 * e#)) / (1 - e# * e#) ^ 10
3410 alpha#(1) = k2#(1) * r#(1) ^ 5 * (15# * f2# * m#(2) / m#(1) + wrwk#(1) * wrwk#(1) * ((1 + m#(2) / m#(1)) / (1 - e# * e#) ^ 2))
3420 alpha#(2) = k2#(2) * r#(2) ^ 5 * (15# * f2# * m#(1) / m#(2) + wrwk#(2) * wrwk#(2) * ((1 + m#(1) / m#(2)) / (1 - e# * e#) ^ 2))
3430 cons# = 365.25# * 360 / p#
3440 dwcldp# = cons# / p# * (alpha#(1) + alpha#(2))
3450 FOR i = 1 TO 2
3460 j = ABS(i - 3)
3470 dadk2#(i) = alpha#(i) / k2#(i)
3480 dadr#(i) = alpha#(i) * 5 / r#(i)
3490 dade#(i) = k2#(i) * r#(i) ^ 5 * (15 * (m#(j) / m#(i)) * f3# + wrwk#(i) ^ 2 * (1 + m#(j) / m#(i)) * (-2) / (1 - e# * e#) ^ 3 * (-2 * e#))
3500 dadmi#(i) = k2#(i) * r#(i) ^ 5 * (15 * f2# * (-m#(j) / m#(i) ^ 2) + wrwk#(i) ^ 2 * (1 - e# * e#) ^ (-2) * (-m#(j) / m#(i) ^ 2))
3510 dadmj#(i) = k2#(i) * r#(i) ^ 5 * (15 * f2# / m#(i) + wrwk#(i) ^ 2 * (1 - e# * e#) ^ (-2) * (-1 / m#(i)))
3520 dadw#(i) = k2#(i) * r#(i) ^ 5 * (2 * wrwk#(i) * (1 - e# * e#) ^ (-2) * (1 + m#(j) / m#(i)))
3530 NEXT i
3540 dwcldk2# = ABS(cons# * (dadk2#(1) + dadk2#(2)))
3550 dwcldr# = ABS(cons# * (dadr#(1) + dadr#(2)))
3560 dwclde# = ABS(cons# * (dade#(1) + dade#(2)))
3570 dwcldmi# = ABS(cons# * (dadmi#(1) + dadmi#(2)))
3580 dwcldmj# = ABS(cons# * (dadmj#(1) + dadmj#(2)))
3590 dwcldw# = ABS(cons# * (dadw#(1) + dadw#(2)))
3600 REM ask for inputs and finish computations
3610 CLS
3620 PRINT "Calculation of Formal Errors"
3630 PRINT "~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
3640 PRINT
3650 PRINT "Input the following values"
3660 PRINT
3670 INPUT "      Delta-e:  ", dele
3680 INPUT "      Delta-p:  ", delp
3690 INPUT "      Delta-m:  ", delm
3700 INPUT "     Delta-k2:  ", delk2
3710 INPUT "      Delta-r:  ", delr
3720 INPUT "  Delta-wr/wk:  ", delw
3730 PRINT : PRINT
3740 PRINT "Formal errors (classical apsidal motion):   "
3750 PRINT USING "    e: ##.####"; dele * dwclde# * 100
3760 PRINT USING "    p: ##.####"; delp * dwcldp# * 100
3770 PRINT USING "    m: ##.####"; SQR((delm * dwcldmi#) ^ 2 + (delm * dwcldmj#) ^ 2 + (delm * dwcldmi#) ^ 2 + (delm * dwcldmj#) ^ 2) * 100
3780 PRINT USING "   k2: ##.####"; delk2 * dwcldk2# * 100
3790 PRINT USING "    r: ##.####"; delr * dwcldr# * 100
3800 PRINT USING "wr/wk: ##.####"; delw * dwcldw# * 100
3810 PRINT
3820 PRINT USING "Apsidal Motion: #####.####"; wcl# * 100;
3830 PRINT " deg/100yr"
3840 term1# = (dele * dwclde#) ^ 2 + (delp * dwcldp#) ^ 2 + (delm * dwcldmi#) ^ 2 + (delm * dwcldmj#) ^ 2 + (delm * dwcldmi#) ^ 2
3850 term2# = (delm * dwcldmj#) ^ 2 + (delk2 * dwcldk2#) ^ 2 + (delr * dwcldr#) ^ 2 + (delw * dwcldw#) ^ 2
3860 answer# = SQR(term1# + term2#) * 100
3870 PRINT USING "Total Errors:   #####.####"; answer#
3880 PRINT
3890 PRINT
3900 PRINT "Formal errors (relativistic apsidal motion):   "
3910 PRINT USING "    e: ##.####"; dele * dwgrde# * 100
3920 PRINT USING "    p: ##.####"; delp * dwgrdp# * 100
3930 PRINT USING "    m: ##.####"; SQR((delm * dwgrdm#) ^ 2 + (delm * dwgrdm#) ^ 2) * 100
3940 PRINT
3950 PRINT USING "Apsidal Motion: #####.####"; wgr# * 100;
3960 PRINT " deg/100yr"
3970 PRINT USING "Total Errors:   #####.####"; SQR((dele * dwgrde#) ^ 2 + (delp * dwgrdp#) ^ 2 + (delm * dwgrdm#) ^ 2 + (delm * dwgrdm#) ^ 2) * 100
3980 PRINT
3990 PRINT "Hit any key to continue. "
4000 WHILE INKEY$ = "": WEND
4010 REM optional save (append) to .apm file
4020 PRINT : INPUT "Do you want to save these results into a file (Y/N)"; yn$
4030 IF (yn$ = "n" OR yn$ = "N") THEN GOTO 4430: REM skip save routine
4040 IF (yn$ <> "y" AND yn$ <> "Y") THEN GOTO 4020: REM must be yes or no
4050 PRINT "The results will be entered into an output file with extension .apm"
4060 IF OFILE$ <> "" THEN F$ = OFILE$ ELSE F$ = IFILE$
4070 PRINT "Enter the name of your data file ["; F$; : INPUT "]  ", i$
4080 IF i$ <> "" THEN F$ = i$
4090 REM
4100 OPEN IFILE$ + ".apm" FOR APPEND AS #4: REM append results to end of file
4110 PRINT #4, ""
4120 PRINT #4, "Calculation of Formal Errors"
4130 PRINT #4, "      Delta-e:  ", dele
4140 PRINT #4, "      Delta-p:  ", delp
4150 PRINT #4, "      Delta-m:  ", delm
4160 PRINT #4, "     Delta-k2:  ", delk2
4170 PRINT #4, "      Delta-r:  ", delr
4180 PRINT #4, "  Delta-wr/wk:  ", delw
4190 PRINT #4, "": PRINT #4, ""
4200 PRINT #4, "Formal errors (classical apsidal motion):   "
4210 PRINT #4, USING "    e: ##.####"; dele * dwclde# * 100
4220 PRINT #4, USING "    p: ##.####"; delp * dwcldp# * 100
4230 PRINT #4, USING "    m: ##.####"; SQR((delm * dwcldmi#) ^ 2 + (delm * dwcldmj#) ^ 2 + (delm * dwcldmi#) ^ 2 + (delm * dwcldmj#) ^ 2) * 100
4240 PRINT #4, USING "   k2: ##.####"; delk2 * dwcldk2# * 100
4250 PRINT #4, USING "    r: ##.####"; delr * dwcldr# * 100
4260 PRINT #4, USING "wr/wk: ##.####"; delw * dwcldw# * 100
4270 PRINT #4, ""
4280 PRINT #4, USING "Apsidal Motion: #####.####"; wcl# * 100;
4290 PRINT #4, " deg/100yr"
4300 PRINT #4, USING "Total Errors:   #####.####"; SQR((dele * dwclde#) ^ 2 + (delp * dwcldp#) ^ 2 + (delm * dwcldmi#) ^ 2 + (delm * dwcldmj#) ^ 2 + (delm * dwcldmi#) ^ 2 + (delm * dwcldmj#) ^ 2 + (delk2 * dwcldk2#) ^ 2 + (delr * dwcldr#) ^ 2 + ( _
4310 * dwcldw#) ^ 2) * 100
4320 PRINT #4, ""
4330 PRINT #4, ""
4340 PRINT #4, "Formal errors (relativistic apsidal motion):   "
4350 PRINT #4, USING "    e: ##.####"; dele * dwgrde# * 100
4360 PRINT #4, USING "    p: ##.####"; delp * dwgrdp# * 100
4370 PRINT #4, USING "    m: ##.####"; SQR((delm * dwgrdm#) ^ 2 + (delm * dwgrdm#) ^ 2) * 100
4380 PRINT #4, ""
4390 PRINT #4, USING "Apsidal Motion: #####.####"; wgr# * 100;
4400 PRINT #4, " deg/100yr"
4410 PRINT #4, USING "Total Errors:   #####.####"; SQR((dele * dwgrde#) ^ 2 + (delp * dwgrdp#) ^ 2 + (delm * dwgrdm#) ^ 2 + (delm * dwgrdm#) ^ 2) * 100
4420 PRINT "Results saved in "; IFILE$; ".apm"
4430 GOTO 1820: REM go back to menu
4440 REM the labels for the parameters
4450 DATA "Period:"
4460 DATA "Mass of the primary:"
4470 DATA "Mass of the secondary:"
4480 DATA "Internal density parameter of primary:"
4490 DATA "Internal density parameter of secondary:"
4500 DATA "Eccentricity:"
4510 DATA "Fractional radius of the primary:"
4520 DATA "Fractional radius of the secondary:"
4530 DATA "Wr/Wk for the primary:"
4540 DATA "Wr/Wk for the secondary:"
4550 DATA "Inclination:"
4560 DATA "V sin(i) for the primary:"
4570 DATA "V sin(i) for the secondary:"
4580 PRINT
4590 PRINT "Glad to be of service.  Give my regards to Albert!"
4600 PRINT
4610 PRINT "Have a nice day!  :)"
4620 END