REM Binary Star Period Study and Apsidal Motion Rate Calculator
REM By James Marshall (class of 1996), summer 1996, Villanova University
REM Last updated on June 6, 1996.
REM
REM This program will take inputs for a system's T_min and period of primary
REM and secondary eclipses, the number of T_min I's and II's, the eccentricty
REM of the orbit, the longitude of periastron, and a set of observations of
REM eclipse timings and make calculations of the O-C's, epochs/cycle numbers,
REM the D's of secondary eclipse, a new ephemeris for primary and secondary
REM eclipse and an observed apsidal motion rate.  It takes into account a
REM relative weight for the data -- this is usually given as a number from
REM one to ten with a ten indicating the best data and a one for the worst.
REM
REM *** Some constants ***
REM
CONST PI = 3.14159265359#
CONST DEGTORAD = PI / 180
CONST RADTODEG = 180 / PI
CONST DAYSPERYEAR = 365.242198781#
wrotein1 = 0: REM assume data comes from file; if not, flag gets set to 1
wrotein2 = 0: REM assume data comes from file; if not, flag gets set to 1
REM
REM *** Instructions and Stuff ***
REM
CLS
PRINT "Welcome to the Binary Star Period Study and Apsidal Motion Rate Calculator!"
PRINT "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
PRINT
PRINT "This program will enable you to analyze the elcipse timings of a binary"
PRINT "star system.  Knowing an intial ephemeris (T_min and period) for both the"
PRINT "primary and secondary eclipses, as well as the eccentricity of the orbit"
PRINT "and the longitude of periastron for the orbit, the program takes a set of"
PRINT "eclipse timings and their relative weights and then calculates the O-C's,"
PRINT "epochs/cycle numbers, D for secondary eclipse [ D = (t2 - t1) - .5 * P ],"
PRINT "a new updated ephemeris for the input data, and the observed apsidal"
PRINT "motion rate for the system.  If you do not have a separate T_min for the"
PRINT "secondary eclipse, use the primary's T_min + 0.5*P and the same period."
PRINT
PRINT "It is suggested that you use input files; however, you may choose to input"
PRINT "the necessary values directly into the program.  If you choose the latter"
PRINT "option, the program will automatically write the data into the proper input"
PRINT "files so that you can use them again later without typing it all in again."
PRINT
PRINT "A note on units:  the normal units are Julian Date for all T_min's, days"
PRINT "for periods, O-C's, D's, eccentricity is unitless, the longitude of"
PRINT "periastron is in degrees, and apsidal motion is in degrees per century."
PRINT
INPUT "Do you require further instructions (Y/N)"; ans$
IF (ans$ <> "Y" AND ans$ <> "y") THEN GOTO skipinst
CLS
PRINT "Instructions on Files"
PRINT "~~~~~~~~~~~~~~~~~~~~~"
PRINT
PRINT "In order for the program to correctly read in the required data from an"
PRINT "input file, you must create the files in a specific format.  The first data"
PRINT "file must be of the following form and it must have a .in1 extension:"
PRINT
PRINT "Star Name"
PRINT "T_min (I)"
PRINT "Period (I)"
PRINT "Number of T_min I's"
PRINT "T_min (II)"
PRINT "Period (II)"
PRINT "Number of T_min II's"
PRINT "eccentricity of the orbit"
PRINT "longitude of periastron (in degrees)"
PRINT "Name of person(s) who determined the data used above"
PRINT
PRINT "Please note that this requires only one item per line and it must be only"
PRINT "a number, or a string in the case of the two names requested."
PRINT
PRINT "<Press any key to continue.>"
DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
CLS
PRINT "More on Files"
PRINT "~~~~~~~~~~~~~"
PRINT
PRINT "The second input file is a little easier to do.  This requires placing"
PRINT "all of your eclipse timings and their relative weights in a file, listing"
PRINT "the primary minima first, followed by the secondary minima.  The relative"
PRINT "weight is an indicator of how accurate the timing is and it tells the"
PRINT "program how many times it should count that particular eclipse timing"
PRINT "when doing the least squares fit to the data to calculate a new ephemeris."
PRINT "A general method is to assign relative weights from 1 to 10, with 10 being"
PRINT "the best data.  For equally weighted data, just use 1 for all timings."
PRINT "The file must have a .in2 extension and it should look somthing like this:"
PRINT
PRINT "T_min (I)"
PRINT "Rel. Wt."
PRINT " . . . "
PRINT "T_min (II)"
PRINT "Rel. Wt."
PRINT " . . . "
PRINT
PRINT "<Press any key to continue.>"
DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
CLS
PRINT "Entering Data by Hand and More Info"
PRINT "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
PRINT
PRINT "In the case of both of these input files, you may choose to input the"
PRINT "data by hand instead.  This is especially good for the first time you"
PRINT "run the program, since the program will automatically write all of the"
PRINT "data you input into the appropriate files for you.  You can then use"
PRINT "the computer-written input files for any future runs."
PRINT
PRINT "All of the calculations will be shown on the screen and written to files."
PRINT "The files will all have the star's name as it's filename, but with"
PRINT "different extension, depending on the data it contains -- .o-c is for"
PRINT "the O-C calculations, .eph is for the new ephemeris, and .cal is for the"
PRINT "calculation of the apsidal motion rate.  The program will remind you of"
PRINT "the file names and extensions later in the program as well."
PRINT
PRINT "<Press any key to begin the program.>"
DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
skipinst:
REM
REM *** Initial Inputs of Data ***
CLS
PRINT "Input Data for the System"
PRINT "~~~~~~~~~~~~~~~~~~~~~~~~~"
PRINT
PRINT "Enter the filename and path of the .in1 input file (or ENTER for none):"
INPUT "   "; inpfile1$
IF inpfile1$ <> "" THEN
 OPEN inpfile1$ + ".in1" FOR INPUT AS #1
 LINE INPUT #1, starname$
 LINE INPUT #1, tmin1$
 tmin1# = VAL(tmin1$)
 LINE INPUT #1, p1$
 p1# = VAL(p1$)
 LINE INPUT #1, numtmin1$
 numtmin1% = VAL(numtmin1$)
 LINE INPUT #1, tmin2$
 tmin2# = VAL(tmin2$)
 LINE INPUT #1, p2$
 p2# = VAL(p2$)
 LINE INPUT #1, numtmin2$
 numtmin2% = VAL(numtmin2$)
 LINE INPUT #1, e$
 e# = VAL(e$)
 LINE INPUT #1, wdeg$
 wdeg# = VAL(wdeg$)
 LINE INPUT #1, ephem$
 CLOSE #1
 REM
 REM *** Define Output File Names ***
 REM
 outfile$ = inpfile1$
