10 REM Program for simplified analysis of partial eclipsing binaries.
20 REM Based on 1992 AJ, 104(1), "A Simple Method for Preliminary Analysis of
30 REM the Light Curves of Partial Eclipsing Binaries" by Nematullah Riazi.
40 REM Programmed at Villanova University Astronomy & Astrophysics Department
50 REM By James Marshall (class of 1996), summer 1995.
60 REM Note an error in equation (6) of the paper, the term (lbp - lvs)
70 REM *should* be (lbp - lvp).  The correct equation is used here.
80 REM print instructions for use
90 CLS
100 PRINT "PRELIMINARY ANALYSIS OF THE LIGHT CURVE OF PARTIAL ECLIPSING BINARIES"
110 PRINT "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
120 PRINT
130 PRINT "This program uses a simple method for analyzing the light curves"
140 PRINT "of partial eclipsing binaries, as given by Riazi (1992 AJ, 104)."
150 PRINT "It assumes spherical components, neglects limb darkening and"
160 PRINT "proximity effects, but is still useful for getting a good estimate"
170 PRINT "of the values needed for more detailed analysis methods."
180 PRINT
190 PRINT "<Press any key for instructions>"
200 DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
210 CLS
220 PRINT "INSTRUCTIONS"
230 PRINT "~~~~~~~~~~~~"
240 PRINT
250 PRINT "You will need two light curves, taken in different filters, in"
260 PRINT "order to run this program.  You must input the system light at"
270 PRINT "primary and secondary minima for each of the two filters."
280 PRINT "These are denoted as l, with a p for primary minimum and an s for"
290 PRINT "secondary minimum.  The filters will be specified as filters 1 and 2"
300 PRINT "instead of by letter."
310 PRINT
320 PRINT "You must also enter a value for the phase angle at external"
330 PRINT "contact; this value is obtainable from the observed light curves"
340 PRINT "and will be denoted as phi_1."
350 PRINT
360 PRINT "You will also be given the opportunity to enter a value for the"
370 PRINT "orbital inclination angle (i), if you so desire.  This need not be"
380 PRINT "close to the final value, so if you choose not to enter a value,"
390 PRINT "a default value will be used."
400 PRINT
410 PRINT "<Press any key to input values>"
420 DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
430 REM get user input
440 CLS
450 PRINT "Please input the following values:"
460 PRINT
470 INPUT "     l_p, filter 1:  ", lp1!
480 INPUT "     l_s, filter 1:  ", ls1!
490 INPUT "     l_p, filter 2:  ", lp2!
500 INPUT "     l_s, filter 2:  ", ls2!
510 INPUT "       phi_1 (deg):  ", phi1!
520 PRINT "The following value is optional."
530 PRINT "Hit enter to use the default value of 80 degrees."
540 INPUT "           i (deg):  ", i!
550 IF (i!) THEN 520 ELSE i! = 80
560 PRINT
570 INPUT "Are all the values correct (Y/N)"; ans$
580 IF (ans$ = "N" OR ans$ = "n") THEN GOTO 430
590 PRINT
600 PRINT "Thank you.  Press any key for the results."
610 DO: anykey$ = INKEY$: LOOP WHILE anykey$ = ""
620 CLS
630 REM calculate L11, L12, L21, L22, k1, k2
640 REM L11 and L12
650 term1! = (lp2! - lp1!) * (1 - ls1!)
660 term2! = (ls1! - ls2!) * (1 - lp1!)
670 term3! = (term1! / term2!) + 1
680 L11! = 1 / term3!
690 term1! = (lp1! - lp2!) * (1 - ls2!)
700 term2! = (ls2! - ls1!) * (1 - lp2!)
710 term3! = (term1! / term2!) + 1
720 L12! = 1 / term3!
730 REM L21 and L22
740 L21! = 1 - L11!
750 L22! = 1 - L12!
760 REM k1 and k2
770 term1! = (1 - lp1!) / (1 - ls1!)
780 term2! = L21! / L11!
790 term3! = term1! * term2!
800 k1! = SQR(term3!)
810 term1! = (1 - lp2!) / (1 - ls2!)
820 term2! = L22! / L12!
830 term3! = term1! * term2!
840 k2! = SQR(term3!)
850 REM do iterative procedure to find r_1, r_2, and i
860 REM define constants and convert from degrees to radians
870 CONST PI = 3.14159265359#
880 CONST DEGTORAD = PI / 180
890 CONST RADTODEG = 180 / PI
900 phi1rad! = phi1! * DEGTORAD
910 irad! = i! * DEGTORAD
920 delta1! = COS(irad!): REM initial value same for both filters
930 delta2! = COS(irad!)
940 REM start of iteration on filter 1
950 REM calculate r1 and r2 for filter 1
960 term1! = (1 / (1 + k1!))
970 term2! = SQR(SIN(phi1rad!) ^ 2 + COS(phi1rad!) ^ 2 * delta1! ^ 2)
980 r11! = term1! * term2!
990 r21! = k1! * r11!
1000 REM calculate new delta value for filter 1
1010 S01! = PI * r11! ^ 2 * (1 - lp1!) / L11!
1020 arg! = (r11! ^ 2 - r21! ^ 2 + delta1! ^ 2) / (2 * r11! * delta1!)
1030 zeta! = ATN(-arg! / SQR(-arg! * arg! + 1)) + (PI / 2): REM arccos(arg)
1040 arg! = (r21! ^ 2 - r11! ^ 2 + delta1! ^ 2) / (2 * r21! * delta1!)
1050 xi! = ATN(-arg! / SQR(-arg! * arg! + 1)) + (PI / 2): REM arccos(arg)
1060 term1! = .5 * r11! ^ 2 * (2 * zeta! - SIN(2 * zeta!))
1070 term2! = .5 * r21! ^ 2 * (2 * xi! - SIN(2 * xi!))
1080 Sdelta1! = term1! + term2!
1090 term1! = -2 * r11! * r21! / delta1!
1100 term2! = SIN(zeta! + xi!)
1110 Sdeltaprime1! = term1! * term2!
1120 newdelta1! = delta1! - ((Sdelta1! - S01!) / Sdeltaprime1!)
1130 delta1! = newdelta1!
1140 IF (ABS(Sdelta1! - S01!) > .0000001) THEN GOTO 1040
1150 REM i1 = arccos(delta1) since delta1 = cos(i1)
1160 i1rad! = ATN(-delta1! / SQR(-delta1! * delta1! + 1)) + (PI / 2)
1170 i1! = i1rad! * RADTODEG
1180 REM r11 and r21 are same as last values calculated
1190 REM filter 1 is now done
1200 REM start of interation on filter 2
1210 REM calculate r1 and r2 for filter 2
1220 term1! = (1 / (1 + k2!))
1230 term2! = SQR(SIN(phi1rad!) ^ 2 + COS(phi1rad!) ^ 2 * delta2! ^ 2)
1240 r12! = term1! * term2!
1250 r22! = k2! * r12!
1260 REM calculate new delta value for filter 2
1270 S02! = PI * r12! ^ 2 * (1 - lp2!) / L12!
1280 arg! = (r12! ^ 2 - r22! ^ 2 + delta2! ^ 2) / (2 * r12! * delta2!)
1290 zeta! = ATN(-arg! / SQR(-arg! * arg! + 1)) + (PI / 2): REM arccos(arg)
1300 arg! = (r22! ^ 2 - r12! ^ 2 + delta2! ^ 2) / (2 * r22! * delta2!)
1310 xi! = ATN(-arg! / SQR(-arg! * arg! + 1)) + (PI / 2): REM arccos(arg)
1320 term1! = .5 * r12! ^ 2 * (2 * zeta! - SIN(2 * zeta!))
1330 term2! = .5 * r22! ^ 2 * (2 * xi! - SIN(2 * xi!))
1340 Sdelta2! = term1! + term2!
1350 term1! = -2 * r12! * r22! / delta2!
1360 term2! = SIN(zeta! + xi!)
1370 Sdeltaprime2! = term1! * term2!
1380 newdelta2! = delta2! - ((Sdelta2! - S02!) / Sdeltaprime2!)
1390 delta2! = newdelta2!
1400 IF (ABS(Sdelta2! - S02!) > .0000001) THEN GOTO 1300
1410 REM i2 = arccos(delta2) since delta2 = cos(i2)
1420 i2rad! = ATN(-delta2! / SQR(-delta2! * delta2! + 1)) + (PI / 2)
1430 i2! = i2rad! * RADTODEG
1440 REM r12 and r22 are same as last values calculated
1450 REM filter 2 is now done
1460 REM print out first set of results
1470 PRINT "Here are the results:"
1480 PRINT
1490 PRINT "  L1, filter 1 = "; USING "##.###"; L11!
1500 PRINT "  L1, filter 2 = "; USING "##.###"; L12!
1540 PRINT "  L2, filter 1 = "; USING "##.###"; L21!
1550 PRINT "  L2, filter 2 = "; USING "##.###"; L22!
1560 PRINT "   k, filter 1 = "; USING "##.###"; k1!
1570 PRINT "   k, filter 2 = "; USING "##.###"; k2!
1580 PRINT
1590 PRINT "  r1, filter 1 = "; USING "##.###"; r11!
1600 PRINT "  r1, filter 2 = "; USING "##.###"; r12!
1610 PRINT "  r2, filter 1 = "; USING "##.###"; r21!
1620 PRINT "  r2, filter 2 = "; USING "##.###"; r22!
1630 PRINT "   i, filter 1 = "; USING "##.###"; i1!
1640 PRINT "   i, filter 2 = "; USING "##.###"; i2!
1650 PRINT
1660 REM closing procedure
1670 PRINT "Thank you for using this preliminary analysis program."
1680 PRINT
1690 INPUT "Would you like to save these values to a file (Y/N)"; ans$
1700 IF (ans$ = "N" OR ans$ = "n") THEN GOTO 1960
1710 PRINT "Enter the filename with path (.ecl extension automatic):  "
1720 PRINT "   "; : INPUT file$
1730 OPEN file$ + ".ecl" FOR OUTPUT AS #1
1740 PRINT #1, "Output of Preliminary Analysis Program"
1750 PRINT #1, ""
1760 PRINT #1, "     l_1p:  "; lp1!
1770 PRINT #1, "     l_1s:  "; ls1!
1780 PRINT #1, "     l_2p:  "; lp2!
1790 PRINT #1, "     l_2s:  "; ls2!
1800 PRINT #1, "    phi_1:  "; phi1!
1810 PRINT #1, ""
1820 PRINT #1, "  L1, filter 1 = "; L11!
1830 PRINT #1, "  L1, filter 2 = "; L12!
1840 PRINT #1, "  L2, filter 1 = "; L21!
1850 PRINT #1, "  L2, filter 2 = "; L22!
1860 PRINT #1, "   k, filter 1 = "; k1!
1870 PRINT #1, "   k, filter 2 = "; k2!
1880 PRINT #1, ""
1890 PRINT #1, "  r1, filter 1 = "; r11!
1900 PRINT #1, "  r1, filter 2 = "; r12!
1910 PRINT #1, "  r2, filter 1 = "; r21!
1920 PRINT #1, "  r2, filter 2 = "; r22!
1930 PRINT #1, "   i, filter 1 = "; i1!
1940 PRINT #1, "   i, filter 2 = "; i2!
1950 CLOSE #1
1960 PRINT
1970 PRINT "Have a nice day!  :)"
1980 END