Two-Body-Simulation/two_body_simulation.py at main · cagintunc/Two-Body-Simulation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
import math
import TwoBodyView

gravity_force = 6.673 * (10**(-11))


class TwoBodyModel:

    def __init__(self, mass_1, mass_2):
        self.mass_1 = mass_1
        self.mass_2 = mass_2
        self.m12 = mass_1 + mass_2
        self.u = [1, 0, 0, 0]
        self.dimensions = [[0, 0], [0, 0]]


class TwoBodyController:

    def __init__(self, mass_1, mass_2, eccentricity = 0.7):
        self.new_model = TwoBodyModel(mass_1, mass_2)
        self.new_model.u[3] = eccentricity
        self.q = mass_1/mass_2

    def derivative(self):
        du = [0, 0, 0, 0]
        cordinate = self.new_model.u[0], self.new_model.u[1]
        distance = math.sqrt(math.pow(cordinate[0], 2) + math.pow(cordinate[1], 2))

        for i in range(2):
            du[i] = self.new_model.u[i + 2]
            du[i + 2] = -(1 + self.q) * cordinate[i] / (math.pow(distance, 3))
        return du

    def update_position(self, time_step,flag):
        if flag == 0:
            # runge-kutta
            self.runge_kutta(time_step)
        else:
            # euler
            self.euler(time_step)
        self.calculate_new_position()

    def runge_kutta(self, h):
        a = [h/2, h/2, h, 0]
        b = [h/6, h/3, h/3, h/6]
        u_0 = []
        u_t = []
        for i in range(len(self.new_model.u)):
            u_0.append(self.new_model.u[i])
            u_t.append(0)

        for i in range(4):
            du = self.derivative()
            for j in range(len(self.new_model.u)):
                self.new_model.u[j] = u_0[j] + (a[i]*du[j])
                u_t[j] = u_t[j] + (b[i]*du[j])
        for i in range(len(self.new_model.u)):
            self.new_model.u[i] = u_0[i] + u_t[i]

    def euler(self, h):
        du = self.derivative()
        for i in range(4):
            self.new_model.u[i] = self.new_model.u[i] + (h*du[i])

    def write_to_file(self):
        string = "f% f% f%\n"%(self.u[0], self.u[1], self.u[2], self.u[3])
        with open("data.txt", "a", encoding = "utf-8") as file:
            file.write(string)

    def calculate_new_position(self):
        r = 1
        a1 = (self.new_model.mass_2 / self.new_model.m12) * r
        a2 = (self.new_model.mass_1 / self.new_model.m12) * r

        self.new_model.dimensions[0][0] = -a2 * self.new_model.u[0]
        self.new_model.dimensions[0][1] = -a2 * self.new_model.u[1]
        self.new_model.dimensions[1][0] = a1 * self.new_model.u[0]
        self.new_model.dimensions[1][1] = a1 * self.new_model.u[1]


def get_data(file_name):
    data_list = list()
    with open(file_name, "r", encoding = "utf-8") as file:
        reader = file.readlines()
        index = 0
        for i in reader:
            reader[index] = i.rstrip().split(" ")
            data_list.append(reader[index])
    return data_list


def get_input_from_console():
    mass_ratio = 0.0
    first_body_cord = 250
    eccentricity = 0
    iteration = 0
    time_step = 0
    method = 0 # runge kutt
    while True:
        mass_ratio = float(input("Mass ratio:\n"))
        if 0 < mass_ratio <= 1:
            break
        else:
            print("Give a number between (0, 1]")
    while True:
        iteration = int(input("Please give a number of iteration:\n"))
        if 100 < iteration < 1000000:
            break
        else:
            print("Give a number between (0, 1000000)")
    while True:
        time_step = float(input("Please give a time step:\n"))
        if 0 < time_step < 1:
            break
        else:
            print("Give a number between (0, 1)")
    while True:
        eccentricity = int(input("Please give an eccentricity between mass (50-150):\n"))
        if 50 <= eccentricity <= 150:
            break
    while True:
        method = int(input("0 -> Runge kutta\n1 -> Euler\nMethod:\n"))
        if method == 0 or method == 1:
            break
    return first_body_cord, eccentricity, mass_ratio, iteration, time_step, method


inputs = get_input_from_console()
with open("data.txt", "w", encoding="utf-8") as file:
    control = TwoBodyController(1, inputs[2])
    control.positions = [inputs[0], 250, inputs[0] + inputs[1], 200]
    step = 0
    while step < inputs[3]:
        control.update_position(inputs[4], inputs[5])
        current_data = control.new_model.dimensions
        string = str(current_data[0][0]) + " " + str(current_data[0][1]) + " " + str(current_data[1][0]) + " " + str(current_data[1][1]) + "\n"
        file.write(string)
        step += 1
model = TwoBodyModel(control.new_model.mass_1, control.new_model.mass_2)
two_body = TwoBodyView.TwoBodyView(inputs[0], 350, inputs[0] + inputs[1], 300, int(15 * inputs[2]), 15)
data = get_data("data.txt")
two_body.circulation(data)