ELSE
 REM *** Enter Data By Hand Instead ***
entry1:
 CLS
 PRINT "Input Data by Hand"
 PRINT "~~~~~~~~~~~~~~~~~~"
 PRINT
 INPUT "           Star Name:  ", starname$
 INPUT "           T_min (I):  ", tmin1#
 INPUT "               P (I):  ", p1#
 INPUT "      # of T_min I's:  ", numtmin1%
 INPUT "          T_min (II):  ", tmin2#
 INPUT "              P (II):  ", p2#
 INPUT "     # of T_min II's:  ", numtmin2%
 INPUT "      Eccentricty, e:  ", e#
 INPUT " Long. peri., w(deg):  ", wdeg#
 PRINT
 INPUT "Data from:  ", ephem$
 PRINT
 INPUT "Are all the values correct (Y/N)"; ans$
 IF (ans$ = "N" OR ans$ = "n") THEN GOTO entry1
 REM
 REM *** Define Output File Names ***
 REM
 star$ = starname$
 FOR i = 1 TO LEN(star$)
  letter$ = LEFT$(star$, 1)
  IF letter$ <> "_" THEN out$ = out$ + letter$
  star$ = MID$(star$, 2)
 NEXT i
 IF LEN(out$) > 8 THEN
  outfile$ = LEFT$(out$, 8)
 ELSE
  outfile$ = out$
 END IF
 REM
 REM *** Write Data to a File for Future Use (.in1) ***
 REM
 wrotein1 = 1
 OPEN outfile$ + ".in1" FOR OUTPUT AS #1
 PRINT #1, starname$
 PRINT #1, tmin1#
 PRINT #1, p1#
 PRINT #1, numtmin1%
 PRINT #1, tmin2#
 PRINT #1, p2#
 PRINT #1, numtmin2%
 PRINT #1, e#
 PRINT #1, wdeg#
 PRINT #1, ephem$
 CLOSE #1
