Application Center - Maplesoft
  Submit your work
Application Center Applications Binary Star or Star with Exoplanet: Fitting a Curve to the Orbital Data Preview

App Preview:

Binary Star or Star with Exoplanet: Fitting a Curve to the Orbital Data

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



Binary Star or Star with Exoplanet: Fitting a Curve to the Orbital Data ***

restart

Two problems are contained in this worksheet. Problem A involves fitting a curve to the redshift variation of the star 51 Pegasi, indicating that a body of significant mass is orbiting it. Problem B involves calculating the mass of the orbiting body to determine whether it is a star or a planet.

 

Problem A: The star 51 Pegasi shows a regular variation in redshift, indicating that it is orbiting around a common barycentre. This indicates the possibility that it is a single-line spectroscopic binary or that it has a planet of significant mass orbiting it. Fit a curve to the given data as a first step in determining which of these possibilities is the case.

 

Hints:

 

Plot the data and fit a curve to it by adjusting the parameters of the fitting function.

Plot the final fitting function.

Find the roots of the fitting function and set them to zero.

Take the difference of alternate roots representing one complete cycle.

 

 

Data:

The data (from Marcy et al., 1997) are given as follows with time given in julian days and velocities in m/s:NULL

 

Velocities := [-20.2, -8.1, 5.6, 56.4, 66.8, -35.1, -42.6, -33.5, -27.5, -22.7, 45.3, 47.6, 56.2, 65.3, 62.5, -22.6, 31.7, -44.1, -37.1, -35.3, 25.1, 35.7, 41.2, 61.3, 56.9, 51, -2.5, -4.6, -38.5, -48.7, 2.7, 17.6]

``

Days := [.6, .7, .8, 1.6, 1.7, 3.6, 3.7, 4.6, 4.7, 4.8, 5.6, 5.7, 5.8, 6.6, 6.7, 7.7, 7.8, 8.6, 8.7, 8.8, 9.6, 9.7, 9.8, 10.6, 10.7, 10.8, 11.7, 11.8, 12.6, 12.7, 13.6, 13.7]

``

Useful Equations:

NULL

FittedLine := fit[leastsquare[[z, x], x = z*RecipFreq+phi]]([ArcsinScaledDays, Velocities])

NULLNULL

NonlinearResiduals := [seq(A0*sin(f0*(Velocities[i]-phi0))+K0-Days[i], i = 1 .. nops(Velocities))]

``

lssoln := LSSolve(NonlinearResiduals, initialpoint = {A0 = A, K0 = K, f0 = f, phi0 = phi})

NULL

NULL

Solution:

 

Combine data into a list of lists and plot.

NULL

Pairs := zip(proc (x, y) options operator, arrow; [x, y] end proc, Days, Velocities)

[[.6, -20.2], [.7, -8.1], [.8, 5.6], [1.6, 56.4], [1.7, 66.8], [3.6, -35.1], [3.7, -42.6], [4.6, -33.5], [4.7, -27.5], [4.8, -22.7], [5.6, 45.3], [5.7, 47.6], [5.8, 56.2], [6.6, 65.3], [6.7, 62.5], [7.7, -22.6], [7.8, 31.7], [8.6, -44.1], [8.7, -37.1], [8.8, -35.3], [9.6, 25.1], [9.7, 35.7], [9.8, 41.2], [10.6, 61.3], [10.7, 56.9], [10.8, 51], [11.7, -2.5], [11.8, -4.6], [12.6, -38.5], [12.7, -48.7], [13.6, 2.7], [13.7, 17.6]]

(1)

 

plot(Pairs, symbol = DIAMOND, style = POINT)

 

Note that the first 14 data points represent one period. Therefore, make lists of the first 14 data points:

 

NewVelocities := [-20.2, -8.1, 5.6, 56.4, 66.8, -35.1, -42.6, -33.5, -27.5, -22.7, 45.3, 47.6, 56.2, 65.3]

[-20.2, -8.1, 5.6, 56.4, 66.8, -35.1, -42.6, -33.5, -27.5, -22.7, 45.3, 47.6, 56.2, 65.3]

(2)

NULL

NewDays := [.6, .7, .8, 1.6, 1.7, 3.6, 3.7, 4.6, 4.7, 4.8, 5.6, 5.7, 5.8, 6.6]

[.6, .7, .8, 1.6, 1.7, 3.6, 3.7, 4.6, 4.7, 4.8, 5.6, 5.7, 5.8, 6.6]

(3)

NULL

Combine these into pairs and plot:

NULL

  NewPairs := zip(proc (x, y) options operator, arrow; [x, y] end proc, NewDays, NewVelocities)

