ECLIPSING BINARY STARS

BASIC Subroutine

Sub SimpleModel (MR, MB, RR, RB, LR, LB, INC, NP, LITE(), XR(), YR(), ZR(), XB(), YB(), ZB())

' ********************************************************
' This subroutine generates a light curve for an eclipsing
' binary star system assuming spherical stars and circular
' orbits.  - Dan Bruton, October 1995, astro@tamu.edu
' ********************************************************
' MR,MB = Star Masses in arbitrary units
' RR,RB = Star Radii in fractions of the semimajor axis length
' LR,LB = Star Luminosities in arbitrary units
' INC = Orbital Inclination
' NP = Number of points to compute for one orbit
' LITE = Array containing the light curve
' XR,YR,ZR,XB,YB,ZB = Star Positions
'
' Note: The second character,"R" or "B", means red or blue star
' which are the colors used for the animations of the stars in orbit.
' ********************************************************

' Require that R1>R2 to simplify calculations
' PHI = Starting angle for orbit
If RR > RB Then
  M1 = MR
  M2 = MB
  R1 = RR
  R2 = RB
  L1 = LR
  L2 = LB
  PHI = Pi
Else
  M1 = MB
  M2 = MR
  R1 = RB
  R2 = RR
  L1 = LB
  L2 = LR
  PHI = 0
End If

' Mass Ratio
Q = M2 / M1

' Semimajor Axis
SMA = 1

' Orbital Inclination
IR = INC * Pi / 180

' Compute Star Positions and Light Curve
For J = 1 To NP
  T = PHI + (J - 1) * 2 * Pi / NP: ' Phase Angle
  
  X = SMA * Cos(T)
  Y = SMA * Cos(IR) * Sin(T)
  Z = SMA * Sin(IR) * Sin(T)
  
  X1 = -1 * X / ((1 / Q) + 1)
  Y1 = -1 * Y / ((1 / Q) + 1)
  Z1 = -1 * Z / ((1 / Q) + 1)
     
  X2 = X / (1 + Q)
  Y2 = Y / (1 + Q)
  Z2 = Z / (1 + Q)

' Apparent Distance Between Stars
  RHO = Sqr((X2 - X1) ^ 2 + (Y2 - Y1) ^ 2)
    
'
' Area of Star Disks (not during an eclipse)
'
  A1 = Pi * (R1 ^ 2)
  A2 = Pi * (R2 ^ 2)
    
'
' Total and Annular Eclipses
'
If RHO < (R1 + R2) Then
    If RHO < (R1 - R2) Then
      If Z1 > Z2 Then
        A2 = 0            ' Total Eclipse
      Else
        A1 = (A1 - A2)    ' Annular Eclipse
      End If
    Else
'
' Partial Eclipses
'   The area covered during partial eclipses is determined
' by considering and area of a circle cut by a line segment:
' Area = r^2 (theta - Sin(theta)) / 2 where r is the radius
' of the circle and sin(theta)=h/r.  See math handbooks.
'
      H = ((R1 ^ 2) * (R2 ^ 2))
      H = H - (.25 * (((RHO ^ 2) - (R1 ^ 2) - (R2 ^ 2)) ^ 2))
      If H <> R2 And H <> R1 And H >= 0 Then
        H = Sqr(H / (RHO ^ 2))
        T1 = 2 * Atn(H / (Sqr((R1 ^ 2) - (H ^ 2))))
        T2 = 2 * Atn(H / (Sqr((R2 ^ 2) - (H ^ 2))))
      End If
      CA1 = (R1 ^ 2) * (T1 - Sin(T1)) / 2
      CA2 = (R2 ^ 2) * (T2 - Sin(T2)) / 2
' Shallow Partial Eclipse
      If (RHO > (Sqr((R1 ^ 2) - (R2 ^ 2)))) Then
          If Z1 > Z2 Then
            A2 = (A2 - CA1 - CA2)
          Else
            A1 = (A1 - CA1 - CA2)
          End If
' Deep Partial Eclipse
      Else
          If Z1 > Z2 Then
            A2 = (CA2 - CA1)
          Else
            A1 = (A1 - A2 + CA2 - CA1)
          End If
      End If
    End If
  End If
' Approximate Light Intensity (no limb darkening)
LITE(J) = ((L1 * A1 / R1 ^ 2) + (L2 * A2 / R2 ^ 2)) / (4 * Pi)

If RR > RB Then
  XR(J) = X1
  YR(J) = Y1
  ZR(J) = Z1
  XB(J) = X2
  YB(J) = Y2
  ZB(J) = Z2
Else
  XR(J) = X2
  YR(J) = Y2
  ZR(J) = Z2
  XB(J) = X1
  YB(J) = Y1
  ZB(J) = Z1
End If

Next J

End Sub