END IF
PRINT
IF wrotein1 = 1 THEN CLS
PRINT "Please note that all calculation will be output to files with the filename"
PRINT outfile$; " and the following extensions:"
PRINT
PRINT "  .o-c for eclipse times, O-C's, epochs/cycle numbers, and D's"
PRINT "  .eph for the new T_min and periods of primary and secondary eclipses"
PRINT "  .cal for the calculate apsidal motion rate"
PRINT
PRINT "You will see this list again when the program is finished calcuating."
PRINT
IF wrotein1 = 1 THEN PRINT "The files will be located in the current working directory."
PRINT
PRINT "<Press any key to continue.>"
DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
REM
REM *** Get the T_min Eclipse Timings ***
REM
CLS
DIM tminprim#(numtmin1%), tminsec#(numtmin2%), e1#(numtmin1%), e2#(numtmin2%)
DIM phi1#(numtmin1%), phi2#(numtmin2%), oc1#(numtmin1%), oc2#(numtmin2%)
DIM e3#(numtmin2%), phi3#(numtmin2%), D#(numtmin2%), wt1%(numtmin1%), wt2%(numtmin2%)
PRINT "Entering T_mins"
PRINT "~~~~~~~~~~~~~~~"
PRINT
PRINT "Enter the filename and path of the .in2 input file (or ENTER for none):"
INPUT "   "; inpfile2$
IF inpfile2$ <> "" THEN
 OPEN inpfile2$ + ".in2" FOR INPUT AS #2
 FOR i = 1 TO numtmin1%
  LINE INPUT #2, temp$
  tminprim#(i) = VAL(temp$)
  LINE INPUT #2, temp$
  wt1%(i) = VAL(temp$)
 NEXT i
 FOR i = 1 TO numtmin2%
  LINE INPUT #2, temp$
  tminsec#(i) = VAL(temp$)
  LINE INPUT #2, temp$
  wt2%(i) = VAL(temp$)
 NEXT i
 CLOSE #2
ELSE
 REM
 REM *** Get Data by Hand Input and Write to a File for Future use (.in2) ***
 REM
 wrotein2 = 1
 OPEN outfile$ + ".in2" FOR OUTPUT AS #1
entry2:
 CLS
 PRINT "Enter T_min(s) for Primary Eclipse"
 PRINT "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
 PRINT
 FOR i = 1 TO numtmin1%
  INPUT " T_min I:  ", tminprim#(i)
  INPUT " Rel. Wt:  ", wt1%(i)
 NEXT i
 INPUT "Are all the values correct (Y/N)"; ans$
 IF (ans$ = "N" OR ans$ = "n") THEN GOTO entry2
 FOR i = 1 TO numtmin1%
  PRINT #1, tminprim#(i)
  PRINT #1, wt1%(i)
 NEXT i
 PRINT "<Press any key to continue.>"
 DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
entry3:
 CLS
 PRINT "Enter T_min(s) for Secondary Eclipse"
 PRINT "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
 PRINT
 FOR i = 1 TO numtmin2%
  INPUT "T_min II:  ", tminsec#(i)
  INPUT " Rel. Wt:  ", wt2%(i)
 NEXT i
 INPUT "Are all the values correct (Y/N)"; ans$
 IF (ans$ = "N" OR ans$ = "n") THEN GOTO entry3
 FOR i = 1 TO numtmin1%
  PRINT #1, tminsec#(i)
  PRINT #1, wt2%(i)
 NEXT i
 PRINT "<Press any key to continue.>"
 DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
 CLOSE #1