[[.6, -20.2], [.7, -8.1], [.8, 5.6], [1.6, 56.4], [1.7, 66.8], [3.6, -35.1], [3.7, -42.6], [4.6, -33.5], [4.7, -27.5], [4.8, -22.7], [5.6, 45.3], [5.7, 47.6], [5.8, 56.2], [6.6, 65.3]]

(4)

NULL

plot(NewPairs, symbol = DIAMOND, style = POINT)

NULL

Rescale the velocities about the mean, convert values to arcsin, fit to a function, and plot.

 

with(stats)

NULL

K := describe[mean](NewVelocities)

10.96428571

(5)

NULL

RescaledVelocities := map(proc (y) options operator, arrow; y-K end proc, NewVelocities)

[-31.16428571, -19.06428571, -5.36428571, 45.43571429, 55.83571429, -46.06428571, -53.56428571, -44.46428571, -38.46428571, -33.66428571, 34.33571429, 36.63571429, 45.23571429, 54.33571429]

(6)

NULL

ScaledVelocityRange := describe[range](RescaledVelocities)

-53.56428571 .. 55.83571429

(7)

NULL

A := max(abs(op(1, ScaledVelocityRange)), abs(op(2, ScaledVelocityRange)))

55.83571429

(8)

ArcsinScaledVelocities := map(proc (y) options operator, arrow; arcsin((y-K)/A) end proc, NewVelocities)

[-.5921454826, -.3484435788, -0.9622107038e-1, .9505568454, 1.570796327, -.9701965406, -1.284581659, -.9212219492, -.7599472418, -.6471520151, .6623121520, .7156844237, .9444204824, 1.338478962]

(9)

NULL

FittedLine := fit[leastsquare[[z, x], x = z*RecipFreq+phi]]([ArcsinScaledVelocities, NewDays])

x = .3718543916*z+3.592206538

(10)

NULL

f := 1/coeff(rhs(FittedLine), z, 1)

2.689224660

(11)

NULLNULL

phi := coeff(rhs(FittedLine), z, 0)

3.592206538

(12)

NULL

DataPlot := plot(NewPairs, symbol = DIAMOND, style = POINT)

``

CurvePlot := plot(A*sin(f*(x-phi))+K, x = NewDays[1] .. NewDays[nops(NewDays)], color = blue)

``

plots[display](DataPlot, CurvePlot)

NULL

NULL

Clearly, the curve does not fit the data. Experiment with different values of f until you find a function that is close to the data points. Try a value of 1.5.

 

f := 1.5

1.5

(13)

NULL

DataPlot := plot(NewPairs, symbol = DIAMOND, style = POINT)

NULL

CurvePlot := plot(A*sin(f*(x-phi))+K, x = NewDays[1] .. NewDays[nops(NewDays)], color = blue)

NULL

plots[display](DataPlot, CurvePlot)

 

NULL

The curve seems to be close to the desired form, except that it is displaced from the points. Try altering the value of phi. Try a value of 9.35.

NULL

phi := 9.35

9.35

(14)

 

DataPlot := plot(NewPairs, symbol = DIAMOND, style = POINT)

``NULL

CurvePlot := plot(A*sin(f*(x-phi))+K, x = NewDays[1] .. NewDays[nops(NewDays)], color = blue)

``

plots[display](DataPlot, CurvePlot)

 

NULL

This is much closer to a true fit. Try optimization.

``

with(Optimization)

``

NonlinearResiduals := [seq(A0*sin(f0*(NewDays[i]-phi0))+K0-NewVelocities[i], i = 1 .. nops(NewDays))]

``

phi := 9.35

9.35

(15)

f := 1.5

1.5

(16)

lssoln := LSSolve(NonlinearResiduals, initialpoint = {A0 = A, K0 = K, f0 = f, phi0 = phi})

[52.4185740754240, [A0 = HFloat(60.03138206441081), K0 = HFloat(9.323116294488084), f0 = HFloat(1.4653063322023352), phi0 = HFloat(9.469852810790444)]]

(17)

NULL

CurvePlot := plot(subs(lssoln[2], A0*sin(f0*(x-phi0))+K0), x = NewDays[1] .. NewDays[nops(NewDays)], color = blue)

NULL

oldplot := plot(NewPairs, symbol = DIAMOND, style = POINT)

plots[display](oldplot, CurvePlot)

NULL

NULL

This is a close fit. Now, show a plot containing all the original data points together with the plot of the fitted  curve.

``

CurvePlot := plot(subs(lssoln[2], A0*sin(f0*(x-phi0))+K0), x = Days[1] .. Days[nops(Days)], color = blue)

``

oldplot := plot(Pairs, symbol = DIAMOND, style = POINT)

plots[display](oldplot, CurvePlot)

``

Show the final fitting function and plot it from -2π to +2π (the default).

``

subs(lssoln[2], A0*sin(f0*(x-phi0))+K0)

