Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- from time import *
- from numba import njit as njit
- import numba.types
- from math import *
- import matplotlib.pyplot as plt
- test=False
- G=6.67e-11 #Постоянная
- #Луна
- M=7.3477e22
- R=1737e3
- @njit
- def integrate_kolenka(m1,h1,v1,dm,t1,U,N=100000,M=7.3477e22,R=1737e3):
- # Параметры Луны по умолчанию вобью
- #dm - масса топлив
- #t1 - время сгорания
- #U - скорость истечения газов
- h=h1
- v=v1
- for c in range(N):
- t=t1*c/N
- dt=t1/N
- ev=U*dm/(m1*t1-t*dm)-G*M/(R+h)**2
- eh=v
- v+=ev*dt
- h+=eh*dt
- #if (c==0):
- # print(f'{ev=} {dt=}')
- if h<0: #Если ракета не взлетает
- h=0
- if (v>1):#Если слишком жёсткое падение
- print("Ракета упала со скоростью: ", abs(v))
- v=max(v,0)
- v2=v
- h2=h
- return h2,v2,m1-dm
- if test:
- print("-----------------------------")
- #Проверим интегратор в невесомости с помощью формулы Цилковского про дельту
- m=10000
- v1=0
- h1=0
- h2,v2,m2=integrate_kolenka(m,h1,v1,dm=7000,t1=20,U=800,M=0)
- v_c2=800*log(10000/3000)
- print(f'{h2=:.1f} {v2=:.3f} {m2=:.1f} Цикловкский: v={v_c2:.3f}')
- h3,v3,m3=integrate_kolenka(m2,h2,v2,dm=2000,t1=20,U=1200,M=0)
- v_c3=v_c2+1200*log(3000/1000)
- print(f'{h3=:.1f} {v3=:.3f} {m3=:.1f} Цикловкский: v={v_c3:.3f}')
- # Выводит это:
- # h2=7744.0 v2=963.172 m2=3000.0 Цикловкский: v=963.178
- # h3=37824.0 v3=2281.498 m3=1000.0 Цикловкский: v=2281.513
- # Что меня вполне устраивает
- def max_h(h,v):
- e=v**2/2 #Энергия кинетическая
- e_h_max= G*M*(1/R)
- if e>e_h_max:
- # 2*e=v**2
- v=((e-e_h_max)*2)**0.5
- return f"Улетит на бесконечность со скорость {v:.1f}"
- else:
- if (v):
- h_max=(R+h)/((2*G*M/((R+h)*v**2))-1)
- return f"Улетит на высоту {h_max:.1f}"
- else:
- return f"Будет падать с высоты {h:.1f}"
- if test:
- print("-----------------------------")
- #Проверю формулу
- for v in [1,10,100,1000,2350,2400,6000]:
- #2380 - вторая космическая на луне
- g=M*G/R/R #Ускорение свободного падения на луне
- h_max=v*v/(2*g) #Высота при постоянно g
- print(f'{v=:.1f} {max_h(0,v)} по формуле: {h_max:.1f}')
- print(f'Остаток скорости при 6000: {(6000**2-2380**2)**0.5:.1f}')
- """ Вывод
- v=1.0 Улетит на высоту 0.3 по формуле: 0.3
- v=10.0 Улетит на высоту 30.8 по формуле: 30.8
- v=100.0 Улетит на высоту 3083.6 по формуле: 3078.2
- v=1000.0 Улетит на высоту 374114.3 по формуле: 307816.9
- v=2350.0 Улетит на высоту 79629015.7 по формуле: 1699918.6
- v=2400.0 Улетит на бесконечность со скорость 342.1 по формуле: 1773025.1
- v=6000.0 Улетит на бесконечность со скорость 5509.7 по формуле: 11081406.6
- Остаток скорости при 6000: 5507.8"""
- # Вроде верно, на низких высотах совпадает, при приближении ко второй космической
- # высота сильно отклоняется, а потом уходит на бесконечность
- if test:
- print("-----------------------------")
- #Проверим мою гипотезу, что время сгорания существенно для итогового результата
- m=10000
- v1=0
- h1=0
- for t in [1,10,100,980,1000,1020,10000]:
- h2,v2,m2=integrate_kolenka(m,h1,v1,dm=7000,t1=t,U=800)
- print(f'Сгорание за {t:5} секунд: {h2= :7.1f} {v2= :7.3f} {m2= :7.1f} {max_h(h2,v2)}')
- """ У меня выводит такое
- Сгорание за 1 секунд: h2= 386.4 v2= 961.548 m2= 3000.0 Улетит на высоту 340533.2
- Сгорание за 10 секунд: h2= 3790.9 v2= 946.950 m2= 3000.0 Улетит на высоту 329742.9
- Сгорание за 100 секунд: h2= 30637.5 v2= 802.396 m2= 3000.0 Улетит на высоту 232197.6
- Сгорание за 980 секунд: h2= 226.5 v2= 9.633 m2= 3000.0 Улетит на высоту 28.6
- Сгорание за 1000 секунд: h2= 153.8 v2= 7.391 m2= 3000.0 Улетит на высоту 16.8
- Сгорание за 1020 секунд: h2= 98.9 v2= 5.470 m2= 3000.0 Улетит на высоту 9.2
- Сгорание за 10000 секунд: h2= 0.0 v2= 0.000 m2= 3000.0 Будет падать с высоты 0.0"""
- #При сгорании за 1-10 секунд ничего не меняется, при 100 скорость меньше, но успевает набраться высота
- #При сгорании за 1000 секунд едва отрывается (это случайно получилось,
- # что при таких параметрах оно едва-едва отрывается)
- #При сгорании за 10000 секунд не взлетает
- if False:
- for t in np.linspace(0.1,0.9,20):
- m=1000
- c1=20
- s=40
- a1=a2=0.6
- U=3326
- m1=m*c1*t
- m2=m*c1*(1-t)
- mm=m+m1+m2
- h=0
- v=0
- mt1=m1*a1
- h, v, mm = integrate_kolenka(mm, h, v, dm=mt1, t1=mt1/s, U=U,N=60000,M=M,R=1737e3)
- mm=m+m2
- mt2=m2*a2
- h, v, mm = integrate_kolenka(mm, h, v, dm=mt2, t1=mt2/s, U=U,N=60000,M=M,R=1737e3)
- print(f'Ракета {m1:.0f}+{m2:.0f}+{m:.0f}: {h= :7.1f} {v= :7.3f} {mm= :7.1f} {max_h(h,v)}')
- exit(0)
- def full_shema(m1_arr,m2_arr,m,a,s,U,M):
- a1=a2=a
- def mass(x,y):
- m1=m*x*y
- m2=m*x*(1-y)
- return m1,m2
- #Первая ступень и вторая
- def f(x,y):
- m1,m2=mass(x,y)
- mm,h,v=m1+m2+m,0,0
- mt1=m1*a1 # Масса топлива первой ступени
- h, v, mm = integrate_kolenka(mm, h, v, dm=mt1, t1=mt1/s, U=U,N=60000,M=M,R=1737e3)
- mm=m+m2
- mt2=m2*a2 # Масса топлива второй ступени
- #print(f'{mt2=:.1f}')
- h, v, mm = integrate_kolenka(mm, h, v, dm=mt2, t1=mt2/s, U=U,N=20000,M=M,R=1737e3)
- e = v**2/2 #Энергия кинетическая
- e_h_max = G*M*(1/(h+R))
- if e>e_h_max:
- # 2*e=v**2
- return ((e-e_h_max)*2)**0.5
- else:
- return -((e_h_max-e)*2)**0.5
- #print(m1_arr)
- xx, yy = np.meshgrid(m1_arr,m2_arr)
- z=np.zeros_like(xx)
- zmax=-124124
- zix,ziy=0,0
- for iy in range(xx.shape[0]):
- for ix in range(yy.shape[1]):
- x=xx[iy,ix]
- y=yy[iy,ix]
- z[iy,ix]=f(x,y)
- if z[iy,ix]>zmax:
- zmax=z[iy,ix]
- zix,ziy=ix,iy
- #m2=yy[ziy,zix]*m
- #m1=xx[ziy,zix]*(m2+m)
- m1, m2 = mass(xx[ziy,zix], yy[ziy,zix])
- print(f"{m=:.1f} {m2=:.1f} {m1=:.1f} {m1/m2=:.2f} {m2/m=:.2f} Остаток скорости: {z[ziy,zix]:.1f}")
- fig,ax=plt.subplots(1, 1)
- extent = np.min(xx), np.max(xx), np.max(yy), np.min(yy)
- def mur(ax, tar, ss=None):
- CS = ax.contour(xx, yy, tar, linestyles=':') #, levels=np.arange(0, +150, 10)
- ax.clabel(CS)
- ax.set_title(ss)
- c = ax.imshow(tar, cmap='RdBu', interpolation='bicubic',extent=extent)
- dy, dx = np.gradient(tar)
- ax.streamplot(xx, yy, dx, dy, color=(0.7, 0.7, 0.7), linewidth=1.0, density=0.5)
- ax.set_aspect("auto")
- ax.set_xlabel("Во сколько раз ракета тяжелее спутника")
- ax.set_ylabel("Процент топлива в первой ступени")
- ax.axis([xx.min(), xx.max(), yy.min(), yy.max()])
- fig.colorbar(c, ax=ax)
- mur(ax, z, "Скорость на бесконечности: ")
- plt.show()
- #0
- full_shema(np.linspace(3,35,90),np.linspace(0.5,0.9,100),m=1000,a=0.6,s=40,U=3500,M=M)
Add Comment
Please, Sign In to add comment