END IF
REM
REM *****************************************************************
REM *** Part 1 -- Calculations of O-C's, Epochs/Cycle #s, and D's ***
REM *****************************************************************
REM
REM *** Primary ***
REM
FOR i = 1 TO numtmin1%
 e1#(i) = FIX((tminprim#(i) - tmin1#) / p1#)
 phi1#(i) = ((tminprim#(i) - tmin1#) / p1#) - e1#(i)
 IF phi1#(i) > .9 THEN
  phi1#(i) = phi1#(i) - 1
  e1#(i) = e1#(i) + 1
 END IF
 IF phi1#(i) < -.9 THEN
  phi1#(i) = phi1#(i) + 1
  e1#(i) = e1#(i) - 1
 END IF
 oc1#(i) = phi1#(i) * p1#
NEXT i
REM
REM *** Secondary ***
REM
FOR i = 1 TO numtmin2%
 e2#(i) = FIX((tminsec#(i) - tmin2#) / p2#)
 phi2#(i) = ((tminsec#(i) - tmin2#) / p2#) - e2#(i)
 IF phi2#(i) > .9 THEN
  phi2#(i) = phi2#(i) - 1
  e2#(i) = e2#(i) + 1
 END IF
 IF phi2#(i) < -.9 THEN
  phi2#(i) = phi2#(i) + 1
  e2#(i) = e2#(i) - 1
 END IF
 oc2#(i) = phi2#(i) * p2#
NEXT i
REM
REM *** D for Secondary ***
REM
FOR i = 1 TO numtmin2%
 e3#(i) = FIX((tminsec#(i) - tmin1#) / p1#)
 phi3#(i) = ((tminsec#(i) - tmin1#) / p1#) - e3#(i)
 IF phi3#(i) < 0 THEN
  phi3#(i) = phi3#(i) + 1
  e3#(i) = e3#(i) - 1
 END IF
 D#(i) = (phi3#(i) - .5) * p1#
NEXT i
REM
REM *** Output from Part 1 ***
REM
REM *** Print to Screen ***
REM
CLS
IF numtmin1% <> 0 THEN
PRINT starname$; "'s O-C Information"
PRINT
PRINT "For Primary Eclipse:"
PRINT
PRINT "JD              O-C         Epoch   Rel.Wt."
FOR i = 1 TO numtmin1%
 PRINT USING "#######.####   "; tminprim#(i);
 PRINT USING "+##.#####   "; oc1#(i);
 PRINT USING "+####"; e1#(i);
 PRINT USING "      ##"; wt1%(i)
NEXT i
PRINT
PRINT "<Press any key to continue.>"
DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
END IF
CLS
IF numtmin2% <> 0 THEN
PRINT starname$; "'s O-C Information"
PRINT
PRINT "For Secondary Eclipse:"
PRINT
PRINT "JD              O-C         D          Epoch     Rel.Wt."
FOR i = 1 TO numtmin2%
 PRINT USING "#######.####   "; tminsec#(i);
 PRINT USING "+##.#####   "; oc2#(i);
 PRINT USING "##.#####   "; D#(i);
 PRINT USING "+####"; e3#(i);
 PRINT USING "       ##"; wt2%(i)