HFloat(60.03138206441081)*sin(HFloat(1.4653063322023352)*x-HFloat(13.87623528867532))+HFloat(9.323116294488084)

(18)

"->"

 

Differentiate the function and find the zero roots.

 

diff(60.0313820644108*sin(1.46530633220234*x-13.8762352886753)+9.32311629448808, x)

87.96436425*cos(1.46530633220234*x-13.8762352886753)

(19)

NULL

with(Student[Calculus1])

NULL

Roots(87.96436425*cos(1.46530633220234*x-13.8762352886753))

[-8.754007258, -6.610023720, -4.466040183, -2.322056645, -.1780731078, 1.965910430, 4.109893967, 6.253877505, 8.397861042]

(20)

NULL

Take the difference of alternate roots representing one complete cycle.

 

8.397861042-4.109893967

4.287967075

(21)

NULL

The orbital period of 51 Pegasi appears to be 4.29 days. The Exoplanet catalogue gives 4.23 days (http://exoplanet.eu/catalog/51_peg_b/).

 

The velocity of 51 Pegasi around its barycentre is 57.75 m/s (Marcy et al.,1997). The shape of the curve indicates a nearly circular orbit of 51 Pegasi's companion. The star itself, through spectral analysis, is known to be a G2 IV main sequence star and has, therefore, approximately 1.11 solar masses. The accepted value is1.06 solar masses.

 

 

 

Problem 2: Use the data given below to determine the mass and the orbital radius of the object orbiting 51 Pegasi.

 

 

Hints:

 

An orbiting object less massive than Jupiter must be a planet.

Make the following assumptions: 1. that the object's orbit is circular, as indicated by the fitting function, above; 2. that the orbit is being viewed edge-on; 3. that the object's mass is an insignificant fraction of the mass of the star.

 

Data:

 

M[s] := (1.06*1.989)*10^30*Unit('kg')

0.2108340000e31*Units:-Unit('kg')

(22)

P := (60*(4.23*24)*60)*Unit('s')

365472.00*Units:-Unit('s')

(23)

G := 6.67*10^(-11)*Unit('N')*Unit('m')^2/Unit('kg')^2

0.6670000000e-10*Units:-Unit('N')*Units:-Unit('m')^2/Units:-Unit('kg')^2

(24)

v[s] := 54.87*Unit('m')/Unit('s')

54.87*Units:-Unit('m')/Units:-Unit('s')

(25)

NULL

Useful Equations:

 

NULL

NULL

NULL

NULL

NULL

NULL

NULL

NULL

Solution:

 

Find the radius of 51 Pegasi's orbit around its barycentre.

NULL

NULL

r[s] := P*v[s]/(2*Pi)

10026724.32*Units:-Unit('m')/Pi

(26)

Find the radius of the orbit of the orbiting object.

 

r[p] := (P^2/(4*Pi^2/(G*M[s])))^(1/3)

6645559888.*4^(2/3)*(Units:-Unit('s')^2*Units:-Unit('N')*Units:-Unit('m')^2/(Pi^2*Units:-Unit('kg')))^(1/3)

(27)

"(->)"

7806774868.*Units:-Unit('m')

(28)

NULL

Convert this value to astronomical units.

 

7.806774868*10^9*Unit('m')/(1.496*10^11*Unit('m'))

0.5218432398e-1

(29)

``

NULL

Find the mass of the orbiting object.

``NULL

``

M[p] := M[s]*r[s]/r[p]

0.7952581990e27*Units:-Unit('kg')*Units:-Unit('m')*4^(1/3)/(Pi*(Units:-Unit('s')^2*Units:-Unit('N')*Units:-Unit('m')^2/(Pi^2*Units:-Unit('kg')))^(1/3))

(30)

"(->)"

0.8619423008e27*Units:-Unit('kg')

(31)

NULL

Convert the mass of the orbiting object to units of Jovian mass.

NULL

8.619423008*10^26*Unit('kg')/(1.898*10^27*Unit('kg'))

.4541318761

(32)

NULLNULL

If these assumptions given earlier do not hold, the equations used above must be modified. For example, the given mass of 0.45 Jupiter masses is a lower limit on the mass. If the orbital inclination is greater than zero, a higher mass will be calculated.

 

    The object is 0.052 AU from 51 Pegasi. Finding a large, Jupiter-type planet so close to its star was surprising. However, its status as a planet was confirmed by Marcy et al., (1997).  It was the first to be found (by Mayor and Queloz, 1995) orbiting a normal star.

 

--------------------------------------------------------------------

References

 

Marcy, G., Butler, R., Williams, E., Bildsten, L., Graham, J., Ghez, A., and Jernigan, J. (1997). The Planet around 51 Pegasi*. ApJ, 481, 2, 926.

 

Mayor, M. and Queloz, D. (1995). A Jupiter-mass companion to a solar-type star. Nature 378, 355-359