Guest User

Untitled

a guest
Apr 9th, 2021
99
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 8.83 KB | None | 0 0
  1.  
  2. import numpy as np
  3. from time import *
  4. from numba import njit as njit
  5. import numba.types
  6. from math import *
  7. import matplotlib.pyplot as plt
  8.  
  9. test=False
  10.  
  11. G=6.67e-11 #Постоянная
  12. #Луна
  13. M=7.3477e22
  14. R=1737e3
  15. @njit
  16. def integrate_kolenka(m1,h1,v1,dm,t1,U,N=100000,M=7.3477e22,R=1737e3):
  17.     #                                             Параметры Луны по умолчанию вобью
  18.     #dm - масса топлив
  19.     #t1 - время сгорания
  20.     #U - скорость истечения газов
  21.     h=h1
  22.     v=v1
  23.     for c in range(N):
  24.         t=t1*c/N
  25.         dt=t1/N
  26.         ev=U*dm/(m1*t1-t*dm)-G*M/(R+h)**2
  27.         eh=v
  28.         v+=ev*dt
  29.         h+=eh*dt
  30.         #if (c==0):
  31.         #    print(f'{ev=} {dt=}')
  32.         if h<0: #Если ракета не взлетает
  33.             h=0
  34.             if (v>1):#Если слишком жёсткое падение
  35.                 print("Ракета упала со скоростью: ", abs(v))
  36.             v=max(v,0)
  37.     v2=v
  38.     h2=h
  39.     return h2,v2,m1-dm
  40.  
  41. if test:
  42.     print("-----------------------------")
  43.     #Проверим интегратор в невесомости с помощью формулы Цилковского про дельту
  44.     m=10000
  45.     v1=0
  46.     h1=0
  47.     h2,v2,m2=integrate_kolenka(m,h1,v1,dm=7000,t1=20,U=800,M=0)
  48.     v_c2=800*log(10000/3000)
  49.     print(f'{h2=:.1f} {v2=:.3f} {m2=:.1f}     Цикловкский: v={v_c2:.3f}')
  50.  
  51.     h3,v3,m3=integrate_kolenka(m2,h2,v2,dm=2000,t1=20,U=1200,M=0)
  52.     v_c3=v_c2+1200*log(3000/1000)
  53.     print(f'{h3=:.1f} {v3=:.3f} {m3=:.1f}     Цикловкский: v={v_c3:.3f}')
  54.  
  55.     # Выводит это:
  56.     # h2=7744.0 v2=963.172 m2=3000.0     Цикловкский: v=963.178
  57.     # h3=37824.0 v3=2281.498 m3=1000.0     Цикловкский: v=2281.513
  58.     # Что меня вполне устраивает
  59.  
  60.  
  61.  
  62.  
  63.  
  64.  
  65.  
  66. def max_h(h,v):
  67.     e=v**2/2 #Энергия кинетическая
  68.     e_h_max= G*M*(1/R)
  69.     if e>e_h_max:
  70.         # 2*e=v**2
  71.         v=((e-e_h_max)*2)**0.5
  72.         return f"Улетит на бесконечность со скорость {v:.1f}"
  73.     else:
  74.         if (v):
  75.             h_max=(R+h)/((2*G*M/((R+h)*v**2))-1)
  76.             return f"Улетит на высоту {h_max:.1f}"
  77.         else:
  78.             return f"Будет падать с высоты {h:.1f}"
  79.  
  80.  
  81. if test:
  82.     print("-----------------------------")
  83.     #Проверю формулу
  84.     for v in [1,10,100,1000,2350,2400,6000]:
  85.         #2380 - вторая космическая на луне
  86.         g=M*G/R/R #Ускорение свободного падения на луне
  87.         h_max=v*v/(2*g) #Высота при постоянно g
  88.         print(f'{v=:.1f}   {max_h(0,v)}   по формуле: {h_max:.1f}')
  89.     print(f'Остаток скорости при 6000:  {(6000**2-2380**2)**0.5:.1f}')
  90.  
  91.     """   Вывод
  92.    v=1.0   Улетит на высоту 0.3   по формуле: 0.3
  93.    v=10.0   Улетит на высоту 30.8   по формуле: 30.8
  94.    v=100.0   Улетит на высоту 3083.6   по формуле: 3078.2
  95.    v=1000.0   Улетит на высоту 374114.3   по формуле: 307816.9
  96.    v=2350.0   Улетит на высоту 79629015.7   по формуле: 1699918.6
  97.    v=2400.0   Улетит на бесконечность со скорость 342.1   по формуле: 1773025.1
  98.    v=6000.0   Улетит на бесконечность со скорость 5509.7   по формуле: 11081406.6
  99.    Остаток скорости при 6000:  5507.8"""
  100.  
  101.     # Вроде верно, на низких высотах совпадает, при приближении ко второй космической
  102.     #  высота сильно отклоняется, а потом уходит на бесконечность
  103.  
  104.  
  105.  
  106. if test:
  107.     print("-----------------------------")
  108.     #Проверим мою гипотезу, что время сгорания существенно для итогового результата
  109.     m=10000
  110.     v1=0
  111.     h1=0
  112.     for t in [1,10,100,980,1000,1020,10000]:
  113.         h2,v2,m2=integrate_kolenka(m,h1,v1,dm=7000,t1=t,U=800)
  114.         print(f'Сгорание за {t:5} секунд:   {h2= :7.1f}   {v2= :7.3f}  {m2= :7.1f}      {max_h(h2,v2)}')
  115.  
  116.  
  117.     """ У меня выводит такое
  118.    Сгорание за     1 секунд:   h2=   386.4   v2= 961.548  m2=  3000.0      Улетит на высоту 340533.2
  119.    Сгорание за    10 секунд:   h2=  3790.9   v2= 946.950  m2=  3000.0      Улетит на высоту 329742.9
  120.    Сгорание за   100 секунд:   h2= 30637.5   v2= 802.396  m2=  3000.0      Улетит на высоту 232197.6
  121.    Сгорание за   980 секунд:   h2=   226.5   v2=   9.633  m2=  3000.0      Улетит на высоту 28.6
  122.    Сгорание за  1000 секунд:   h2=   153.8   v2=   7.391  m2=  3000.0      Улетит на высоту 16.8
  123.    Сгорание за  1020 секунд:   h2=    98.9   v2=   5.470  m2=  3000.0      Улетит на высоту 9.2
  124.    Сгорание за 10000 секунд:   h2=     0.0   v2=   0.000  m2=  3000.0      Будет падать с высоты 0.0"""
  125.  
  126.     #При сгорании за 1-10 секунд ничего не меняется, при 100 скорость меньше, но успевает набраться высота
  127.     #При сгорании за 1000 секунд едва отрывается (это случайно получилось,
  128.     #                                             что при таких параметрах оно едва-едва отрывается)
  129.     #При сгорании за 10000 секунд не взлетает
  130.  
  131. if False:
  132.     for t in np.linspace(0.1,0.9,20):
  133.         m=1000
  134.         c1=20
  135.         s=40
  136.         a1=a2=0.6
  137.         U=3326
  138.         m1=m*c1*t
  139.         m2=m*c1*(1-t)
  140.         mm=m+m1+m2
  141.         h=0
  142.         v=0
  143.         mt1=m1*a1
  144.         h, v, mm = integrate_kolenka(mm, h, v, dm=mt1, t1=mt1/s, U=U,N=60000,M=M,R=1737e3)
  145.         mm=m+m2
  146.         mt2=m2*a2
  147.         h, v, mm = integrate_kolenka(mm, h, v, dm=mt2, t1=mt2/s, U=U,N=60000,M=M,R=1737e3)
  148.         print(f'Ракета {m1:.0f}+{m2:.0f}+{m:.0f}:   {h= :7.1f}   {v= :7.3f}  {mm= :7.1f}      {max_h(h,v)}')
  149.     exit(0)
  150.  
  151.  
  152. def full_shema(m1_arr,m2_arr,m,a,s,U,M):
  153.     a1=a2=a
  154.  
  155.     def mass(x,y):
  156.         m1=m*x*y
  157.         m2=m*x*(1-y)
  158.         return m1,m2
  159.     #Первая ступень и вторая
  160.     def f(x,y):
  161.         m1,m2=mass(x,y)
  162.         mm,h,v=m1+m2+m,0,0
  163.         mt1=m1*a1 # Масса топлива первой ступени
  164.         h, v, mm = integrate_kolenka(mm, h, v, dm=mt1, t1=mt1/s, U=U,N=60000,M=M,R=1737e3)
  165.         mm=m+m2
  166.         mt2=m2*a2 # Масса топлива второй ступени
  167.         #print(f'{mt2=:.1f}')
  168.         h, v, mm = integrate_kolenka(mm, h, v, dm=mt2, t1=mt2/s, U=U,N=20000,M=M,R=1737e3)
  169.         e = v**2/2  #Энергия кинетическая
  170.         e_h_max = G*M*(1/(h+R))
  171.         if e>e_h_max:
  172.             # 2*e=v**2
  173.             return ((e-e_h_max)*2)**0.5
  174.         else:
  175.             return -((e_h_max-e)*2)**0.5
  176.  
  177.     #print(m1_arr)
  178.     xx, yy = np.meshgrid(m1_arr,m2_arr)
  179.     z=np.zeros_like(xx)
  180.  
  181.     zmax=-124124
  182.     zix,ziy=0,0
  183.     for iy in range(xx.shape[0]):
  184.         for ix in range(yy.shape[1]):
  185.                 x=xx[iy,ix]
  186.                 y=yy[iy,ix]
  187.                 z[iy,ix]=f(x,y)
  188.                 if z[iy,ix]>zmax:
  189.                     zmax=z[iy,ix]
  190.                     zix,ziy=ix,iy
  191.     #m2=yy[ziy,zix]*m
  192.     #m1=xx[ziy,zix]*(m2+m)
  193.  
  194.     m1, m2 = mass(xx[ziy,zix], yy[ziy,zix])
  195.     print(f"{m=:.1f} {m2=:.1f} {m1=:.1f}  {m1/m2=:.2f}  {m2/m=:.2f}  Остаток скорости: {z[ziy,zix]:.1f}")
  196.  
  197.  
  198.     fig,ax=plt.subplots(1, 1)
  199.  
  200.     extent = np.min(xx), np.max(xx), np.max(yy), np.min(yy)
  201.     def mur(ax, tar, ss=None):
  202.         CS = ax.contour(xx, yy, tar, linestyles=':') #, levels=np.arange(0, +150, 10)
  203.         ax.clabel(CS)
  204.         ax.set_title(ss)
  205.         c = ax.imshow(tar, cmap='RdBu', interpolation='bicubic',extent=extent)
  206.         dy, dx = np.gradient(tar)
  207.         ax.streamplot(xx, yy, dx, dy, color=(0.7, 0.7, 0.7), linewidth=1.0, density=0.5)
  208.         ax.set_aspect("auto")
  209.         ax.set_xlabel("Во сколько раз ракета тяжелее спутника")
  210.         ax.set_ylabel("Процент топлива в первой ступени")
  211.         ax.axis([xx.min(), xx.max(), yy.min(), yy.max()])
  212.         fig.colorbar(c, ax=ax)
  213.     mur(ax, z, "Скорость на бесконечности: ")
  214.     plt.show()
  215.  
  216.  
  217.  
  218. #0
  219. 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)
  220.  
Add Comment
Please, Sign In to add comment