NEXT i
PRINT
PRINT "<Press any key to continue.>"
DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
END IF
REM
REM *** Print to File ***
REM
OPEN outfile$ + ".o-c" FOR OUTPUT AS #1
PRINT #1, starname$; "'s O-C Information"
PRINT #1, ""
PRINT #1, "        T_min (I) = "; USING "#######.#####"; tmin1#
PRINT #1, "            P (I) = "; USING "##.#########"; p1#
PRINT #1, "   # of T_min I's = "; USING "##"; numtmin1%
PRINT #1, "       T_min (II) = "; USING "#######.#####"; tmin2#
PRINT #1, "           P (II) = "; USING "##.#########"; p2#
PRINT #1, "  # of T_min II's = "; USING "##"; numtmin2%
PRINT #1, ""
PRINT #1, "Ephemeris calculated by:  "; ephem$
PRINT #1, ""
PRINT #1, "For Primary Eclipse:"
PRINT #1, ""
PRINT #1, "Julian Date         O-C (day)       Epoch    Rel.Wt."
FOR i = 1 TO numtmin1%
 PRINT #1, USING "#######.########   "; tminprim#(i);
 PRINT #1, USING "+##.#########   "; oc1#(i);
 PRINT #1, USING "+####"; e1#(i);
 PRINT #1, USING "      ##"; wt1%(i)
NEXT i
PRINT #1, ""
PRINT #1, "For Secondary Eclipse:"
PRINT #1, ""
PRINT #1, "Julian Date         O-C (day)       D (day)        Epoch     Rel.Wt."
FOR i = 1 TO numtmin2%
 PRINT #1, USING "#######.########   "; tminsec#(i);
 PRINT #1, USING "+##.#########   "; oc2#(i);
 PRINT #1, USING "##.#########   "; D#(i);
 PRINT #1, USING "+####"; e3#(i);
 PRINT #1, USING "       ##"; wt2%(i)
NEXT i
CLOSE #1
REM
REM *********************
REM *** End of Part 1 ***
REM *********************
REM
REM ****************************************
REM *** Part 2 -- Ephemeris Calculations ***
REM ****************************************
REM
REM *** Find the new ephemeris (do the liner least squares fit) ***
REM *** epoch = x, tmin = y; period = slope, newtmin = y intercept ***
REM
REM *** Primary ***
REM
sumxy1# = sumx1# = sumy1# = sumx21# = sumy21# = 0!
sumwt1% = 0
FOR i = 1 TO numtmin1%
  sumwt1% = sumwt1% + wt1%(i)
  sumxy1# = sumxy1# + (e1#(i) * tminprim#(i)) * wt1%(i)
  sumx1# = sumx1# + e1#(i) * wt1%(i)
  sumy1# = sumy1# + tminprim#(i) * wt1%(i)
  sumx21# = sumx21# + (e1#(i) * e1#(i)) * wt1%(i)
  sumy21# = sumy21# + (tminprim#(i) * tminprim#(i)) * wt1%(i)
NEXT i
delta1# = (sumwt1% * sumx21#) - (sumx1# * sumx1#)
newtmin1# = (sumx21# * sumy1# - sumx1# * sumxy1#) / delta1#
period1# = (sumwt1% * sumxy1# - sumx1# * sumy1#) / delta1#
sum2line1# = 0!
FOR i = 1 TO numtmin1%
 term1# = (tminprim#(i) - (newtmin1# + period1# * e1#(i))) * wt1%(i)
 sum2line1# = sum2line1# + (term1# * term1#)
NEXT i
sigmay21# = (1 / (sumwt1% - 2)) * sum2line1#
sigmatmin1# = SQR(sigmay21# * sumx21# / delta1#)
sigmaper1# = SQR(sumwt1% * sigmay21# / delta1#)
REM
REM *** Secondary ***
REM
sumxy2# = sumx2# = sumy2# = sumx22# = sumy22# = 0!
sumwt2% = 0
FOR i = 1 TO numtmin2%
  sumwt2% = sumwt2% + wt2%(i)
  sumxy2# = sumxy2# + (e2#(i) * tminsec#(i)) * wt2%(i)
  sumx2# = sumx2# + e2#(i) * wt2%(i)
  sumy2# = sumy2# + tminsec#(i) * wt2%(i)
  sumx22# = sumx22# + (e2#(i) * e2#(i)) * wt2%(i)
  sumy22# = sumy22# + (tminsec#(i) * tminsec#(i)) * wt2%(i)
