I want the solution for assignment in c++ and also I want graph in Matlab with coding

Ensuring the proper orbital insertion for a satellite is critical if the satellite is to function as intended. While space scientists learn a great deal about orbital mechanics, it is also important for aerospace engineers to have an understanding of the basics. In fact, you will learn these basics when you take MAE 4263Space Flight Mechanics.

As an introduction to orbital mechanics, you will be focusing on planetary orbits for this project. Elementary orbital mechanics begins with the study of the two-body problem, in which it is assumed that there are only two masses in the universe, with one orbiting about the other. In this case, orbital motion is described by the following ordinary differential equation:

r^→ = d²r/dt² = -(G(M + m)/r²).(r^→/r) = -(G(M + m)/r³)r^→

where r^→ is the vector between the two bodies, r^→ is the acceleration vector, M and m are the masses of the two bodies, and G is the universal gravitational constant. This equation may be expanded into four equations so as to be able to calculate the velocity of the orbiting body:

dr_x/dt = v_x, dr_y/dt = v_y

dv_x/dt = -G(M + m)r_x/r³, dv_y/dt = -G(M + m)r_y/r³

where r_x, r_y, v_x are the x- and y- components of the position and velocity vectors, respectively.

You are to solve these equations using the 4th-order Runge-Kutta scheme. You must write a computer program that accepts as input the aphelion (farthest) distance the planet is from the sun in astronomical units (AU), the eccentricity of the orbit, , and how many orbits, 0 ≤ e ≤ 1 are to be calculated. With the help of the equations on the back, your code will solve the four ODEs above, plotting the position of the orbiting body (starting at the aphelion) as a function of time as it orbits counterclockwise about the sun.

Your code should continue until a user-specified number of orbital periods (i.e. complete orbits) have been calculated. The time increment you use should be equal to 0.01 years. Also, when using the units of AU and years, G(M + m) = μ = 4Π² irrespective of the masses of the planetary bodies involved. Note that the position of the aphelion should be along the negative x-axis, and the initial velocity should only be in the negative y direction.

Case 1: r_a = 1.01671 AU, e = 0.0167 (Earth)
Case 2: r_a = 1.666 AU, e = 0.0934 (Mars)
Case 3: r_a = 35.100 AU, e = 0.967 (Halley's Comet, Dec. 2023)
Case 4: r_a = 27.530 AU, e = 0.744 (Some random asteroid I made up)

where r_a is the radius of aphelion. For each of these cases, compute and plot a single orbit.

Plot each orbit on its own graph. Be sure to scale your axes equally to see the true shape of the orbits. Also put a small symbol at (x,y) = (0,0) to denote the sun. Make sure axes have labels, the plot has a title, and that you use all other good plotting techniques.

Up to 10% Extra Credit: We have all heard the expression that "the gravitational force between two bodies is inversely-proportional to the distance between them squared." This is embodied in the governing ODE above, and dictates the motion of the planets and other orbiting bodies. But what if the laws of physics were different? What if the force weren't inversely-proportional to the distance squared, but instead to the distance raised to the 2.1 power? How would planets orbit the sun then? As extra credit, modify your code to explore these fantasy physics. Using the same data above for the asteroid's orbit (case 4), plot 6 periods of the new orbit for the asteroid if the force were inversely-proportional to the distance raised to the 2.1 power. Considering the number of data points this will generate, you may want to increase the time step to 0.1 years from 0.01 years for this problem only.

You may choose any computer language you prefer for your code provided you are not using a commercial or other pre-written application to do the work for you (MATLAB, Mathematica, Excel, and other available software are not acceptable). The plots may be created by importing your data into any available software. While working in small groups to set up the problem and determine how to solve it is acceptable, all code must be original work. Code that is copied or substantially similar to other students' work (past or present) will be rejected and no credit for the assignment received.

To complete this project, you must turn in both a hard copy of the code listing, the plots, and numeric values for rx, ry, vx, and vy for Case 1 (only) out to exactly six decimal places for each iteration using all proper formatting techniques for tabulated data. Everything must also be submitted electronically on Canvas. Please do not email your assignment to me. This project must be turned in by the end of class on the due date to receive credit. No late assignments will be accepted.