NEXT i
delta2# = (sumwt2% * sumx22#) - (sumx2# * sumx2#)
newtmin2# = (sumx22# * sumy2# - sumx2# * sumxy2#) / delta2#
period2# = (sumwt2% * sumxy2# - sumx2# * sumy2#) / delta2#
sum2line2# = 0!
FOR i = 1 TO numtmin2%
 term2# = (tminsec#(i) - (newtmin2# + period2# * e2#(i))) * wt2%(i)
 sum2line2# = sum2line2# + (term2# * term2#)
NEXT i
sigmay22# = (1 / (sumwt2% - 2)) * sum2line2#
sigmatmin2# = SQR(sigmay22# * sumx22# / delta2#)
sigmaper2# = SQR(sumwt2% * sigmay22# / delta2#)
REM
REM *** Output from Part 2 ***
REM
REM *** Print to Screen ***
REM
CLS
PRINT
PRINT starname$; "'s Ephemeris Information"
PRINT
PRINT "For primary eclipse:"
PRINT "   T_min (I) = "; USING "#######.#######"; newtmin1#
PRINT " Uncertainty = "; USING "#######.#######"; sigmatmin1#
PRINT "  Period (I) = "; USING "##.#########"; period1#
PRINT " Uncertainty = "; USING "##.#########"; sigmaper1#
PRINT
PRINT "For secondary eclipse:"
PRINT "  T_min (II) = "; USING "#######.#######"; newtmin2#
PRINT " Uncertainty = "; USING "#######.#######"; sigmatmin2#
PRINT " Period (II) = "; USING "##.#########"; period2#
PRINT " Uncertainty = "; USING "##.#########"; sigmaper2#
PRINT
PRINT "<Press any key to continue.>"
DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
REM END IF
REM
REM *** Print to File ***
REM
OPEN outfile$ + ".eph" FOR OUTPUT AS #1
PRINT #1, starname$; "'s Ephemeris Information"
PRINT #1, ""
PRINT #1, "For primary eclipse:"
PRINT #1, "   T_min (I) = "; USING "#######.#######"; newtmin1#
PRINT #1, " Uncertainty = "; USING "#######.#######"; sigmatmin1#
PRINT #1, "  Period (I) = "; USING "##.#########"; period1#
PRINT #1, " Uncertainty = "; USING "##.#########"; sigmaper1#
PRINT #1, ""
PRINT #1, "For secondary eclipse:"
PRINT #1, "  T_min (II) = "; USING "#######.#######"; newtmin2#
PRINT #1, " Uncertainty = "; USING "#######.#######"; sigmatmin2#
PRINT #1, " Period (II) = "; USING "##.#########"; period2#
PRINT #1, " Uncertainty = "; USING "##.#########"; sigmaper2#
CLOSE #1
REM
REM *********************
REM *** End of Part 2 ***
REM *********************
REM
REM *****************************************************************
REM *** Part 3 -- Calculation of the Observed Apsidal Motion Rate ***
REM *****************************************************************
REM
wrad# = wdeg# * DEGTORAD
pbar# = (period1# + period2#) / 2
term1# = pbar# / (period1# - period2#)
term2# = (1 - (e# ^ 2)) ^ (3 / 2)
term3# = 4 * e# * SIN(wrad#)
term4# = (1 - ((e# ^ 2) * ((SIN(wrad#)) ^ 2))) ^ 2
Uday# = ((term2# * term3# * pbar#) / term4#) * term1#
Uyr# = Uday# / DAYSPERYEAR
omegadot# = (360 / Uyr#) * 100
REM
REM *** Uncertainty Calculation ***
REM
REM Step by step propagation for U = K1 * (p1 + p2)^2 * (p1 - p2)^(-1)
REM
psum# = period1# + period2#
pdif# = period1# - period2#
delpsumdif# = sigmaper1# + sigmaper2#: REM same error for sum and difference
fracpsum# = delpsumdif# / ABS(psum#)
fracpdif# = delpsumdif# / ABS(pdif#)
fracpsum2# = fracpsum# ^ 2
fracpdif2# = fracpdif# ^ 2
REM Next step uses power law, error in x^n = n * frac x
fracU2# = 2 ^ 2 * fracpsum2# + (-1) ^ 2 * fracpdif2#
fracU# = SQR(fracU2#)
REM frac U = sigmaU / abs(U), so sigmaU = U * fracU
sigmaUyr# = fracU# * ABS(Uyr#)
REM fracw is really the same as fracU, but since w = K2 / U, the powers are
REM different, so the equation is explicitly shown below (result is same)
fracw2# = (-2) ^ 2 * fracpsum2# + 1 ^ 2 * fracpdif2#
fracw# = SQR(fracw2#)
sigmaw# = fracw# * omegadot#
REM
REM *** Output from Part 3 ***
REM
REM *** Print to Screen ***
REM
CLS
PRINT
PRINT "Observed apsidal motion rate calculation for "; starname$
PRINT
PRINT "       Eccentricty of the orbit, e:  "; e#
PRINT "  Longitude of periastron, w (deg):  "; wdeg#
PRINT "   Period of primary eclipse (day):  "; USING "###.##########"; period1#
PRINT "            Uncertainty in primary:  "; USING "###.##########"; sigmaper1#
PRINT " Period of secondary eclipse (day):  "; USING "###.##########"; period2#
PRINT "          Uncertainty in secondary:  "; USING "###.##########"; sigmaper2#
PRINT
PRINT "           U_aps (years) = "; USING "######.#######"; Uyr#
PRINT "     Uncertainty (years) = "; USING "######.#######"; sigmaUyr#
PRINT "   omega dot (deg/100yr) = "; USING "###.##########"; omegadot#
PRINT " Uncertainty (deg/100yr) = "; USING "###.##########"; sigmaw#
PRINT
PRINT "<Press any key to continue.>"
DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
OPEN outfile$ + ".cal" FOR OUTPUT AS #1
PRINT #1, "Observed apsidal motion rate calculation for "; starname$
PRINT #1, ""
PRINT #1, "       Eccentricty of the orbit, e:  "; e#
PRINT #1, "  Longitude of periastron, w (deg):  "; wdeg#
PRINT #1, "   Period of primary eclipse (day):  "; USING "###.##########"; period1#
PRINT #1, "            Uncertainty in primary:  "; USING "###.##########"; sigmaper1#
PRINT #1, " Period of secondary eclipse (day):  "; USING "###.##########"; period2#
PRINT #1, "          Uncertainty in secondary:  "; USING "###.##########"; sigmaper2#
PRINT #1, ""
PRINT #1, "           U_aps (years) = "; USING "######.#######"; Uyr#
PRINT #1, "     Uncertainty (years) = "; USING "######.#######"; sigmaUyr#
PRINT #1, "   omega dot (deg/100yr) = "; USING "###.##########"; omegadot#
PRINT #1, " Uncertainty (deg/100yr) = "; USING "###.##########"; sigmaw#
CLOSE #1
REM
REM *********************
REM *** End of Part 3 ***
REM *********************
REM
REM *** Finish ***
CLS
PRINT
PRINT "Thank you for using this period study program / apsidal motion calculator."
PRINT "Here are the names of the all the files that I have saved for you:"
PRINT
IF wrotein1 = 1 THEN
 PRINT "  "; outfile$ + ".in1"; " -- all the intial inputs for this program"
END IF
IF wrotein2 = 2 THEN
 PRINT "  "; outfile$ + ".in2"; " -- all the eclipse timings (I and II)"
END IF
PRINT "  "; outfile$ + ".o-c"; " -- eclipse times, O-C's, epochs/cycle #s, D's"
PRINT "  "; outfile$ + ".eph"; " -- two new ephemeris for primary and secondary"
PRINT "  "; outfile$ + ".cal"; " -- the observed apsidal motion rate"
PRINT
PRINT "Have a nice day!  :)"
END