VERSION 5.00
Begin VB.Form Form2 
   AutoRedraw      =   -1  'True
   Caption         =   "Моделирование и имитирование движения планет    Solsys7 - форма 2"
   ClientHeight    =   8295
   ClientLeft      =   60
   ClientTop       =   315
   ClientWidth     =   11880
   FillStyle       =   0  'Заливка
   Icon            =   "Form2.frx":0000
   LinkTopic       =   "Form2"
   Picture         =   "Form2.frx":0442
   ScaleHeight     =   8295
   ScaleWidth      =   11880
   Begin VB.CheckBox Check23 
      Caption         =   "text dR"
      Height          =   195
      Left            =   6480
      TabIndex        =   125
      Top             =   7560
      Width           =   975
   End
   Begin VB.CheckBox Check18 
      Caption         =   "graf dR"
      Height          =   195
      Left            =   7440
      TabIndex        =   124
      Top             =   7560
      Width           =   855
   End
   Begin VB.CommandButton Command9 
      Caption         =   "дополнительные параметры"
      Height          =   255
      Left            =   6840
      TabIndex        =   123
      Top             =   2400
      Width           =   2415
   End
   Begin VB.CheckBox Check21 
      Caption         =   "текст теста"
      Height          =   195
      Left            =   4440
      TabIndex        =   122
      Top             =   8020
      Width           =   1215
   End
   Begin VB.CheckBox Check19 
      Caption         =   "Xzem=0"
      Height          =   195
      Left            =   7440
      TabIndex        =   121
      Top             =   4080
      Width           =   865
   End
   Begin VB.CommandButton Command8 
      Caption         =   "Закрыть форму"
      Height          =   255
      Left            =   9600
      TabIndex        =   120
      Top             =   8000
      Width           =   1455
   End
   Begin VB.CommandButton Command10 
      Caption         =   "К плану"
      Height          =   255
      Left            =   8680
      TabIndex        =   119
      Top             =   8000
      Width           =   855
   End
   Begin VB.TextBox Text53 
      Height          =   285
      Left            =   6240
      TabIndex        =   118
      Text            =   "Text53"
      Top             =   7440
      Visible         =   0   'False
      Width           =   615
   End
   Begin VB.CheckBox Check17 
      Caption         =   "Xsol=0"
      Height          =   195
      Left            =   7440
      TabIndex        =   114
      Top             =   3840
      Value           =   1  'Отмечено
      Width           =   800
   End
   Begin VB.TextBox Text51 
      Height          =   285
      Left            =   7200
      TabIndex        =   112
      Text            =   "Text51"
      Top             =   5760
      Visible         =   0   'False
      Width           =   615
   End
   Begin VB.TextBox Text50 
      Height          =   285
      Left            =   11280
      TabIndex        =   111
      Text            =   "Text50"
      Top             =   7660
      Width           =   615
   End
   Begin VB.TextBox Text49 
      Height          =   285
      Left            =   8760
      TabIndex        =   110
      Text            =   "Text49"
      Top             =   6360
      Width           =   900
   End
   Begin VB.TextBox Text48 
      Height          =   285
      Left            =   7800
      TabIndex        =   109
      Text            =   "Text48"
      Top             =   6360
      Width           =   900
   End
   Begin VB.TextBox Text47 
      Height          =   285
      Left            =   8760
      TabIndex        =   108
      Text            =   "Text47"
      Top             =   6040
      Width           =   900
   End
   Begin VB.TextBox Text46 
      Height          =   285
      Left            =   7800
      TabIndex        =   106
      Text            =   "Text46"
      Top             =   6040
      Width           =   900
   End
   Begin VB.Frame Frame4 
      Caption         =   "модель"
      Height          =   455
      Left            =   6800
      TabIndex        =   103
      Top             =   440
      Width           =   2415
      Begin VB.OptionButton Option8 
         Caption         =   "ОТО"
         Enabled         =   0   'False
         Height          =   195
         Left            =   1320
         TabIndex        =   105
         Top             =   200
         Width           =   645
      End
      Begin VB.OptionButton Option7 
         Caption         =   "Классика"
         Height          =   195
         Left            =   60
         TabIndex        =   104
         Top             =   200
         Value           =   -1  'True
         Width           =   1095
      End
   End
   Begin VB.CheckBox Check16 
      Caption         =   "dPeriod, дни"
      Height          =   195
      Left            =   9720
      TabIndex        =   102
      Top             =   6400
      Width           =   1215
   End
   Begin VB.CheckBox Check15 
      Caption         =   "dRsr,тыс.км"
      Height          =   195
      Left            =   9720
      TabIndex        =   101
      Top             =   6090
      Width           =   1215
   End
   Begin VB.TextBox Text45 
      Height          =   285
      Left            =   7840
      TabIndex        =   100
      Text            =   "2"
      Top             =   1800
      Width           =   1335
   End
   Begin VB.CheckBox Check14 
      Caption         =   "эксперимент закончить в"
      Height          =   195
      Left            =   8280
      TabIndex        =   98
      Top             =   3440
      Value           =   1  'Отмечено
      Width           =   2415
   End
   Begin VB.Frame Frame3 
      Caption         =   "имитатор"
      Height          =   455
      Left            =   6800
      TabIndex        =   94
      Top             =   860
      Width           =   2415
      Begin VB.OptionButton Option9 
         Caption         =   "JPL0"
         Height          =   195
         Left            =   880
         TabIndex        =   115
         Top             =   200
         Value           =   -1  'True
         Width           =   735
      End
      Begin VB.OptionButton Option6 
         Caption         =   "JPL2"
         Height          =   195
         Left            =   1640
         TabIndex        =   96
         Top             =   200
         Width           =   720
      End
      Begin VB.OptionButton Option5 
         Caption         =   "New0"
         Height          =   195
         Left            =   40
         TabIndex        =   95
         Top             =   200
         Width           =   720
      End
   End
   Begin VB.Frame Frame2 
      Caption         =   "использовать"
      Height          =   455
      Left            =   6800
      TabIndex        =   91
      Top             =   0
      Width           =   2415
      Begin VB.OptionButton Option4 
         Caption         =   "имитатор"
         Height          =   195
         Left            =   1200
         TabIndex        =   93
         Top             =   200
         Width           =   1080
      End
      Begin VB.OptionButton Option3 
         Caption         =   "модель"
         Height          =   195
         Left            =   60
         TabIndex        =   92
         Top             =   200
         Value           =   -1  'True
         Width           =   875
      End
   End
   Begin VB.TextBox Text44 
      Height          =   285
      Left            =   9960
      TabIndex        =   90
      Text            =   "Text44"
      Top             =   4400
      Width           =   495
   End
   Begin VB.TextBox Text43 
      Height          =   285
      Left            =   9600
      TabIndex        =   89
      Text            =   "Text43"
      Top             =   4400
      Width           =   305
   End
   Begin VB.TextBox Text42 
      Height          =   285
      Left            =   9120
      TabIndex        =   88
      Text            =   "Text42"
      Top             =   4400
      Width           =   420
   End
   Begin VB.TextBox Text14 
      Height          =   285
      Left            =   10920
      TabIndex        =   87
      Text            =   "Text14"
      Top             =   4400
      Width           =   975
   End
   Begin VB.CheckBox Check12 
      Caption         =   "dEsys, dMsys, %"
      Height          =   255
      Left            =   8300
      TabIndex        =   86
      Top             =   3760
      Width           =   1505
   End
   Begin VB.TextBox Text41 
      Height          =   285
      Left            =   10920
      TabIndex        =   85
      Text            =   "Text41"
      Top             =   3740
      Width           =   975
   End
   Begin VB.TextBox Text40 
      Height          =   285
      Left            =   10680
      TabIndex        =   84
      Text            =   "2000"
      Top             =   3400
      Width           =   615
   End
   Begin VB.CommandButton Command7 
      Caption         =   "Тест гравитации*"
      Height          =   255
      Left            =   5640
      TabIndex        =   83
      Top             =   8000
      Width           =   1605
   End
   Begin VB.Frame Frame1 
      Caption         =   "время идет"
      Height          =   455
      Left            =   6800
      TabIndex        =   80
      Top             =   1300
      Width           =   2415
      Begin VB.OptionButton Option2 
         Caption         =   "Назад"
         Height          =   195
         Left            =   1320
         TabIndex        =   82
         Top             =   200
         Width           =   855
      End
      Begin VB.OptionButton Option1 
         Caption         =   "Вперед"
         Height          =   195
         Left            =   60
         TabIndex        =   81
         Top             =   200
         Value           =   -1  'True
         Width           =   975
      End
   End
   Begin VB.TextBox Text39 
      Height          =   285
      Left            =   8880
      TabIndex        =   77
      Text            =   "1000"
      Top             =   3080
      Width           =   495
   End
   Begin VB.TextBox Text38 
      Height          =   285
      Left            =   8120
      TabIndex        =   76
      Text            =   "3600"
      Top             =   3080
      Width           =   615
   End
   Begin VB.TextBox Text37 
      Height          =   285
      Left            =   10800
      TabIndex        =   75
      Text            =   "1000"
      Top             =   3080
      Width           =   495
   End
   Begin VB.CheckBox Check11 
      Caption         =   "dEks"
      Height          =   255
      Left            =   9720
      TabIndex        =   68
      Top             =   5760
      Width           =   735
   End
   Begin VB.TextBox Text36 
      Height          =   285
      Left            =   6240
      TabIndex        =   66
      Text            =   "0"
      Top             =   7680
      Width           =   495
   End
   Begin VB.CommandButton Command3 
      Caption         =   "К статистике"
      Height          =   255
      Left            =   7300
      TabIndex        =   65
      Top             =   8000
      Width           =   1320
   End
   Begin VB.TextBox Text35 
      Height          =   285
      Left            =   7840
      TabIndex        =   64
      Text            =   "1"
      Top             =   2120
      Width           =   735
   End
   Begin VB.TextBox Text22 
      Height          =   285
      Left            =   8760
      TabIndex        =   63
      Text            =   "Text22"
      Top             =   5400
      Width           =   900
   End
   Begin VB.TextBox Text20 
      Height          =   285
      Left            =   7800
      TabIndex        =   62
      Text            =   "Text20"
      Top             =   5720
      Width           =   900
   End
   Begin VB.CheckBox Check20 
      Caption         =   "Записать параметры орбиты в файл"
      Height          =   195
      Left            =   8300
      TabIndex        =   61
      Top             =   4080
      Width           =   3140
   End
   Begin VB.TextBox Text23 
      Height          =   285
      Left            =   10920
      TabIndex        =   60
      Text            =   "Text23"
      Top             =   5080
      Width           =   975
   End
   Begin VB.TextBox Text34 
      Height          =   285
      Left            =   8280
      TabIndex        =   59
      Text            =   "Text34"
      Top             =   7680
      Width           =   615
   End
   Begin VB.TextBox Text33 
      Height          =   285
      Left            =   8880
      TabIndex        =   57
      Text            =   "Text33"
      Top             =   7000
      Width           =   495
   End
   Begin VB.TextBox Text32 
      Height          =   285
      Left            =   8520
      TabIndex        =   56
      Text            =   "Text32"
      Top             =   7000
      Width           =   305
   End
   Begin VB.TextBox Text31 
      Height          =   285
      Left            =   8040
      TabIndex        =   55
      Text            =   "Text31"
      Top             =   7000
      Width           =   420
   End
   Begin VB.CheckBox Check13 
      Caption         =   "Vgr=Vsv *"
      Height          =   195
      Left            =   6840
      TabIndex        =   53
      Top             =   2160
      Width           =   1125
   End
   Begin VB.TextBox Text30 
      Height          =   285
      Left            =   10560
      TabIndex        =   52
      Text            =   "Text30"
      Top             =   7320
      Width           =   1335
   End
   Begin VB.TextBox Text29 
      Height          =   285
      Left            =   10920
      TabIndex        =   51
      Text            =   "Text29"
      Top             =   5720
      Width           =   975
   End
   Begin VB.CheckBox Check10 
      Caption         =   "dBetta, сек"
      Height          =   195
      Left            =   9720
      TabIndex        =   50
      Top             =   5440
      Width           =   1140
   End
   Begin VB.CheckBox Check9 
      Caption         =   "dAlfaU, сек"
      Height          =   255
      Left            =   9720
      TabIndex        =   49
      Top             =   5120
      Width           =   1140
   End
   Begin VB.CheckBox Check8 
      Caption         =   "dAlfaP, сек"
      Height          =   255
      Left            =   9720
      TabIndex        =   48
      Top             =   4800
      Width           =   1125
   End
   Begin VB.TextBox Text28 
      Height          =   285
      Left            =   8880
      TabIndex        =   45
      Text            =   "Text28"
      Top             =   6720
      Width           =   495
   End
   Begin VB.TextBox Text27 
      Height          =   285
      Left            =   8520
      TabIndex        =   44
      Text            =   "Text27"
      Top             =   6720
      Width           =   305
   End
   Begin VB.TextBox Text26 
      Height          =   285
      Left            =   8040
      TabIndex        =   43
      Text            =   "Text26"
      Top             =   6720
      Width           =   420
   End
   Begin VB.TextBox Text25 
      Height          =   285
      Left            =   10920
      TabIndex        =   42
      Text            =   "Text25"
      Top             =   4760
      Width           =   975
   End
   Begin VB.TextBox Text24 
      Height          =   285
      Left            =   8760
      TabIndex        =   41
      Text            =   "Text24"
      Top             =   5080
      Width           =   900
   End
   Begin VB.CommandButton Command2 
      Caption         =   "NC0 шагов"
      Height          =   255
      Left            =   3240
      TabIndex        =   40
      Top             =   8000
      Width           =   1065
   End
   Begin VB.CheckBox Check7 
      Caption         =   "Z, тыс.км"
      Height          =   195
      Left            =   9480
      TabIndex        =   39
      Top             =   7330
      Width           =   1095
   End
   Begin VB.CheckBox Check6 
      Caption         =   "Y, тыс.км"
      Height          =   195
      Left            =   9480
      TabIndex        =   38
      Top             =   7030
      Width           =   1095
   End
   Begin VB.CheckBox Check5 
      Caption         =   "X, тыс.км"
      Height          =   255
      Left            =   9480
      TabIndex        =   37
      Top             =   6700
      Width           =   1095
   End
   Begin VB.TextBox Text21 
      Height          =   285
      Left            =   3280
      TabIndex        =   35
      Text            =   "10000"
      Top             =   7680
      Width           =   735
   End
   Begin VB.CheckBox Check4 
      Caption         =   "Найти показатели планеты №"
      Height          =   255
      Left            =   6840
      TabIndex        =   34
      Top             =   2790
      Value           =   1  'Отмечено
      Width           =   2715
   End
   Begin VB.CheckBox Check3 
      Caption         =   "Анимация"
      Height          =   195
      Left            =   6840
      TabIndex        =   33
      Top             =   3360
      Value           =   1  'Отмечено
      Width           =   1060
   End
   Begin VB.TextBox Text19 
      Height          =   285
      Left            =   9840
      TabIndex        =   32
      Text            =   "Text19"
      Top             =   3740
      Width           =   975
   End
   Begin VB.TextBox Text18 
      Height          =   285
      Left            =   8760
      TabIndex        =   29
      Text            =   "Text18"
      Top             =   5720
      Width           =   900
   End
   Begin VB.TextBox Text17 
      Height          =   285
      Left            =   11400
      TabIndex        =   28
      Text            =   "1"
      Top             =   2760
      Width           =   495
   End
   Begin VB.TextBox Text16 
      Height          =   285
      Left            =   10920
      TabIndex        =   27
      Text            =   "Text16"
      Top             =   5400
      Width           =   975
   End
   Begin VB.TextBox Text15 
      Height          =   285
      Left            =   10920
      TabIndex        =   25
      Text            =   "Text15"
      Top             =   6360
      Width           =   975
   End
   Begin VB.TextBox Text13 
      Height          =   285
      Left            =   9560
      TabIndex        =   24
      Text            =   "1"
      Top             =   2760
      Width           =   375
   End
   Begin VB.TextBox Text12 
      Height          =   285
      Left            =   10560
      TabIndex        =   23
      Text            =   "Text12"
      Top             =   7000
      Width           =   1335
   End
   Begin VB.TextBox Text11 
      Height          =   285
      Left            =   10560
      TabIndex        =   22
      Text            =   "Text11"
      Top             =   6680
      Width           =   1335
   End
   Begin VB.TextBox Text10 
      Height          =   285
      Left            =   10920
      TabIndex        =   21
      Text            =   "Text10"
      Top             =   6040
      Width           =   975
   End
   Begin VB.CheckBox Check2 
      Caption         =   "Записать координаты  X, Y, Z"
      Height          =   195
      Left            =   6840
      TabIndex        =   20
      Top             =   7340
      Width           =   2625
   End
   Begin VB.TextBox Text9 
      Height          =   285
      Left            =   5040
      TabIndex        =   18
      Text            =   "100"
      Top             =   7680
      Width           =   495
   End
   Begin VB.TextBox Text8 
      Height          =   285
      Left            =   11410
      TabIndex        =   17
      Text            =   "1"
      Top             =   4060
      Width           =   460
   End
   Begin VB.TextBox Text7 
      Height          =   285
      Left            =   7800
      TabIndex        =   16
      Text            =   "Text7"
      Top             =   5080
      Width           =   900
   End
   Begin VB.TextBox Text6 
      Height          =   285
      Left            =   7800
      TabIndex        =   14
      Text            =   "Text6"
      Top             =   5400
      Width           =   900
   End
   Begin VB.TextBox Text5 
      Height          =   285
      Left            =   1280
      TabIndex        =   13
      Text            =   "100"
      Top             =   7680
      Width           =   735
   End
   Begin VB.TextBox Text4 
      Height          =   285
      Left            =   8760
      TabIndex        =   12
      Text            =   "Text4"
      Top             =   4760
      Width           =   900
   End
   Begin VB.TextBox Text3 
      Height          =   285
      Left            =   7800
      TabIndex        =   11
      Text            =   "Text3"
      Top             =   4760
      Width           =   900
   End
   Begin VB.HScrollBar HScroll1 
      Height          =   195
      Left            =   0
      Max             =   10000
      Min             =   1
      SmallChange     =   100
      TabIndex        =   10
      Top             =   7440
      Value           =   3000
      Width           =   3495
   End
   Begin VB.CheckBox Check1 
      Caption         =   "Траектории "
      Height          =   195
      Left            =   6840
      TabIndex        =   9
      Top             =   3600
      Width           =   1335
   End
   Begin VB.TextBox Text2 
      Height          =   285
      Left            =   8040
      TabIndex        =   8
      Text            =   "Text2"
      Top             =   4400
      Width           =   615
   End
   Begin VB.TextBox Text1 
      Height          =   285
      Left            =   9840
      TabIndex        =   6
      Text            =   "Text1"
      Top             =   7680
      Width           =   565
   End
   Begin VB.PictureBox Picture2 
      AutoRedraw      =   -1  'True
      BackColor       =   &H00000000&
      BeginProperty Font 
         Name            =   "Arial"
         Size            =   8.25
         Charset         =   204
         Weight          =   400
         Underline       =   0   'False
         Italic          =   0   'False
         Strikethrough   =   0   'False
      EndProperty
      ForeColor       =   &H80000005&
      Height          =   2670
      Left            =   9240
      ScaleHeight     =   174
      ScaleMode       =   3  'Пиксель
      ScaleWidth      =   174
      TabIndex        =   5
      Top             =   30
      Width           =   2670
   End
   Begin VB.CommandButton Command6 
      Caption         =   "Выход"
      Height          =   255
      Left            =   11120
      TabIndex        =   4
      Top             =   8000
      Width           =   760
   End
   Begin VB.CommandButton Command5 
      Caption         =   "Продолжить"
      Height          =   255
      Left            =   2040
      TabIndex        =   3
      Top             =   8000
      Width           =   1120
   End
   Begin VB.CommandButton Command4 
      Caption         =   "Остановить"
      Height          =   255
      Left            =   840
      TabIndex        =   2
      Top             =   8000
      Width           =   1080
   End
   Begin VB.CommandButton Command1 
      Caption         =   "Начать"
      Height          =   255
      Left            =   0
      TabIndex        =   1
      Top             =   8000
      Width           =   735
   End
   Begin VB.PictureBox Picture1 
      AutoRedraw      =   -1  'True
      BackColor       =   &H00000000&
      ForeColor       =   &H80000005&
      Height          =   7425
      Left            =   30
      ScaleHeight     =   7365
      ScaleWidth      =   6705
      TabIndex        =   0
      Top             =   0
      Width           =   6765
   End
   Begin VB.Label Label24 
      Caption         =   "дата"
      Height          =   255
      Left            =   8720
      TabIndex        =   117
      Top             =   4440
      Width           =   375
   End
   Begin VB.Label Label23 
      Caption         =   "dTreal, час"
      Height          =   255
      Left            =   10440
      TabIndex        =   116
      Top             =   7720
      Width           =   855
   End
   Begin VB.Label Label20 
      Caption         =   "году"
      Height          =   255
      Left            =   11480
      TabIndex        =   113
      Top             =   3440
      Width           =   495
   End
   Begin VB.Label Label16 
      Caption         =   "шагов"
      Height          =   255
      Left            =   11400
      TabIndex        =   107
      Top             =   3120
      Width           =   495
   End
   Begin VB.Label Label21 
      Caption         =   "F=G*m*M/R^"
      Height          =   255
      Left            =   6840
      TabIndex        =   99
      Top             =   1860
      Width           =   975
   End
   Begin VB.Label Label28 
      Caption         =   "JD="
      Height          =   255
      Left            =   10560
      TabIndex        =   97
      Top             =   4440
      Width           =   255
   End
   Begin VB.Label Label18 
      Caption         =   "Шаг решения, с"
      Height          =   255
      Left            =   6820
      TabIndex        =   79
      Top             =   3120
      Width           =   1335
   End
   Begin VB.Label Label15 
      Caption         =   "/"
      Height          =   255
      Left            =   8760
      TabIndex        =   78
      Top             =   3120
      Width           =   135
   End
   Begin VB.Label Label3 
      Caption         =   "с выводом через"
      Height          =   255
      Left            =   9440
      TabIndex        =   74
      Top             =   3120
      Width           =   1335
   End
   Begin VB.Label Label12 
      Caption         =   "AlfaEklip, град"
      Height          =   255
      Left            =   6840
      TabIndex        =   73
      Top             =   4440
      Width           =   1095
   End
   Begin VB.Label Label10 
      Caption         =   "Eks"
      Height          =   255
      Left            =   6840
      TabIndex        =   72
      Top             =   5745
      Width           =   495
   End
   Begin VB.Label Label7 
      Caption         =   "Betta, град"
      Height          =   255
      Left            =   6840
      TabIndex        =   71
      Top             =   5415
      Width           =   855
   End
   Begin VB.Label Label6 
      Caption         =   "AlfaU, град"
      Height          =   255
      Left            =   6840
      TabIndex        =   70
      Top             =   5115
      Width           =   855
   End
   Begin VB.Label Label5 
      Caption         =   "AlfaP, град"
      Height          =   255
      Left            =   6840
      TabIndex        =   69
      Top             =   4800
      Width           =   855
   End
   Begin VB.Label Label4 
      Caption         =   "dY, см"
      Height          =   255
      Left            =   5640
      TabIndex        =   67
      Top             =   7720
      Width           =   495
   End
   Begin VB.Label Label14 
      Caption         =   "К-во оборотов, шт"
      Height          =   255
      Left            =   6840
      TabIndex        =   58
      Top             =   7720
      Width           =   1455
   End
   Begin VB.Label Label13 
      Caption         =   "AlfaU пройден"
      Height          =   255
      Left            =   6840
      TabIndex        =   54
      Top             =   7020
      Width           =   1215
   End
   Begin VB.Image Image1 
      Height          =   285
      Left            =   6795
      Picture         =   "Form2.frx":0754
      Top             =   3920
      Width           =   540
   End
   Begin VB.Label Label27 
      Caption         =   "AlfaP пройден"
      Height          =   255
      Left            =   6840
      TabIndex        =   47
      Top             =   6720
      Width           =   1215
   End
   Begin VB.Label Label26 
      Caption         =   " dTmod, лет"
      Height          =   255
      Left            =   8880
      TabIndex        =   46
      Top             =   7720
      Width           =   975
   End
   Begin VB.Label Label22 
      Caption         =   "Mx, тыс.км/см "
      Height          =   255
      Left            =   2080
      TabIndex        =   36
      Top             =   7720
      Width           =   1215
   End
   Begin VB.Label Label19 
      Caption         =   "Period, дни"
      Height          =   255
      Left            =   6840
      TabIndex        =   31
      Top             =   6410
      Width           =   960
   End
   Begin VB.Label Label8 
      Caption         =   "Rsr, тыс.км."
      Height          =   255
      Left            =   6840
      TabIndex        =   30
      Top             =   6105
      Width           =   975
   End
   Begin VB.Label Label17 
      Caption         =   "при +/- Alfa, град="
      Height          =   255
      Left            =   9960
      TabIndex        =   26
      Top             =   2790
      Width           =   1455
   End
   Begin VB.Label Label11 
      Caption         =   "MT, лет/см"
      Height          =   255
      Left            =   4080
      TabIndex        =   19
      Top             =   7720
      Width           =   975
   End
   Begin VB.Label Label2 
      Caption         =   "MdAlfa, сек/см"
      Height          =   255
      Left            =   80
      TabIndex        =   15
      Top             =   7725
      Width           =   1215
   End
   Begin VB.Label Label1 
      Caption         =   "Масштаб орбит = 300 млн. км/см"
      Height          =   255
      Left            =   3600
      TabIndex        =   7
      Top             =   7440
      Width           =   2895
   End
End
Attribute VB_Name = "Form2"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Option Explicit

Private Sub Form_Load()
GR = Text36.Text: MLmax = 1.8 * MK * MK
Picture2.Scale (-3.9, 3.9)-(3.9, -3.9)
For i = 0 To Ne
X(i) = X0(i): Y(i) = Y0(i): Z(i) = Z0(i)
Xold(i) = X(i): Yold(i) = Y(i): Zold(i) = Z(i)
Next i
Risunok21
Risunok22
Text14.Text = JD02
Call JD_DMY(JD02) '                      по JDut находим дату для времени UT
Text44.Text = IGnew: Text43.Text = IMnew: Text42.Text = IDnew
N1 = 1: N2 = 2 'это номера планет для теста гравитации
kodPgr = 0: kodDX = 0 'не учитывать давление гравитации и не учитывать запаздывание по координатам
kod4diff = 0 'не уточнятьть dV и dR в Diffur
'kodD3R = 0 ' не уточнять dR/dt, d2R/dt2, dVi, dVj и не определять в param d3R/dt3 и d2V/dt
kPgr1 = 0: kPgr2 = 0 'так давление гравитации не будет учитываться в любом случае
End Sub

Private Sub Form_Activate() 'при активизации формы2 масштаб ML может не совпасть с масштабом ML на форме1 и по этому надо выяснить какой масштаб сейчас на этой форме
'kodDMY = 0: HourNewDay = 0 'Григорианский календарь от 0 часов
Xmouse = 0: Ymouse = 0
ML = HScroll1.Value '           Min.Value=1   Max.Value=10000
ML = 100000 * ML '         переводим в километры в сантиметре
ML = 1000 * ML '           переводим в метры в сантиметре
kodDMY = 0: HourNewDay = 0 'все делать по Григорианскому календарю от 0 часов
End Sub

Private Sub Command9_Click() 'задать дополнительные параметры взаимодействия
Form23.Show
End Sub

'ЗДЕСЬ ФАЙЛЫ С ПАРАМЕТРАМИ И С КООРДИНАТАМИ ОТКРЫВАЮТСЯ ДЛЯ ЗАПИСИ И ТУТ ЖЕ ЗАКРЫВАЮТСЯ
'
'А ФАЙЛЫ JPL ЗАКРЫВАЮТСЯ ПРИ ПЕРЕХОДЕ НА ДРУГИЕ ФОРМЫ If kodFile = 10 Then Close NfileJPL

Private Sub Command1_Click() '      начать эксперимент
On Error GoTo 30
Ipl = Text13.Text
If Ipl < 0 Or Ipl > 10 Then
Info = "номер объекта должен быть от 0 до 10"
MsgBox Info, , "Solsys7.0": Exit Sub
End If
'kodMoon = 0 или 1 для имитаторов должен прийти с формы 1 , а при моделирование
'              Луны с формы 15 по плану эксперимента он нужен для Ipl2 в NP1max
kodEnd = 0:  kodAlfaU = 0: kodAlfaP = 0 'вообще то kodAlfaP не определен
kodRezult = 2: 'чтобы при определение 1-го перигелия планета вошла в интервал, а не попала туда в самом начале
kodModel = 0: If Option3.Value = True Then kodModel = 1 'при kodModel=1 координаты определяются по моделям, а при kodModel=0 по имитаторам
kodOTO = 0: If Option8.Value = True Then kodOTO = 1 'при kodOTO=0 классическая модель а при kodОТО=1 ОТОшная модель
kodJPL2 = 0: If Option6.Value = True Then kodJPL2 = 1 'при kodJPL2=0 использовать имитаторы New0 и JPL0, а при kodJPL=1 имитатор DE405
kodTime = 1: If Option2.Value = True Then kodTime = -1 'при kodTime = -1 время течет назад от настоящего в прошлое
kodVgr = 0: If Check13.Value = 1 Then kodVgr = 1 'при kodVgr = 1 скорость гравитации не равна бесконечности

'kodPgr = 0: If Check18.Value = 1 Then kodPgr = 1 'при kodPgr = 1 учитывать давление гравитации
'If Check13.Value = 0 And Check18.Value = 1 Then
'MsgBox "Давление гравитации не может быть задано без скорости гравитации", , "Solsys7": Exit Sub
'End If
'kodDX = 0: If Check22.Value = 1 Then kodDX = 1 'при kodDX = 1 НЕ учитывать запаздывание по координатам


P00 = Text38.Text: P0 = P00: NC0 = Text37.Text: Kp0 = Text39.Text: P0old = P0
Vgr = Vsv * Text35.Text: PStR = Text45.Text '                       а может не задавать если kodVgr = 0
MdAlfa = Text5.Text: MT = Text9.Text: Mx = Text21.Text: GR = Text36.Text
Nfile = Text8.Text ' должен быть уже задан автоматически, если эксперимент по плану
NP1 = 0: Text34.Text = NP1:  NC = 0: kColor = 1

Risunok22 'Risunok21 не включать, т.к. если идет наложение графиков, то старые сотрутся

Text14.Text = JD02: JD = JD02  'расчитано по начальной дате UT на форме 1
'   ЗДЕСЬ МЕНЯ НЕ ДОЛЖНО ИНТЕРЕСОВАТЬ ПО ЭФЕМЕРИДНОМУ ИЛИ ПО СОЛНЕЧНОМУ ВРЕМЕНИ
'задано JD, а главное, что координаты и скорости планет заданы на одно и тоже время
'(пусть оно даже отличается на несколько дней, т.к. меня интересуют параметры орбит
'планет, а они значительно измениться за несколько дней не могут - даже для Луны)
'ДАЖЕ ЕСЛИ РАССЧИТАНО ПО JD ДЛЯ ИМИТАТОРА SER0, А ДЛЯ ВСЕХ ОСТАЛЬНЫХ ПО JDet1
Call JD_DMY(JD02) '                      по JDut находим дату для времени UT
IGn = IGnew: IMn = IMnew: IDn = IDnew '      вообще то они заданыо для даты UT на форме 1,
Text44.Text = IGn: Text43.Text = IMn: Text42.Text = IDn ' но могут испортиться на форме 6
IGk = Text40.Text 'год окончания эксперимента по UT (ET)
TJD = JD  '        текущая дата по UT (ET)
Text29.Text = 0: Text25.Text = 0: Text16.Text = 0: Text23.Text = 0
T0 = 0: T1 = 0: T3 = 0 'T3 используется для записи в File_AlfaP (для обработки данных) и в File_dat (XYZ)

If Check4.Value = 0 Then GoTo 20 'если не вычислять параметры орбиты планеты
If Ipl < 0.5 Or Ipl > 10.5 Then
MsgBox "Номер планеты для вычисления параметров должен быть от 1 до 10", , "Solsys7": Exit Sub
End If

AlfaP = kUgol * YP0(1, Ipl): Text3.Text = AlfaP: Text4.Text = 0
AlfaU = kUgol * YP0(2, Ipl): Text7.Text = AlfaU: Text24.Text = 0
BettaMax = YP0(3, Ipl): Text6.Text = Fix(MK * kUgol * BettaMax) / MK: Text22.Text = 0
Eks = YP0(4, Ipl): Text20.Text = Fix(MK * Eks) / MK: Text18.Text = 0
Rsr(Ipl) = YP0(5, Ipl): Text46.Text = Rsr(Ipl) / MK: Text47.Text = 0
Lmin = Rsr(Ipl) * (1 - Eks)
Lmax = Rsr(Ipl) * (1 + Eks)

X(15) = 0: Y(15) = 0: Z(15) = 0 'это координаты начала системы координат (фиктивного центра)
Ipl0 = 15 '                         расчет расстояний и углов от начала координат
If Check17.Value = 1 Then Ipl0 = 0 'расчет расстояний и углов от центра Солнца
If Check19.Value = 1 Then Ipl0 = 3 'расчет расстояний и углов от центра Земли
DZold = Z(Ipl) - Z(Ipl0) ' чтобы не сработало условие 180 при первом шаге если DZold будет равно нулю

dAlfa = Text17.Text 'при первом пролете dAlfa увеличиваем в 10 раз потому что при Vsys=300 и Vgr=Vsv очень большое смещение
AlfaPplus = AlfaP + 10 * dAlfa: If AlfaPplus > 360 Then AlfaPplus = AlfaPplus - 360 '   это интервал углов, где
AlfaPminus = AlfaP - 10 * dAlfa: If AlfaPminus < 0 Then AlfaPminus = AlfaPminus + 360 ' будем искать перигелий
AlfaUplus = AlfaU + 10 * dAlfa: If AlfaUplus > 360 Then AlfaUplus = AlfaUplus - 360 ' а это интервал углов, где
AlfaUminus = AlfaU - 10 * dAlfa: If AlfaUminus < 0 Then AlfaUminus = AlfaUminus + 360 ' будем искать угол узла восхождения
AlfaPPplus = AlfaP + 180 + 10 * dAlfa: If AlfaPPplus > 360 Then AlfaPPplus = AlfaPPplus - 360 '   это интервал углов, где
AlfaPPminus = AlfaP + 180 - 10 * dAlfa: If AlfaPPminus > 360 Then AlfaPPminus = AlfaPPminus - 360 ' будем искать афелий
AlfaUUplus = AlfaU + 180 + 10 * dAlfa: If AlfaUUplus > 360 Then AlfaUUplus = AlfaUUplus - 360 ' а это интервал углов, где
AlfaUUminus = AlfaU + 180 - 10 * dAlfa: If AlfaUUminus > 360 Then AlfaUUminus = AlfaUUminus - 360 ' будем искать угол узла нисхождения

20: 'сюда приходим, если не вычислять параметры орбит
If kodModel = 1 Then GoTo 25 'если будет модель то переход
If kodMoon = 1 And Check12.Value = 1 Then
Info = "В имитаторах энергия системы считается только для Луны вместе с Землей" & vbCrLf
Info = Info & "      (если надо то задайте на форме 1 движение Луны и Земли совместное)"
MsgBox Info, , "Solsys7.0": Exit Sub
End If

RasMassJPL0 'задаем систему масс для имитаторов JPL0 и JPL2 (если Луна, то на 27 найдем m3 и m10)
Ipl2 = Ipl: If kodMoon = 1 Then Ipl2 = 11
If kodJPL2 = 1 Then GoTo 22 ' если имитатор JPL2 то переход
If kodMoon = 1 Or Ipl = 0 Or Ipl = 10 Then
Info = "     В имитаторах New0 и JPL0 не может быть номер объекта 0 или 10" & vbCrLf
Info = Info & "и не может быть задано на форме 1 отдельно движение Луны и Земли"
MsgBox Info, , "Solsys7": Exit Sub
End If
If Option5.Value = True Then
RasMassNew0 'задаем систему масс для имитатора New0 (если Луна, то на 27 найдем m3 и m10)
InfoModel = "NEW_0" ' имитатор Ньюкома без периодических членов
FilePlan = "\NEW\Solsys" & Ipl2 & "Nab" & ".txt" '
End If
If Option9.Value = True Then
InfoModel = "JPL_0" ' имитатор JPL без периодических членов
FilePlan = "\JPL\Solsys" & Ipl2 & "Nab" & ".txt" '
End If

22: ' имитатор JPL2 (эфемериды DE405 с периодическими членами)
If kodJPL2 = 1 Then ' если файл не открыт
If kodMoon = 0 And Ipl = 10 Then
Info = "     В имитаторе JPL2 объектом 10 может быть только Луна" & vbCrLf
Info = Info & "(если надо то задайте на форме 1 отдельно движение Луны и Земли)"
MsgBox Info, , "Solsys7": Exit Sub
End If
InfoModel = "JPL_DE405 "
FilePlan = "\JPL\Solsys" & Ipl2 & "Nab" & ".txt" 'FileJPL = 405 это теперь не надо
If kodFile = 10 Then Close NfileJPL ' если файл уже был открыт
FileJPL2 'ищем первый файл данных JPL для начальной даты и читаем блок данных
End If 'файл остается открытым (а если просто остановим программу, то закроется при выходе)
GoTo 27

25: 'вид модели с различными параметрами орбит и начальными данными
InfoModel = "КЛАССИКА  "
If kodOTO = 1 Then InfoModel = "OTO  "
If Form1.Option4.Value = True Then InfoModel = InfoModel & " NEW0"
If Form1.Option5.Value = True Then InfoModel = InfoModel & " JPL0"
If Form1.Option6.Value = True Then InfoModel = InfoModel & " FILE"
If Form1.Option7.Value = True Then InfoModel = InfoModel & " SCREEN"
If Form1.Option8.Value = True Then InfoModel = InfoModel & " SER0"
If Form1.Option9.Value = True Then InfoModel = InfoModel & " JPL2"

27: If kodMoon = 1 Then '    если изменили массы для имитаторов
m(3) = m(13): m(10) = m(11) 'то восстанавливаем их для моделирования
End If

If NPlan > 0 Then Risunok21 '        если по плану, то надо очистить экран
For i = 0 To Ne
Picture1.FillColor = RGB(0, 0, 0) 'стираем точки от предыдущего эксперимента если Ipl0 не изменилось
Picture1.Circle ((X(i) - X(Ipl0)) / ML, (Y(i) - Y(Ipl0)) / ML), 0.06, RGB(0, 0, 0)
Picture1.Circle ((X(i) - X(Ipl0)) / ML, (Z(i) - Z(Ipl0)) / ML + 11), 0.06, RGB(0, 0, 0)
Picture2.Circle (X(i) / MLmax, Y(i) / MLmax), 0.04, RGB(0, 0, 0)
'координаты и скорости должны быть рассчитаны для всех объектов даже если
'kodIpl(i)=0 по тому что при расчете энергии системы нужны все R0(i), т.к.
'в моделях расчет идет всегда от 1 до 10 (в ДИФФУР используется уже kodF(i,j))
X(i) = X0(i): VX(i) = kodTime * (VX0(i) + 1000 * VXsys)
Y(i) = Y0(i): VY(i) = kodTime * (VY0(i) + 1000 * VYsys)
Z(i) = Z0(i): VZ(i) = kodTime * (VZ0(i) + 1000 * VZsys)
Xold(i) = X(i): Yold(i) = Y(i): Zold(i) = Z(i)
For j = 0 To Ne ' находим начальные расстояния между планетами и задаем относительную скорость равной нулю
If i = j Then GoTo 28
DX(i, j) = X(j) - X(i)
DY(i, j) = Y(j) - Y(i)
DZ(i, j) = Z(j) - Z(i)
Rxy(i, j) = Sqr(DX(i, j) * DX(i, j) + DY(i, j) * DY(i, j)) '  находим расстояния между объектами
R(i, j) = Sqr(Rxy(i, j) * Rxy(i, j) + DZ(i, j) * DZ(i, j))
Rold(i, j) = R(i, j): dRold(i, j) = 0: d2Rold(i, j) = 0
dViOld(i, j) = 0: d2Vi(i, j) = 0: dVjOld(i, j) = 0: d2Vj(i, j) = 0
28: Next j
Next i

'AlfaGR(Ipl) = YP0(7, Ipl) * kUgol ' для определения 1-го AlfaOld, т.е. значения на предыдущем шаге решения
Call AlfaBetta(Ipl0, Ipl) 'YP0(7, Ipl)=AlfaL0 всегда с большим минусом и давало ошибку

E0sys = 0: MKD0sys = 0
Ipl2 = 9: If kodModel = 1 Then Ipl2 = 10 'в имитаторах не может быть 10-го объекта
For i = 0 To Ipl2
'For i = 0 To 10' отладка до 10
Q(i) = X(i): Q(i + Ne + 1) = Y(i): Q(i + 2 * (Ne + 1)) = Z(i)
Q(i + 3 * (Ne + 1)) = VX(i): Q(i + 4 * (Ne + 1)) = VY(i): Q(i + 5 * (Ne + 1)) = VZ(i)
DX(0, i) = X(i) - X(0): DY(0, i) = Y(i) - Y(0): DZ(0, i) = Z(i) - Z(0)
R0(i) = Sqr(DX(0, i) * DX(0, i) + DY(0, i) * DY(0, i) + DZ(0, i) * DZ(0, i))
V(i) = Sqr(VX(i) * VX(i) + VY(i) * VY(i) + VZ(i) * VZ(i))
MKD0sys = MKD0sys + kodIpl(i) * m(i) * ((VY(i) - VY(0)) * DX(0, i) - (VX(i) - VX(0)) * DY(0, i))
'MKD0sys - момент количества движения системы относительно Солнца, но только в плоскости XY
E0sys = E0sys + kodIpl(i) * m(i) * V(i) * V(i) / 2 'энергия системы вместе с Солнцем и когда за
'Picture1.FillColor = RGB(250, 250, 250) 'нулевой уровень для потенциальной энергии планет принято R0(I)
'Info1 = i & "   " & Fix(X(i)) & "     " & Fix(Y(i)) & "     " & Fix(Z(i)) & "     " & Fix(m(i))
'Picture1.CurrentX = -9: Picture1.CurrentY = 1 - i / 2: Picture1.Print Info1
'Info1 = i & "   " & Fix(VX(i)) & "     " & Fix(VY(i)) & "     " & Fix(VZ(i)) & "     " & Fix(R0(i))
'Picture1.CurrentX = -9: Picture1.CurrentY = -4.5 - i / 2: Picture1.Print Info1
Next i
If Check12.Value = 1 Then 'выводим энергию системы и момент при заданных Xi и VXi
Picture1.FillColor = RGB(250, 250, 250) 'на форме 1 и системе масс заданных здесь
Info1 = "                                               для имитаторов"
Picture1.CurrentX = 3: Picture1.CurrentY = -8: Picture1.Print Info1
Info1 = "E0sys (Дж)=              " & E0sys & "      X0i и VX0i заданны на форме 1"
Picture1.CurrentX = -9: Picture1.CurrentY = -8.5: Picture1.Print Info1
Info1 = "M0sys (кг*м^2/с) =   " & MKD0sys & "      а массы планет заданны здесь"
Picture1.CurrentX = -9: Picture1.CurrentY = -9: Picture1.Print Info1
End If


Ipl2 = Ipl: If kodMoon = 1 Then Ipl2 = 11
Period = Period0(Ipl2): Text48.Text = Period: Text49.Text = 0
'NP1max = Fix((1.1 + Vsv / Vgr) * 2 * (IGk - IGn) * 365.25 / Period)'  kodVgr
NP1max = Fix(2.5 * (IGk - IGn) * 365.25 / Period)

'а может быть массивы AlfaP и AlfaU пригодятся
'ReDim Yi(NP1max): ReDim XD(NP1max) ': ReDim Ves(2) 'это будут массивы AlfaP и AlfaU
'Ves(1) = 0: Ves(2) = 0 'это будут полные обороты AlfaP и AlfaU (обязательно для Луны)

If Check20.Value = 1 Then '              если включена запись параметров
If kodModel = 1 Then FilePlan = "\Plan" & NPlan & "\Txt" & Ipl2 & "\Solsys" & Ipl2 & "Ras" & Nfile & ".txt"
NFileP = FreeFile
Open App.Path + FilePlan For Output As #NFileP
Info2 = "эксперимент № " & Nfile & " для объекта № " & Ipl2 & " от  " & Date & " плоскость эклиптики в эпоху J2000 начало " & Format(IDn, "00.00") & "." & Format(IMn, "00") & "." & Format(IGn, "0000") & vbCrLf
Info2 = Info2 & "шаг решения " & P00 & " сек  уменьшен у перигелия в " & Kp0 & " раз  " & "   dAlfa= " & dAlfa & "   kodBary= " & kodBary & "   kodMoon= " & kodMoon & "   kodDX= " & kodDX & "   kodPgr= " & kodPgr & vbCrLf
Info2 = Info2 & InfoModel & "    VXsys=  " & VXsys & "   VYsys= " & VYsys & "   VZsys= " & VZsys & "   Vgr= " & Fix(kodVgr * 100 * Vgr / Vsv) / 100 & " Vsv " & "  k1= " & kPgr1 & "  k2= " & kPgr2 & "  m2=  " & Fix(m(0) / m(2)) & "     m5= " & Fix(1000 * m(0) / m(5)) / 1000 & vbCrLf
Info2 = Info2 & " № замера     прошло дней   перигелий     восхождение     наклон   эксцентриситет    полуось     угл.скорость"
Print #NFileP, Info2
Print #NFileP, NP1max, Fix(100 * JD0)
Close NFileP
Info = InfoModel & "      План №  " & NPlan & "       Эксперимент №  " & Nfile & "   kodDX=  " & kodDX & "   kodPgr=  " & kodPgr
Picture1.CurrentX = -9: Picture1.CurrentY = 10: Picture1.Print Info
Info = " VXSys = " & Fix(VXsys) & "    VYSys = " & Fix(VYsys) & "    VZSys = " & Fix(VZsys) & "  (км/с)   Vgr = " & Fix(kodVgr * 100 * Vgr / Vsv) / 100 & " (Vsv)" & "  k1= " & kPgr1 & "  k2= " & kPgr2
Picture1.CurrentX = -10: Picture1.CurrentY = 9: Picture1.Print Info
End If

If Check2.Value = 1 Then '              если включено   записать координаты XYZ
Info = vbCrLf & "эксперимент для объекта № " & Ipl2 & " выполнен  " & Date & " плоскость эклиптики в эпоху J2000    начало " & Format(IDold, "00") & "." & Format(IMold, "00") & "." & Format(IGold, "0000") & vbCrLf
Info = Info & "шаг решения " & P00 & " сек  вблизи перигелия и узла восхождения уменьшен в " & Kp0 & " раз  " & "   kodMoon=  " & kodMoon & "   kodDX=  " & kodDX & "   kodPgr=  " & kodPgr & vbCrLf
Info = Info & InfoModel & "             VXsys=  " & VXsys & "   VYsys= " & VYsys & "   VZsys= " & VZsys & "    Vgr= " & Fix(kodVgr * 100 * Vgr / Vsv) / 100 & " Vsv " & "  k1= " & kPgr1 & "  k2= " & kPgr2 & vbCrLf
Info = Info & "время в секундах/1000000     X, Y, Z и Rpl в тыс. км     VX, VY, VZ и Vpl в м/с    энергия системы/10^18" & vbCrLf & vbCrLf
InfoFile = "\Dat\Solsys" & Ipl2 & "XYZ" & Nfile & ".txt"
NFileDat = FreeFile
Open App.Path + InfoFile For Output As #NFileDat
Print #NFileDat, Info
Close NFileDat
End If

'             не понял это при экпериментах по плану что ли
Do While DoEvents() ' идем на 40, после обработки сообщений в очереди
GoTo 40 '             операционной системы
Loop '
30: MsgBox "Задайте правильно данные", , "Solsys7": Exit Sub
40: Timer0 = Timer ' реальное время начала эксперимента в милисекундах
param
End Sub

Private Sub param()
On Error GoTo 330
10: 'сюда возвращаемся после каждого шага с 285 метки или после NC0 шагов с
'DoEvents имеем T1 и TJD текущего времени (в имитаторах данные будут для T1,
'                                           а в модели для T2, т.е. T1 + P0

If kodModel = 0 Then GoTo 50 '                идем на имитаторы

NVgr = 0 'счетчик обращений к Diffur (при разных обращениях по разному считается dR)
If Check23.Value = 1 Then '              если вывести текст   dR   dVj-dVi
InfoM(16) = " ДО      " & NVgr & "  dVj-dVi= " & dVj(N1, N2) - dVi(N1, N2) & "    d2Vij=   " & d2Vj(N1, N2) - d2Vi(N1, N2)
InfoM(6) = " ДО       " & NVgr & "    dR=      " & dR(N1, N2) & "    d2R=   " & d2R(N1, N2)
End If

T2 = T1 '                            метод Рунге-Кутта по 4-м коэффициентам
For k = 0 To N4: O(k) = Q(k): Next k
If kodOTO = 0 Then Diffur21
'If kodOTO = 1 Then Diffur22
For k = 0 To N4: C1(k) = P0 * W(k): Q(k) = O(k) + C1(k) / 2: Next k

T2 = T1 + P0 / 2
If kodOTO = 0 Then Diffur21
'If kodOTO = 1 Then Diffur22
For k = 0 To N4: C2(k) = P0 * W(k): Q(k) = O(k) + C2(k) / 2: Next k
If kodOTO = 0 Then Diffur21
'If kodOTO = 1 Then Diffur22
For k = 0 To N4: C3(k) = P0 * W(k): Q(k) = O(k) + C3(k): Next k

T2 = T1 + P0
If kodOTO = 0 Then Diffur21
'If kodOTO = 1 Then Diffur22
For k = 0 To N4: C4(k) = P0 * W(k): Next k
For k = 0 To N4: Q(k) = O(k) + (C1(k) + 2 * C2(k) + 2 * C3(k) + C4(k)) / 6: Next k

For i = 0 To Ne
X(i) = Q(i): Y(i) = Q(i + Ne + 1): Z(i) = Q(i + 2 * (Ne + 1)) '         делаем расшифровку переменных для более
VX(i) = Q(i + 3 * (Ne + 1)): VY(i) = Q(i + 4 * (Ne + 1)): VZ(i) = Q(i + 5 * (Ne + 1))
Next i '                                                                 понятного их восприятия
'VX(Ipl) = Q(Ipl + 3 * (Ne + 1)): VY(Ipl) = Q(Ipl + 4 * (Ne + 1)): VZ(Ipl) = Q(Ipl + 5 * (Ne + 1))
Vpl = Sqr(VX(Ipl) * VX(Ipl) + VY(Ipl) * VY(Ipl) + VZ(Ipl) * VZ(Ipl)) 'скорость выбранной планеты нужна для записи XYZ


'ВЫВОДИМ ДАННЫЕ ТАКИМИ КАК ОНИ БЫЛИ В 4-м ОБРАЩЕНИИ К DIFFUR ДЛЯ СРАВНЕНИЯ С ТЕСТОМ (form7)
If Check21.Value = 1 Then ' если вывести на экран расчет давления гравитации
Risunok21 ' очищаем текст
Picture1.CurrentX = -10: Picture1.CurrentY = 11.8: Picture1.Print InfoM(1)
Picture1.CurrentX = -10: Picture1.CurrentY = 11: Picture1.Print InfoM(2)
Info = "FX12= " & Fix(FX(N1, N2)) & "     FY= " & Fix(FY(N1, N2)) & "     FZ= " & Fix(FZ(N1, N2))
Picture1.CurrentX = -10: Picture1.CurrentY = 10.2: Picture1.Print Info
Info = "VX1=    " & VX(N1) & "     VY1=    " & VY(N1) & "     VZ1=    " & VZ(N1)
Picture1.CurrentX = -10: Picture1.CurrentY = 9.4: Picture1.Print Info
Info = "Betta12=  " & Fix(1000 * (Betta(N1, N2)) * kUgol) / 1000 & "               AlfaVX1=  " & Fix(AlfaVX(N1) * kUgol) & _
"          AlfaVY1=  " & Fix(AlfaVY(N1) * kUgol) & "          BettaVZ1=  " & Fix(BettaVZ(N1) * kUgol)
Picture1.CurrentX = -10: Picture1.CurrentY = 8.6: Picture1.Print Info
Info = "Alfa12=   " & Fix(1000 * Alfa12(N1, N2) * kUgol) / 1000 & "             dAlfaVX1=  " & Fix(10 * dAlfaVXi(N1, N2) * kUgol) / 10 _
& "        dAlfaVY1=  " & Fix(10 * dAlfaVYi(N1, N2) * kUgol) / 10 & "        dBettaVZ1=  " & Fix(10 * dBettaVZi(N1, N2) * kUgol) / 10
Picture1.CurrentX = -10: Picture1.CurrentY = 7.8: Picture1.Print Info
Info = "dV1=  " & Fix(10 * dVi(N1, N2)) / 10 & "            dVX1=  " & Fix(10 * dVXi(N1, N2)) / 10 _
& "        dVY1=  " & Fix(10 * dVYi(N1, N2)) / 10 & "        dVZ1=  " & Fix(10 * dVZi(N1, N2)) / 10
Picture1.CurrentX = -10: Picture1.CurrentY = 7: Picture1.Print Info
Info = "                                dAlfaVX2=  " & Fix(10 * dAlfaVXj(N1, N2) * kUgol) / 10 _
& "        dAlfaVY2=  " & Fix(10 * dAlfaVYj(N1, N2) * kUgol) / 10 & "        dBettaVZ2=  " & Fix(10 * dBettaVZj(N1, N2) * kUgol) / 10
Picture1.CurrentX = -10: Picture1.CurrentY = 6.2: Picture1.Print Info
Info = "dV2=  " & Fix(10 * dVj(N1, N2)) / 10 & "              dVX2=  " & Fix(10 * dVXj(N1, N2)) / 10 _
& "        dVY2=  " & Fix(10 * dVYj(N1, N2)) / 10 & "        dVZ2=  " & Fix(10 * dVZj(N1, N2)) / 10
Picture1.CurrentX = -10: Picture1.CurrentY = 5.4: Picture1.Print Info
Info = "F12=   " & F(N1, N2) & "   dV2-dV1=    " & dVj(N1, N2) - dVi(N1, N2) & "     dR12=    " & dR(N1, N2)
Picture1.CurrentX = -10: Picture1.CurrentY = 4.6: Picture1.Print Info
Picture1.CurrentX = -10: Picture1.CurrentY = 3.8: Picture1.Print InfoM(3)

Picture1.CurrentX = -10: Picture1.CurrentY = 3: Picture1.Print "N1=  " & N1
Picture1.CurrentX = -10: Picture1.CurrentY = 2: Picture1.Print InfoM(5)
Picture1.CurrentX = -10: Picture1.CurrentY = 1: Picture1.Print InfoM(4)
Picture1.CurrentX = -10: Picture1.CurrentY = 0: Picture1.Print "N2=  " & N2

Picture1.CurrentX = -10: Picture1.CurrentY = -1.3: Picture1.Print Info3
Info = "F21=   " & F(N2, N1) & "   dV1-dV2=    " & dVj(N2, N1) - dVi(N2, N1) & "     dR21=    " & dR(N2, N1)
Picture1.CurrentX = -10: Picture1.CurrentY = -2.1: Picture1.Print Info
Info = "dV1=  " & Fix(10 * dVj(N2, N1)) / 10 & "                   dVX1=  " & Fix(10 * dVXj(N2, N1)) / 10 _
& "        dVY1=  " & Fix(10 * dVYj(N2, N1)) / 10 & "        dVZ1=  " & Fix(10 * dVZj(N2, N1)) / 10
Picture1.CurrentX = -10: Picture1.CurrentY = -2.9: Picture1.Print Info
Info = "                                   dAlfaVX1=  " & Fix(10 * dAlfaVXj(N2, N1) * kUgol) / 10 _
& "        dAlfaVY1=  " & Fix(10 * dAlfaVYj(N2, N1) * kUgol) / 10 & "        dBettaVZ1=  " & Fix(10 * dBettaVZj(N2, N1) * kUgol) / 10
Picture1.CurrentX = -10: Picture1.CurrentY = -3.7: Picture1.Print Info
Info = "dV2=  " & Fix(10 * dVi(N2, N1)) / 10 & "                 dVX2=  " & Fix(10 * dVXi(N2, N1)) / 10 _
& "        dVY2=  " & Fix(10 * dVYi(N2, N1)) / 10 & "        dVZ2=  " & Fix(10 * dVZi(N2, N1)) / 10
Picture1.CurrentX = -10: Picture1.CurrentY = -4.5: Picture1.Print Info
Info = "Alfa21=   " & Fix(1000 * Alfa12(N2, N1) * kUgol) / 1000 & "                 dAlfaVX2=  " & Fix(10 * dAlfaVXi(N2, N1) * kUgol) / 10 _
& "        dAlfaVY2=  " & Fix(10 * dAlfaVYi(N2, N1) * kUgol) / 10 & "        dBettaVZ2=  " & Fix(10 * dBettaVZi(N2, N1) * kUgol) / 10
Picture1.CurrentX = -10: Picture1.CurrentY = -5.3: Picture1.Print Info
Info = "Betta21=  " & Fix(1000 * (Betta(N2, N1)) * kUgol) / 1000 & "                   AlfaVX2=  " & Fix(AlfaVX(N2) * kUgol) & _
"          AlfaVY2=  " & Fix(AlfaVY(N2) * kUgol) & "          BettaVZ2=  " & Fix(BettaVZ(N2) * kUgol)
Picture1.CurrentX = -10: Picture1.CurrentY = -6.1: Picture1.Print Info
Info = "VX2=    " & VX(N2) & "     VY2=    " & VY(N2) & "     VZ2=    " & VZ(N2)
Picture1.CurrentX = -10: Picture1.CurrentY = -6.9: Picture1.Print Info
Info = "FX21= " & Fix(FX(N2, N1)) & "     FY= " & Fix(FY(N2, N1)) & "     FZ= " & Fix(FZ(N2, N1))
Picture1.CurrentX = -10: Picture1.CurrentY = -7.7: Picture1.Print Info
Picture1.CurrentX = -10: Picture1.CurrentY = -8.5: Picture1.Print Info2
Picture1.CurrentX = -10: Picture1.CurrentY = -9.3: Picture1.Print Info1
For i = 0 To Ne
X1(i + 11) = dVX(i): X2(i + 11) = dVY(i): X3(i + 11) = dVZ(i) ' СОХРАНЯЕМ ПОСЛЕ DIFFUR для ТЕСТА, но до изменения ниже
Next i
'X1(i) = X(i): X2(i) = Y(i): X3(i) = Z(i)                     ЭТО СОХРАНЕНО В DIFFUR
'X4(i) = VX(i): X5(i) = VY(i): X6(i) = VZ(i)
End If

If Check23.Value = 1 Then '              если вывести текст   dR   dVj-dVi
Risunok21 ' очищаем текст
Picture1.CurrentX = -10: Picture1.CurrentY = 9: Picture1.Print InfoM(6)
Picture1.CurrentX = -10: Picture1.CurrentY = 8: Picture1.Print InfoM(16)
Picture1.CurrentX = -10: Picture1.CurrentY = 6: Picture1.Print InfoM(7)
Picture1.CurrentX = -10: Picture1.CurrentY = 5: Picture1.Print InfoM(17)
Picture1.CurrentX = -10: Picture1.CurrentY = 3: Picture1.Print InfoM(8)
Picture1.CurrentX = -10: Picture1.CurrentY = 2: Picture1.Print InfoM(18)
Picture1.CurrentX = -10: Picture1.CurrentY = 0: Picture1.Print InfoM(9)
Picture1.CurrentX = -10: Picture1.CurrentY = -1: Picture1.Print InfoM(19)
Picture1.CurrentX = -10: Picture1.CurrentY = -3: Picture1.Print InfoM(10)
Picture1.CurrentX = -10: Picture1.CurrentY = -4: Picture1.Print InfoM(20)
End If

If kodPgr = 0 Then GoTo 100 'если НЕ учитывать давление гравитации, то не нужет и весь этот расчет dR и dV

For i = 0 To Ne
BettaVZ(i) = pi / 2: If VZ(i) < 0 Then BettaVZ(i) = -pi / 2 ' находим углы скоростей i-х и j-х тел
AlfaVX(i) = 0: If VX(i) < 0 Then AlfaVX(i) = pi
AlfaVY(i) = pi / 2: If VY(i) < 0 Then AlfaVY(i) = 3 * pi / 2
dVX(i) = W(i + 3 * (Ne + 1)):
dVY(i) = W(i + 4 * (Ne + 1)):
dVZ(i) = W(i + 5 * (Ne + 1))
Next i

For i = 0 To Ne
For j = 0 To Ne '                                     Alfa12(i, j) и Betta(i, j) это направление скорости гравитации
If i = j Then GoTo 28 '                              при ее распространении от i-того тела до j-того тела

'If kodVgr = 0 Then '               НАВЕРНОЕ НАДО СЧИТАТЬ ПО НОВОЙ ТУТ    НО НЕ ОБЯЗАТЕЛЬНО
DX(i, j) = X(j) - X(i) '                       можно T(i,j) взять из 4-го обращения к Diffur
DY(i, j) = Y(j) - Y(i)
DZ(i, j) = Z(j) - Z(i)
Rxy(i, j) = Sqr(DX(i, j) * DX(i, j) + DY(i, j) * DY(i, j)) '  находим расстояния между объектами
R(i, j) = Sqr(Rxy(i, j) * Rxy(i, j) + DZ(i, j) * DZ(i, j))
T(i, j) = R(i, j) / Vgr '            время необходимое для распространения гравитации
'End If

If kodVgr = 1 And kodDX = 1 Then  ' запаздывания по координатам учитывать
For k = 1 To 2 '               количество итераций равно 2
DX(i, j) = X(j) - X(i) + VX(i) * T(i, j) 'расстояние между следом 1-го тела воздействующим на 2-е тело по оси X
DY(i, j) = Y(j) - Y(i) + VY(i) * T(i, j)
DZ(i, j) = Z(j) - Z(i) + VZ(i) * T(i, j)
Rxy(i, j) = Sqr(DX(i, j) * DX(i, j) + DY(i, j) * DY(i, j)) ' расстояния в плоскости XY между следом 1-го тела воздействующим на 2-е тело
R(i, j) = Sqr(Rxy(i, j) * Rxy(i, j) + DZ(i, j) * DZ(i, j)) ' общее расстояния между следом 1-го тела воздействующим на 2-е тело
T(i, j) = R(i, j) / Vgr '                                    время необходимое для распространения гравитации
Next k
End If
'         УГЛЫ , РАДИУС  И dVi(i, j) ПОЛУЧИЛИСЬ С УЧЕТОМ ЗАПАЗДЫВАНИЯ А dVj(i, j) БУДУТ ДЛЯ ТЕКУЩЕГО МОМЕНТА

Betta(i, j) = Atn(DZ(i, j) / Rxy(i, j)) '                    на 2-е тело с учетом наличия i и j планет, т.е. kodF(i, j)
Alfa12(i, j) = pi / 2: If DY(i, j) < 0 Then Alfa12(i, j) = 3 * pi / 2 'угол пригодиться при расчете сил давления
If DX(i, j) = 0 Then GoTo 22
Alfa12(i, j) = Atn(DY(i, j) / DX(i, j))
22: '                        ПЕРЕВОДИМ Alfa12(i,j) В РАЗМЕРНОСТЬ ОТ 0 ДО 2*pi
If DX(i, j) > 0 Then GoTo 24 '                     переход в 1-ю и 4-ю четверти
Alfa12(i, j) = Alfa12(i, j) + pi: GoTo 26 '        2-я и 3-я четверти (во 2-ой четверти Alfa12(i,j) был с минусом)
24:  If DY(i, j) > 0 Then GoTo 26 '                значит была 1-я четверть
Alfa12(i, j) = Alfa12(i, j) + 2 * pi '             это 4-я четверть (Alfa12(i,j) был с минусом)
26:

dAlfaVXi(i, j) = AlfaVX(i) - Alfa12(i, j) '      dAlfaVXi(i, j)- угол между скоростью по оси X i-го тела и проекцией вектора
If dAlfaVXi(i, j) > pi Then dAlfaVXi(i, j) = dAlfaVXi(i, j) - 2 * pi 'скорости гравитации от i-го тела к j-му на плоскость XY
If dAlfaVXi(i, j) < -pi Then dAlfaVXi(i, j) = dAlfaVXi(i, j) + 2 * pi '
dAlfaVYi(i, j) = AlfaVY(i) - Alfa12(i, j)
If dAlfaVYi(i, j) > pi Then dAlfaVYi(i, j) = dAlfaVYi(i, j) - 2 * pi
If dAlfaVYi(i, j) < -pi Then dAlfaVYi(i, j) = dAlfaVYi(i, j) + 2 * pi
dBettaVZi(i, j) = BettaVZ(i) - Betta(i, j)

dAlfaVXj(i, j) = AlfaVX(j) - Alfa12(i, j) '      dAlfaVX(i, j)- угол между скоростью по оси X j-го тела и проекцией вектора
If dAlfaVXj(i, j) > pi Then dAlfaVXj(i, j) = dAlfaVXj(i, j) - 2 * pi 'скорости гравитации от i-го тела к j-му на плоскость XY
If dAlfaVXj(i, j) < -pi Then dAlfaVXj(i, j) = dAlfaVXj(i, j) + 2 * pi '
dAlfaVYj(i, j) = AlfaVY(j) - Alfa12(i, j)
If dAlfaVYj(i, j) > pi Then dAlfaVYj(i, j) = dAlfaVYj(i, j) - 2 * pi
If dAlfaVYj(i, j) < -pi Then dAlfaVYj(i, j) = dAlfaVYj(i, j) + 2 * pi
dBettaVZj(i, j) = BettaVZ(j) - Betta(i, j)

VXold(i, j) = VX(i): VYold(i, j) = VY(i): VZold(i, j) = VZ(i)
If kodVgr = 1 And kodDX = 1 Then  ' запаздывания по координатам учитывать
VXold(i, j) = VX(i) - dVX(i) * T(i, j)
VYold(i, j) = VY(i) - dVY(i) * T(i, j)
VZold(i, j) = VZ(i) - dVZ(i) * T(i, j)
End If

dVZi(i, j) = Cos(dBettaVZi(i, j)) * Abs(VZold(i, j)) '                        проекции скоростей i-го тела по осям координат
dVXi(i, j) = Abs(Cos(Betta(i, j))) * Cos(dAlfaVXi(i, j)) * Abs(VXold(i, j)) ' на вектор скорости гравитации от i-го тела
dVYi(i, j) = Abs(Cos(Betta(i, j))) * Cos(dAlfaVYi(i, j)) * Abs(VYold(i, j))
dVZj(i, j) = Cos(dBettaVZj(i, j)) * Abs(VZ(j)) '                        проекции скоростей j-го тела по осям координат
dVXj(i, j) = Abs(Cos(Betta(i, j))) * Cos(dAlfaVXj(i, j)) * Abs(VX(j)) ' на вектор скорости гравитации от i-го тела
dVYj(i, j) = Abs(Cos(Betta(i, j))) * Cos(dAlfaVYj(i, j)) * Abs(VY(j))
dVi(i, j) = dVZi(i, j) + dVXi(i, j) + dVYi(i, j) ' рассчитываем относительные скорости и ускорения
dVj(i, j) = dVZj(i, j) + dVXj(i, j) + dVYj(i, j)
dR(i, j) = dVj(i, j) - dVi(i, j)
28: Next j
Next i

For i = 0 To Ne ' грубый и уточненый расчет dR(i, j) и d2R(i, j) и уточненный dVi(i, j) и dVj(i, j)
For j = 0 To Ne '                                                        для 4-х обращений к Diffur
If i = j Then GoTo 40
If P0 = P0old Then 'рассчитывать новые d2Vi(i, j) и d2Vj(i, j) только если шаг решения не менялся
d2Vi(i, j) = (dVi(i, j) - dViOld(i, j)) / P0
d2Vj(i, j) = (dVj(i, j) - dVjOld(i, j)) / P0
End If
dViOld(i, j) = dVi(i, j)
dVjOld(i, j) = dVj(i, j)
dVi(i, j) = dVi(i, j) + d2Vi(i, j) * P0 / 2 'уточняем dVi(i, j) и dVj(i, j) для середины текущего интервала P0
dVj(i, j) = dVj(i, j) + d2Vj(i, j) * P0 / 2 'от начала этого интервала
40: Next j
Next i
P0old = P0

If Check23.Value = 1 Then '              если вывести текст   dR   dVj-dVi
Info = "ПОСЛЕ   " & NVgr + 1 & "    dR=      " & dR(N1, N2) & "    d2R=   " & d2R(N1, N2)
Picture1.CurrentX = -10: Picture1.CurrentY = -6: Picture1.Print Info
Info = "ПОСЛЕ   " & NVgr + 1 & "  dVj-dVi=  " & dVj(N1, N2) - dVi(N1, N2) & "    d2Vij=   " & d2Vj(N1, N2) - d2Vi(N1, N2)
Picture1.CurrentX = -10: Picture1.CurrentY = -7: Picture1.Print Info
End If

GoTo 100 '                         пропускаем расчет по имитаторам

50:
If kodJPL2 = 0 Then GoTo 70 ' если не имитатор JPL2, т.е. не эфемериды DE405
If TJD >= JDn And TJD < JDk Then GoTo 65 'если НЕ вышли за интервал считанного блока данных
'если сделать TJD<=JDk, то при TJD=JDk получается интервал на единицу больше
'максимального и надо переходить к 1-у подинтервалу данных в следующем блоке
'где первая дата равна последней дате в предыдущем блоке данных
If IBlok < IBlokK Then GoTo 61 '     если не дошли до конца открытого файла

Close NfileJPL '                      если дошли до конца открытого файла
IGfile = IGfile + 20: IBlokK = 229 '  открываем следующий по порядку файл
If IGfile = 1680 Or IGfile = 1760 Then IBlokK = 230
If IGfile = 1820 Or IGfile = 1900 Then IBlokK = 230
If IGfile = 1980 Or IGfile = 2040 Then IBlokK = 230
If IGfile = 2120 Or IGfile = 2180 Then IBlokK = 230
NfileJPL = FreeFile
kodFile = 0
Open App.Path + "\JPL\ascp" & IGfile & ".405" For Input As #NfileJPL
kodFile = 10

61: Input #NfileJPL, IBlok, L ' номер блока и количество данных в блоке
For i = 1 To NJPL - 2 Step 3
Input #NfileJPL, Xi(i), Xi(i + 1), Xi(i + 2) ' читаем блок данных
Next i
If JDn = Xi(1) Then GoTo 61 'это если в новом файле данные 1-го блока повторяют
'данные последнего блока предыдущего файла, чтобы не было повтора тех же вычислений
JDn = Xi(1): JDk = Xi(2) 'начальная и конечная дата для данных в блоке

65: If Check12.Value = 0 Then 'если не считать энергию и момент системы
Call UrRegresJPL2(0) 'находим ГЕЛИО БАРИцентрические координаты объектов по JPL2
Call UrRegresJPL2(3) 'Земля нужна чтобы получить ГЕЛИО Луны и ГЕО планеты
If Ipl < 10 Then Call UrRegresJPL2(Ipl)
If kodMoon = 1 And Ipl = 10 Then ' если это Луна находим ГЕО координаты
Call UrRegresJPL2(11) '                 ЭКЛ ГЕО Луны
X(10) = X(11) + X(3): Y(10) = Y(11) + Y(3): Z(10) = Z(11) + Z(3) 'ЭКЛ ГЕЛИО БАРИ Луны
'здесь нет смысла уточнять ГЕОбари координаты Земли, т.к. при расчете параметров Луны
'из X(10) будет вычтено X(3) и получатся ГЕО Луны X(11), что и надо (а точные ГЕЛИО Луны не нужны)
End If
End If
If Check12.Value = 1 Then 'если считать энергию и момент системы (будет Земля+Луна)
For Ipl2 = 0 To 9 'находим БАРИцентрические координаты объектов по JPL2
Call UrRegresJPL2(Ipl2) '   и скорости
Next Ipl2 'РАЗДЕЛЕНИЕ ЭНЕРГИИ ЗЕМЛИ И ЛУНЫ И ПОТОМ СУММИРОВАНИЕ НИЧЕГО НЕ ДАДУТ
End If
If kodEnd = 1 Then GoTo 350
GoTo 100

70:  '                 для аналитически полученных имитаторов
X(0) = 0: Y(0) = 0: Z(0) = 0 'координаты Солнца всегда равны нулю
'                массы заданы при начале эксперимента (m3=m13+m11 т.к. Луны нет)
For i = 1 To 9
If Option5.Value = True Then RasParamNew0 '           теория Ньюкома
If Option9.Value = True Then RasParamJPL0 '           теория JPL0
RasKoordinat
RasSkorostey
Next i

100: 'ЗАМЕНИЛ в Solsys6 КОД ОТСЮДА ДО 285 МЕТКИ (в чем было дело так и не понял)
'                ПОСТАРАЮСЬ В 7-ОЙ ВЕРСИИ СЮДА ПОМЕНЬШЕ ВНОСИТЬ ИЗМЕНЕНИЙ
'проверяем не изменил ли я точку отсчета расстояний до планеты в расчете и рисовании
Ipl0 = 15 '                         от начала координат
If Check17.Value = 1 Then Ipl0 = 0 'от центра Солнца
If Check19.Value = 1 Then Ipl0 = 3 'от центра Земли

If Check4.Value = 0 Then GoTo 285 'если не надо определять показатели движения обекта
AlfaOld = AlfaGR(Ipl) 'AlfaGR(Ipl) при предыдущем шаге до того как Rpl стало меньше Rmax, т.е. до перигелия для AlfaP и до пересечения эклиптики для AlfaU
Call AlfaBetta(Ipl0, Ipl) 'Betta1 и Rpl будут рассчитаны при любом AlfaGR(Ipl)
'Betta1 и Rpl можно и не вычислять, т.к. все эти значения уже посчитаны в Diffur2
Betta1 = Abs(Betta1 / kUgol)
If BettaMax < Betta1 Then BettaMax = Betta1 'находим максимальный угол наклона орбиты


If AlfaGR(Ipl) <= AlfaOld Then 'прохождение 0 градусов
If NP1 > 0 Then Period = (T1 - T0) / sekS  'находим период обращения за оборот
If NP1 = 2 Then Period00 = Period:
Text48.Text = Fix(MK * Period) / MK
T0 = T1  '                    будет время предыдущего прохождение 0 градусов
End If


If Sgn(DZold) = Sgn(DZ(Ipl0, Ipl)) Then GoTo 185 'угол узла не определяем
If kodAlfaU = 1 Then GoTo 185 '      если после перигелия угол уже определен
AlfaU = (AlfaGR(Ipl) + AlfaOld) / 2 'находим один из углов линии узлов как полусумму углов до пересечения эклиптики и после
If Abs(AlfaGR(Ipl) - AlfaOld) > 355 Then AlfaU = AlfaU - 180 'если одно значение больше 0, а другое меньше, т.е. 359
If Sgn(DZold) > 0 Then GoTo 182
AlfaUplus = AlfaU + dAlfa: If AlfaUplus > 360 Then AlfaUplus = AlfaUplus - 360 ' а это интервал углов, где будем
AlfaUminus = AlfaU - dAlfa: If AlfaUminus < 0 Then AlfaUminus = AlfaUminus + 360 ' искать угол восходящего узла
Call JD_DMY(JD) 'определяем дату пролета узла восхождения для времени ET (UT)
Text33.Text = IGnew: Text32.Text = IMnew: Text31.Text = IDnew
GoTo 183 '                             это был угол восходящего узла

182: '                                 а здесь будет угол нисходящего узла
AlfaUUplus = AlfaU + dAlfa: If AlfaUUplus > 360 Then AlfaUUplus = AlfaUUplus - 360 ' а это интервал углов, где будем
AlfaUUminus = AlfaU - dAlfa: If AlfaUUminus < 0 Then AlfaUUminus = AlfaUUminus + 360 ' искать угол нисходящего узла
AlfaU = AlfaU + 180 'делаем из угла нисходящего узла угол восходящего узла
If AlfaU > 360 Then AlfaU = AlfaU - 360

183:
Text24.Text = Fix(MK * AlfaU) / MK: Text23.Text = Fix(MK * 3600 * (AlfaU - AlfaU00)) / MK
kodAlfaU = 1 '       теперь если AlfaU уже определен то в интервале поиска P0=P00

185: DZold = DZ(Ipl0, Ipl) ' высота над эклиптикой при предыдущем пролете, которая будет использована при сравнение с текущим пролетом на следующем шаге решения


'                     вычисление текущего угла в перигелии или афелии
If AlfaPminus > 360 - 2 * dAlfa Then GoTo 190 'если угол близок к 0 или 360 градусам
If AlfaGR(Ipl) > AlfaPplus Or AlfaGR(Ipl) < AlfaPminus Then GoTo 194
190: If AlfaGR(Ipl) > AlfaPplus And AlfaGR(Ipl) < AlfaPminus Then GoTo 194
kodAlfaP = 1: GoTo 200 'если не ушли на 194 значит находимся в интервале угла перигелия
194:
If AlfaPPminus > 360 - 2 * dAlfa Then GoTo 196 'если угол близок к 0 или 360 градусам
If AlfaGR(Ipl) > AlfaPPplus Or AlfaGR(Ipl) < AlfaPPminus Then GoTo 240
196: If AlfaGR(Ipl) > AlfaPPplus And AlfaGR(Ipl) < AlfaPPminus Then GoTo 240
kodAlfaP = 0: '         если не ушли на 240, значит находимся в интервале угла афелия

200: If kodRezult = 2 Then GoTo 250 'переход, т.к. расчет в этой точке уже зделан
P0 = P00 / Kp0: kodRezult = 0 ' шаг решения в перигелии или афелии уменьшаем в kP0 раз
If kodAlfaP = 0 Then GoTo 204 ' переход к афелию
If RplOld < Rpl Then GoTo 206 ' если расстояние стало увеличиваться
kodRezult = 1: RplOld = Rpl: Lmin = Rpl 'если расстояние в перигелии уменьшается
GoTo 206
204: '                                   здесь афелий
If RplOld > Rpl Then GoTo 206 '          если расстояние стало уменьшаться
kodRezult = 1: RplOld = Rpl: Lmax = Rpl 'если расстояние в афелии увеличивается

206: If kodRezult = 1 Then GoTo 239 ' если до перигелия или афелия еще не дошли, т.к. Lmin > Rpl и по этому kodRezult = 1
kodRezult = 2 ' первый шаг решения после прохождения перигелия или афелия, т.к. kodRezult осался равным нулю и значит надо определить все показатели на этом шаге
AlfaP = (AlfaGR(Ipl) + AlfaOld) / 2 '          это угол в плоскости эклиптики
If Abs(AlfaGR(Ipl) - AlfaOld) > 355 Then AlfaP = AlfaP - 180 'если одно значение больше 0, а другое меньше, т.е. 359
If kodAlfaP = 1 Then GoTo 210 '                 если находимся в перигелии
AlfaP = AlfaP + 180: If AlfaP > 360 Then AlfaP = AlfaP - 360 'если находимся в афелии

210: UgolW = AlfaP - AlfaU 'угол от узла до перигелия в плоскости эклиптики, если и он будет от 0 до +180 или от 0 до -180
If UgolW > 180 Then UgolW = UgolW - 360 ' приводим угол к размерности от 0 до +/- pi
If UgolW < -180 Then UgolW = UgolW + 360
If UgolW = -90 Or UgolW = 90 Then GoTo 216 'при -90 и 90 градусах тангенс будет равен бесконечности, по этому пишем старое значение AlfaP
UgolWreal = kUgol * Abs(Atn(Tan(UgolW / kUgol) / Cos(BettaMax))) 'Abs дает  UgolWreal от 0 до 90 градусов
If Abs(UgolW) > 90 Then UgolWreal = 180 - UgolWreal 'если UgolWreal от 90 до 180 градусов
If UgolW < 0 Then UgolWreal = -UgolWreal 'если UgolW и UgolWreal отрицательные, т.е. от 0 до -180
AlfaP = AlfaU + UgolWreal 'сумма углов в двух плоскостях если перигелий над эклиптикой
If AlfaP < 0 Then AlfaP = AlfaP + 360 ' приводим угол к размерности от 0 до 360
If AlfaP > 360 Then AlfaP = AlfaP - 360
If NP1 > 0 Then GoTo 216 ' это самый первый пролет именно перигелия когда NP1 не подсчитано
AlfaP00 = AlfaP ' AlfaP00 - точка отсчета при расчете смещения перигелия
Text3.Text = Fix(MK * AlfaP) / MK: Text4.Text = Fix(MK * AlfaP) / MK: GoTo 220

216:
Rsr(Ipl) = (Lmin + Lmax) / 2 ' находим большую полуось планеты за этот пролет
Eks = 1 - Lmin / Rsr(Ipl) '    находим эксцентриситет за этот пролет

If NP1 = 1 Then 'будет выполняться до второго пролета перигелия, что обязательно для имитаторов
AlfaU00 = AlfaU: Text7.Text = Fix(MK * AlfaU) / MK ' может определится как до AlfaP
'так и после AlfaP и по этому определяется выше только если NP1>0, а в противном случае будет заданное в начале эксперимента значение
Betta00 = BettaMax: Text6.Text = Fix(MK * kUgol * BettaMax) / MK
Eks00 = Eks: Text20.Text = Fix(MK * Eks) / MK 'определяется только если NP1>0'вместе с периодом, но
Rsr00 = Rsr(Ipl): Text46.Text = Fix(Rsr(Ipl)) / MK
End If

Text4.Text = Fix(MK * AlfaP) / MK: Text25.Text = Fix(MK * 3600 * (AlfaP - AlfaP00)) / MK
Text22.Text = Fix(MK * BettaMax * kUgol) / MK: Text16.Text = Fix(MK * 3600 * (BettaMax - Betta00) * kUgol) / MK
Text18.Text = Fix(MK * Eks) / MK: Text29.Text = Fix(MK * (Eks - Eks00)) / MK
Text47.Text = Fix(Rsr(Ipl)) / MK: Text10.Text = Fix((Rsr(Ipl) - Rsr00)) / MK
Text49.Text = Fix(MK * Period) / MK: Text15.Text = Fix(MK * (Period - Period00)) / MK

If Check8.Value = 1 Then Picture1.PSet (T1 / sekG / MT - 10, 3600 * (AlfaP - AlfaP00) / MdAlfa + GR), RGB(ColorR(Ipl), ColorG(Ipl), ColorB(Ipl))
If Check9.Value = 1 Then Picture1.PSet (T1 / sekG / MT - 10, 3600 * (AlfaU - AlfaU00) / MdAlfa + GR), RGB(ColorR(Ipl), ColorG(Ipl), ColorB(Ipl))
If Check10.Value = 1 Then Picture1.PSet (T1 / sekG / MT - 10, 3600 * kUgol * (BettaMax - Betta00) / MdAlfa + GR), RGB(ColorR(Ipl), ColorG(Ipl), ColorB(Ipl))
If Check11.Value = 1 Then Picture1.PSet (T1 / sekG / MT - 10, Mx * (Eks - Eks00) + GR), RGB(ColorR(Ipl), ColorG(Ipl), ColorB(Ipl))
If Check15.Value = 1 Then Picture1.PSet (T1 / sekG / MT - 10, (Rsr(Ipl) - Rsr00) / MK / Mx + GR), RGB(ColorR(Ipl), ColorG(Ipl), ColorB(Ipl))
If Check16.Value = 1 Then Picture1.PSet (T1 / sekG / MT - 10, Period - Period00 + GR), RGB(ColorR(Ipl), ColorG(Ipl), ColorB(Ipl))

If kodAlfaP = 1 Then GoTo 220
AlfaPPplus = AlfaP + 180 + dAlfa: If AlfaPPplus > 360 Then AlfaPPplus = AlfaPPplus - 360 ' а это интервал углов, где будем
AlfaPPminus = AlfaP + 180 - dAlfa: If AlfaPPminus > 360 Then AlfaPPminus = AlfaPPminus - 360 'искать афелий при следующем пролете
GoTo 225
220: If NP1 > 0 Then NP1 = NP1 + 1 ' здесь перигелий
If NP1 = 0 Then NP1 = 1
Text34.Text = NP1 ' суммарное количество пролетов именно в перигелии (теперь начиная с самого первого)
AlfaPplus = AlfaP + dAlfa: If AlfaPplus > 360 Then AlfaPplus = AlfaPplus - 360 ' а это интервал углов, где будем
AlfaPminus = AlfaP - dAlfa: If AlfaPminus < 0 Then AlfaPminus = AlfaPminus + 360 'искать перигелий при следующем пролете
Call JD_DMY(JD) 'определяем дату пролета именно перигелия для времени ET (UT)
Text28.Text = IGnew: Text27.Text = IMnew: Text26.Text = IDnew

225: If Check20.Value = 1 And NP1 > 2 Then File_AlfaP 'если включена запись параметров
BettaMax = 0 'обнуляем значение угла наклона после прохождения перигелий или афелия
239: kodAlfaU = 0: GoTo 285 'если до перигелия или афелия еще не дошли находясь в интервале, т.к. kodRezult = 1
'                          а если узел находится вблизи перигелия так делать наверное нельзя
240: kodRezult = 0 'сюда пришли когда угол перигелия вышел из нужного интервала при нахождение перигелия или афелия
250: P0 = P00: RplOld = Rpl 'сюда пришли находясь в интнервале уже после нахождения перигелия или афелия по kodRezult=2

If kodAlfaU = 1 Then GoTo 285 ' если угол узла уже нашли работая в интервале углов для узла, т.е. эклиптика уже пройдена, то P0 не надо уменьшать
If AlfaUminus > 360 - 2 * dAlfa Then GoTo 260 ' если угол узла восхождения близок к 0 или 360 градусам
If AlfaGR(Ipl) > AlfaUplus Or AlfaGR(Ipl) < AlfaUminus Then GoTo 265
260: If AlfaGR(Ipl) > AlfaUplus And AlfaGR(Ipl) < AlfaUminus Then GoTo 265
P0 = P00 / Kp0: GoTo 285 ' планета находится в интервале поиска узла восхождения и P0 надо уменьшить
265: If AlfaUUminus > 360 - 2 * dAlfa Then GoTo 270 '   если угол узла восхождения близок к 0 или 360 градусам
If AlfaGR(Ipl) > AlfaUUplus Or AlfaGR(Ipl) < AlfaUUminus Then GoTo 285
270: If AlfaGR(Ipl) > AlfaUUplus And AlfaGR(Ipl) < AlfaUUminus Then GoTo 285
P0 = P00 / Kp0 ' планета находится в интервале поиска узла восхождения и P0 надо уменьшить


285: If AlfaGR(Ipl) > 359 Then P0 = P00 / Kp0 'это для периода
NC = NC + 1 'увеличиваем число общих шагов решения (не частных 4-x шагов в методе Рунге-Кутта)
T1 = T1 + P0 ' увеличиваем эфемеридное время на 1 шаг решения
'(при расчете Period и при решение дифф. уравнений время T1 всегда течет вперед)
T3 = T3 + kodTime * P0 'эфемеридное время T3 может течь и назад и используется
'для записи в File_AlfaP (для обработки данных) и File_dat (XYZ), а надо бы
'для формы 6 записать параметры в File_AlfaP по JDut от JD0 (вообще то большой
'разницы для нахождения k для функции вековых смещений Param=Param0+k*(JD-JD0) по UT
'или ET нет по этому оставлю только эфемеридное время и JDet0, т.е. TJD0, приму равным JDut
'(для расчета мгновенных значений AlfaL0 разница будет - А ГДЕ Я ИХ ТУТ ВОЗЬМУ AlfaL0)
TJD = TJD + kodTime * P0 / sekS: JD = TJD
Call JD_DMY(JD) '                                   по JD находим дату



'                           КОНЕЦ ЦИКЛА РАСЧЕТА И ВЫВОД ДАННЫХ НА ЭКРАН
If Check14.Value = 1 And IGk + 1 <= IGnew Then GoTo 350 ' эксперимент закончен по времени
If NC < NC0 Then GoTo 10 ' если не сделали NC0 шагов решения, то данные не выводим на экран
NC = 0
Text1.Text = Fix(100 * (T1 - P0) / sekG) / 100 'прошло времени в годах с начала эксперимента
Text50.Text = Fix(1000 * (Timer - Timer0) / 3600) / 1000 'прошло компьютерного времени в часах
Text44.Text = IGnew: Text43.Text = IMnew: Text42.Text = IDnew
Text14.Text = JD
Text2.Text = Fix(1000 * AlfaGR(Ipl)) / 1000 'текущее положение планеты по углу

MdAlfa = Text5.Text: MT = Text9.Text: Mx = Text21.Text: Nfile = Text8.Text
P00 = Text38.Text: NC0 = Text37.Text: Kp0 = Text39.Text: GR = Text36.Text

kodVgr = 0: If Check13.Value = 1 Then kodVgr = 1 'при kodVgr = 1 учитывать скорость гравитации



'ДЛЯ ГРАФИКОВ МОЖНО ВЫВЕСТИ И ДАННЫЕ РАССЧИТАННЫЕ ЗДЕСЬ, А НЕ В DIFFUR
If Check18.Value = 1 Then 'отладка расчета скоростей гравитации
Picture1.PSet (T1 / sekG / MT - 10, (dVj(N1, N2) - dVi(N1, N2)) / MdAlfa + 6), RGB(ColorR(N2), ColorG(N2), ColorB(N2))
Picture1.PSet (T1 / sekG / MT - 10, dR(N1, N2) / MdAlfa + 3), RGB(ColorR(N2), ColorG(N2), ColorB(N2))

Picture1.PSet (T1 / sekG / MT - 10, (dVj(N2, N1) - dVi(N2, N1)) / MdAlfa - 3), RGB(ColorR(N1), ColorG(N1), ColorB(N1))
Picture1.PSet (T1 / sekG / MT - 10, dR(N2, N1) / MdAlfa - 6), RGB(ColorR(N1), ColorG(N1), ColorB(N1))
End If




If Check2.Value = 1 Then File_dat ' если надо сохранить текущие результаты в файле

If Check3.Value = 1 Then '          если рисовать планеты
For i = 0 To Ne
Picture1.FillColor = RGB(0, 0, 0)
'                         если нужны траектории, то не стирать предыдущие точки
If Check1.Value = 0 Then Picture1.Circle ((Xold(i) - Xold(Ipl0)) / ML, (Yold(i) - Yold(Ipl0)) / ML), 0.06, RGB(0, 0, 0)
If Check1.Value = 0 Then Picture1.Circle ((Xold(i) - Xold(Ipl0)) / ML, (Zold(i) - Zold(Ipl0)) / ML + 11), 0.06, RGB(0, 0, 0)
If Check1.Value = 0 Then Picture2.Circle (Xold(i) / MLmax, Yold(i) / MLmax), 0.05, RGB(0, 0, 0)
Next i
For i = 0 To Ne
Xold(i) = X(i): Yold(i) = Y(i): Zold(i) = Z(i)
kColor = 1: If Z(i) < 0 Then kColor = 0.5
Picture1.FillColor = RGB(kColor * ColorR(i), kColor * ColorG(i), kColor * ColorB(i))
Picture1.Circle ((X(i) - X(Ipl0)) / ML, (Y(i) - Y(Ipl0)) / ML), 0.06, RGB(kColor * ColorR(i), kColor * ColorG(i), kColor * ColorB(i))
Picture1.Circle ((X(i) - X(Ipl0)) / ML, (Z(i) - Z(Ipl0)) / ML + 11), 0.06, RGB(kColor * ColorR(i), kColor * ColorG(i), kColor * ColorB(i))
Picture2.Circle ((X(i) - X(0)) / MLmax, (Y(i) - Y(0)) / MLmax), 0.05, RGB(kColor * ColorR(i), kColor * ColorG(i), kColor * ColorB(i))
Next i
End If

Text11.Text = (X(Ipl) - X(Ipl0)) / MK '                   если выводить координату X
If Check5.Value = 1 Then Picture1.PSet (T1 / sekG / MT - 10, (X(Ipl) - X(Ipl0)) / MK / Mx + GR), RGB(255, 0, 0)
Text12.Text = (Y(Ipl) - Y(Ipl0)) / MK '                   если выводить координату Y
If Check6.Value = 1 Then Picture1.PSet (T1 / sekG / MT - 10, (Y(Ipl) - Y(Ipl0)) / MK / Mx + GR), RGB(0, 0, 255)
Text30.Text = (Z(Ipl) - Z(Ipl0)) / MK '                   если выводить координату Z
If Check7.Value = 1 Then Picture1.PSet (T1 / sekG / MT - 10, (Z(Ipl) - Z(Ipl0)) / MK / Mx + GR), RGB(0, 250, 0)

'расчет ошибки решения по энергии системы и моменту количества движения при его проекции на плоскость эклиптики
If Check12.Value = 1 Then ' здесь Луна тоже относительно Солнца
    If kodModel = 1 Then ' делаем расшифровку переменных для более понятного их восприятия
For i = 0 To Ne '          от Солнца до Луны или до фиктивного объекта с №10
X(i) = Q(i): Y(i) = Q(i + Ne + 1): Z(i) = Q(i + 2 * (Ne + 1))
VX(i) = Q(i + 3 * (Ne + 1)): VY(i) = Q(i + 4 * (Ne + 1)): VZ(i) = Q(i + 5 * (Ne + 1))
Next i
    End If
V(0) = Sqr(VX(0) * VX(0) + VY(0) * VY(0) + VZ(0) * VZ(0))
MKDsys = 0 '   это возможно только относительно Солнца и только в плоскости XY
Esys = kodIpl(0) * m(0) * V(0) * V(0) / 2 'энергия системы вместе с Солнцем и когда
Ipl2 = 9: If kodModel = 1 Then Ipl2 = 10 '     в имитаторах не может быть 10-го объекта
For i = 1 To Ipl2 '  за нулевой уровень для потенциальной энергии планет принято R0(I)
DX(0, i) = X(i) - X(0): DY(0, i) = Y(i) - Y(0): DZ(0, i) = Z(i) - Z(0)
R(0, i) = Sqr(DX(0, i) * DX(0, i) + DY(0, i) * DY(0, i) + DZ(0, i) * DZ(0, i))
V(i) = Sqr(VX(i) * VX(i) + VY(i) * VY(i) + VZ(i) * VZ(i))
MKDsys = MKDsys + kodIpl(i) * m(i) * ((VY(i) - VY(0)) * DX(0, i) - (VX(i) - VX(0)) * DY(0, i))
Esys = Esys + kodIpl(i) * m(i) * V(i) * V(i) / 2 + kodIpl(i) * gamma * m(0) * m(i) * (1 / R0(i) - 1 / R(0, i))
Next i
'                масштаб энергии соответствует масштабу MX (процентов/см)
Picture1.Circle (T1 / sekG / MT - 10, 100 * (Esys - E0sys) / E0sys / Mx), 0.02, RGB(255, 0, 0)
Picture1.Circle (T1 / sekG / MT - 10, 100 * (MKDsys - MKD0sys) / MKD0sys / Mx), 0.02, RGB(0, 0, 255)
Text19.Text = Fix(10000 * MK * (Esys - E0sys) / E0sys) / MK / 100
Text41.Text = Fix(10000 * MK * (MKDsys - MKD0sys) / MKD0sys) / MK / 100
End If ' P.S.   можно и для одной планеты посчитать, т.к. используется kodIpl(i)

300: If kodEnd = 1 Then GoTo 350
Do While DoEvents() 'идем на 10, после обработки сообщений в очереди
GoTo 10 '            операционной системы, а если будет сообщение от кнопки <Остановить>,
Loop '               то после его обработки, kodEnd будет равен 1 и будет GoTo 350
GoTo 350
330: MsgBox "Ошибка решения уравнений", , "Solsys7": Exit Sub
350: If NPlan > 0 Then Form15.Text3.Text = Timer 'по плану, если NPlan > 0
NPlan = 0
End Sub

Private Sub Text53_Change() 'выполняем эксперименты по плану с формы 15
Command1_Click
End Sub





Private Sub Diffur21() 'классические дифференциальные уравнения для второй формы, т.е. для всей
For i = 0 To Ne '    Солнечной системы при движение всех планет по трем координатам и с учетом скорости гравитации
X(i) = Q(i): Y(i) = Q(i + Ne + 1): Z(i) = Q(i + 2 * (Ne + 1)) ' для удобства работы производим замену обозначений
VX(i) = Q(i + 3 * (Ne + 1)): VY(i) = Q(i + 4 * (Ne + 1)): VZ(i) = Q(i + 5 * (Ne + 1))
Next i

For i = 0 To Ne '  NE=10'           первое приближение - без запаздывания по координатам
For j = 0 To Ne
If i = j Then GoTo 200
DX(i, j) = X(j) - X(i)
DY(i, j) = Y(j) - Y(i)
DZ(i, j) = Z(j) - Z(i)
Rxy(i, j) = Sqr(DX(i, j) * DX(i, j) + DY(i, j) * DY(i, j)) '  находим расстояния между объектами
R(i, j) = Sqr(Rxy(i, j) * Rxy(i, j) + DZ(i, j) * DZ(i, j))
T(i, j) = R(i, j) / Vgr '            время необходимое для распространения гравитации
F(i, j) = kodF(i, j) * m(i) * m(j) * gamma / R(i, j) ^ PStR 'находим силу притяжения между следом 1-го тела воздействующим
200: Next j
Next i
If Check21.Value = 1 Then InfoM(5) = "F12=   " & F(N1, N2) & "    R12=   " & R(N1, N2) & "    T12=   " & T(N1, N2) 'сохраняем для теста

If kodVgr = 0 Then GoTo 310 'скорость равна бесконечности и запаздывания по координатам не будет
If kodDX = 0 Then GoTo 310 ' запаздывания по координатам не учитывать

For k = 1 To 2 '               количество итераций равно 2
For i = 0 To Ne '  NE=10
For j = 0 To Ne
If i = j Then GoTo 300 '                  расчет координат следов планет
DX(i, j) = X(j) - X(i) + VX(i) * T(i, j) 'расстояние между следом 1-го тела воздействующим на 2-е тело по оси X
DY(i, j) = Y(j) - Y(i) + VY(i) * T(i, j)
DZ(i, j) = Z(j) - Z(i) + VZ(i) * T(i, j)
Rxy(i, j) = Sqr(DX(i, j) * DX(i, j) + DY(i, j) * DY(i, j)) ' расстояния в плоскости XY между следом 1-го тела воздействующим на 2-е тело
R(i, j) = Sqr(Rxy(i, j) * Rxy(i, j) + DZ(i, j) * DZ(i, j)) ' общее расстояния между следом 1-го тела воздействующим на 2-е тело
T(i, j) = R(i, j) / Vgr '                                    время необходимое для распространения гравитации
300: Next j
Next i
InfoM(4) = " число итераций DX=  " & k
Next k

310:
For i = 0 To Ne '  NE=10                                      расчет сил
For j = 0 To Ne '       и если задана скорость гравитации и kodDX = 0, то и с учетом запаздывания по координатам
If i = j Then GoTo 400
F(i, j) = kodF(i, j) * m(i) * m(j) * gamma / R(i, j) ^ PStR 'находим силу притяжения между следом 1-го тела воздействующим
Betta(i, j) = Atn(DZ(i, j) / Rxy(i, j)) '                    на 2-е тело с учетом наличия i и j планет, т.е. kodF(i, j)
Alfa12(i, j) = pi / 2: If DY(i, j) < 0 Then Alfa12(i, j) = 3 * pi / 2 'угол пригодиться при расчете сил давления
If DX(i, j) = 0 Then GoTo 390
Alfa12(i, j) = Atn(DY(i, j) / DX(i, j))
390:
FX(i, j) = F(i, j) * Abs(Cos(Alfa12(i, j)) * Cos(Betta(i, j))) 'находим проекции силы притяжения на оси X и Y
If DX(i, j) > 0 Then FX(i, j) = -FX(i, j)
FY(i, j) = F(i, j) * Abs(Sin(Alfa12(i, j)) * Cos(Betta(i, j))) 'сила действующая со стороны i-того тела на j-тое тело
If DY(i, j) > 0 Then FY(i, j) = -FY(i, j)
FZ(i, j) = F(i, j) * Abs(Sin(Betta(i, j))) 'при расчете FX(i, j) и FY(i, j) можно без Abs, т.к. там косинус
If DZ(i, j) > 0 Then FZ(i, j) = -FZ(i, j)
400: Next j
Next i

If Check21.Value = 1 Then ' если вывести на экран расчет теста для давления гравитации
'InfoM(1) = "DX12=    " & Fix(DX(1, 2) / MK) & "     DY12=    " & Fix(DY(1, 2) / MK) & "     DZ12=    " & Fix(DZ(1, 2) / MK)
InfoM(1) = "X1-dX=  " & X(N1) - VX(N1) * T(N1, N2) & "     Y1-dY=  " & Y(N1) - VY(N1) * T(N1, N2) & "     Z1-dZ=  " & Z(N1) - VZ(N1) * T(N1, N2)
InfoM(2) = "FX12= " & Fix(FX(N1, N2)) & "     FY= " & Fix(FY(N1, N2)) & "     FZ= " & Fix(FZ(N1, N2))
InfoM(3) = "F12=   " & F(N1, N2) & "    R12=   " & R(N1, N2) & "    T12=   " & T(N1, N2)
'Info1 = "DX21=    " & Fix(DX(2, 1) / MK) & "     DY21=    " & Fix(DY(2, 1) / MK) & "     DZ21=    " & Fix(DZ(2, 1) / MK)
Info1 = "X2-dX=  " & X(N2) - VX(N2) * T(N2, N1) & "     Y2-dY=  " & Y(N2) - VY(N2) * T(N2, N1) & "     Z2-dZ=  " & Z(N2) - VZ(N2) * T(N2, N1)
Info2 = "FX21= " & Fix(FX(N2, N1)) & "     FY= " & Fix(FY(N2, N1)) & "     FZ= " & Fix(FZ(N2, N1))
Info3 = "F21=   " & F(N2, N1) & "    R21=   " & R(N2, N1) & "    T21=   " & T(N2, N1)
For i = 0 To 10
X1(i) = X(i): X2(i) = Y(i): X3(i) = Z(i)
X4(i) = VX(i): X5(i) = VY(i): X6(i) = VZ(i)
Next i
End If

If kodPgr = 0 Then GoTo 500 ' если не учитывать давление гравитации

NVgr = NVgr + 1 'счетчик обращений к Diffur (при разных обращениях по разному считается dR)

If kod4diff = 0 Then GoTo 492 '            если не выполнять уточненый расчет dVi(i, j)  dVj(i, j)

For i = 0 To Ne
BettaVZ(i) = pi / 2: If VZ(i) < 0 Then BettaVZ(i) = -pi / 2 ' находим углы скоростей i-х и j-х тел
AlfaVX(i) = 0: If VX(i) < 0 Then AlfaVX(i) = pi
AlfaVY(i) = pi / 2: If VY(i) < 0 Then AlfaVY(i) = 3 * pi / 2
Next i

For i = 0 To Ne
For j = 0 To Ne '                                     Alfa12(i, j) и Betta(i, j) это направление скорости гравитации
If i = j Then GoTo 490 '                              при ее распространении от i-того тела до j-того тела
'DX(i, j) = X(j) - X(i) - VX(i) * T(i, j) 'расстояние между следом 1-го тела (i) воздействующим на 2-е тело (j) по оси
'X для расчета НАПРАВЛЕНИЯ СКОРОСТИ ГРАВИТАЦИИ а для расчета силы действующей со стороны i-того тела на j-тое тело
'направление будет обратное по этому выше при расчете сил было получено Alfa12(i, j) = Atn(DY(i, j) / DX(i, j))
'а использовано было Abs(Sin(Alfa12(i,j))) и Abs(Cos(Alfa12(i,j))) а потом If DX(i,j) > 0 Then FX(i,j)= -FX(i,j)
'                        А ТЕПЕРЬ ПЕРЕВОДИМ Alfa12(i,j) ПОЛУЧЕННУЮ ВЫШЕ В РАЗМЕРНОСТЬ ОТ 0 ДО 2*pi
If DX(i, j) > 0 Then GoTo 420 '                     переход в 1-ю и 4-ю четверти
Alfa12(i, j) = Alfa12(i, j) + pi: GoTo 480 '        2-я и 3-я четверти (во 2-ой четверти Alfa12(i,j) был с минусом)
420:  If DY(i, j) > 0 Then GoTo 480 '               значит была 1-я четверть
Alfa12(i, j) = Alfa12(i, j) + 2 * pi '              это 4-я четверть (Alfa12(i,j) был с минусом)
480:

dAlfaVXi(i, j) = AlfaVX(i) - Alfa12(i, j) '      dAlfaVXi(i, j)- угол между скоростью по оси X i-го тела и проекцией вектора
If dAlfaVXi(i, j) > pi Then dAlfaVXi(i, j) = dAlfaVXi(i, j) - 2 * pi 'скорости гравитации от i-го тела к j-му на плоскость XY
If dAlfaVXi(i, j) < -pi Then dAlfaVXi(i, j) = dAlfaVXi(i, j) + 2 * pi '
dAlfaVYi(i, j) = AlfaVY(i) - Alfa12(i, j)
If dAlfaVYi(i, j) > pi Then dAlfaVYi(i, j) = dAlfaVYi(i, j) - 2 * pi
If dAlfaVYi(i, j) < -pi Then dAlfaVYi(i, j) = dAlfaVYi(i, j) + 2 * pi
dBettaVZi(i, j) = BettaVZ(i) - Betta(i, j)

dAlfaVXj(i, j) = AlfaVX(j) - Alfa12(i, j) '      dAlfaVX(i, j)- угол между скоростью по оси X j-го тела и проекцией вектора
If dAlfaVXj(i, j) > pi Then dAlfaVXj(i, j) = dAlfaVXj(i, j) - 2 * pi 'скорости гравитации от i-го тела к j-му на плоскость XY
If dAlfaVXj(i, j) < -pi Then dAlfaVXj(i, j) = dAlfaVXj(i, j) + 2 * pi '
dAlfaVYj(i, j) = AlfaVY(j) - Alfa12(i, j)
If dAlfaVYj(i, j) > pi Then dAlfaVYj(i, j) = dAlfaVYj(i, j) - 2 * pi
If dAlfaVYj(i, j) < -pi Then dAlfaVYj(i, j) = dAlfaVYj(i, j) + 2 * pi
dBettaVZj(i, j) = BettaVZ(j) - Betta(i, j)

'                                                    СТАРЫЙ ВАРИАНТ
'If kodDX = 1 Then dVi(i, j) = dVi(i, j) - d2Vi(i, j) * T(i, j) ' а если было запаздывание по координатам, то надо
'найти dVi(i, j) в следе планеты, а это уточнение будет только если в param будет расчитано d2Vi(i, j)
VXold(i, j) = VX(i): VYold(i, j) = VY(i): VZold(i, j) = VZ(i) ' НАХОЖУ СКОРОСТИ В СЛЕДАХ
If kodVgr = 1 And kodDX = 1 Then  ' запаздывания по координатам учитывать
VXold(i, j) = VX(i) - dVX(i) * T(i, j)
VYold(i, j) = VY(i) - dVY(i) * T(i, j)
VZold(i, j) = VZ(i) - dVZ(i) * T(i, j)
End If

dVZi(i, j) = Cos(dBettaVZi(i, j)) * Abs(VZold(i, j)) '                        проекции скоростей i-го тела по осям координат
dVXi(i, j) = Abs(Cos(Betta(i, j))) * Cos(dAlfaVXi(i, j)) * Abs(VXold(i, j)) ' на вектор скорости гравитации от i-го тела
dVYi(i, j) = Abs(Cos(Betta(i, j))) * Cos(dAlfaVYi(i, j)) * Abs(VYold(i, j))
dVZj(i, j) = Cos(dBettaVZj(i, j)) * Abs(VZ(j)) '                        проекции скоростей j-го тела по осям координат
dVXj(i, j) = Abs(Cos(Betta(i, j))) * Cos(dAlfaVXj(i, j)) * Abs(VX(j)) ' на вектор скорости гравитации от i-го тела
dVYj(i, j) = Abs(Cos(Betta(i, j))) * Cos(dAlfaVYj(i, j)) * Abs(VY(j))
dVi(i, j) = dVZi(i, j) + dVXi(i, j) + dVYi(i, j)
dVj(i, j) = dVZj(i, j) + dVXj(i, j) + dVYj(i, j)
dR(i, j) = dVj(i, j) - dVi(i, j)
490: Next j
Next i

492:
For i = 0 To Ne
For j = 0 To Ne '
If i = j Then GoTo 495
If kodPgr = 2 Then F(i, j) = F(i, j) * (1 - kPgr1 * dR(i, j) * dR(i, j) / Vgr / Vgr + kPgr2 * R(i, j) * d2R(i, j) / Vgr / Vgr)
If kodPgr = 3 Then F(i, j) = F(i, j) * (1 - kPgr1 * dR(i, j) * dR(i, j) / (Vgr + dVi(i, j)) / (Vgr + dVi(i, j)) + kPgr2 * R(i, j) * d2R(i, j) / (Vgr + dVi(i, j)) / (Vgr + dVi(i, j)))

If kodPgr = 5 Then F(i, j) = F(i, j) * (1 - kPgr1 * dVj(i, j) / Vgr) / (1 - kPgr2 * dVi(i, j) / Vgr)
If kodPgr = 6 Then F(i, j) = F(i, j) * (1 - kPgr1 * dR(i, j) / Vgr)
If kodPgr = 7 Then '                                                 Сухоруков
Vpl = Sqr(VX(j) * VX(j) + VY(j) * VY(j) + VZ(j) * VZ(j)) 'скорость j
F(i, j) = F(i, j) * Sqr(Vgr * Vgr + Vpl * Vpl - 2 * Vgr * dVj(i, j)) '/ Vgr
Vpl = Sqr(VX(i) * VX(i) + VY(i) * VY(i) + VZ(i) * VZ(i)) 'скорость i
F(i, j) = F(i, j) / Sqr(Vgr * Vgr + Vpl * Vpl - 2 * Vgr * dVi(i, j)) ' / Vgr
End If

FZ(i, j) = Sgn(FZ(i, j)) * F(i, j) * Abs(Sin(Betta(i, j)))
FX(i, j) = Sgn(FX(i, j)) * F(i, j) * Abs(Cos(Alfa12(i, j)) * Cos(Betta(i, j)))
FY(i, j) = Sgn(FY(i, j)) * F(i, j) * Abs(Sin(Alfa12(i, j)) * Cos(Betta(i, j)))
495: Next j
Next i

500: 'for i=0 to 32 :w(i)=q(33+i):next i
W(0) = Q(33): W(1) = Q(34): W(2) = Q(35): W(3) = Q(36): W(4) = Q(37): W(5) = Q(38) ' скорости по оси Х
W(6) = Q(39): W(7) = Q(40): W(8) = Q(41): W(9) = Q(42): W(10) = Q(43)

W(11) = Q(44): W(12) = Q(45): W(13) = Q(46): W(14) = Q(47): W(15) = Q(48): W(16) = Q(49) ' скорости по оси Y
W(17) = Q(50): W(18) = Q(51): W(19) = Q(52): W(20) = Q(53): W(21) = Q(54)

W(22) = Q(55): W(23) = Q(56): W(24) = Q(57): W(25) = Q(58): W(26) = Q(59): W(27) = Q(60) ' скорости по оси Z
W(28) = Q(61): W(29) = Q(62): W(30) = Q(63): W(31) = Q(64): W(32) = Q(65)

W(33) = (FX(1, 0) + FX(2, 0) + FX(3, 0) + FX(4, 0) + FX(5, 0) + FX(6, 0) + FX(7, 0) + FX(8, 0) + FX(9, 0) + FX(10, 0)) / m(0)  ' ускорение Солнца по оси Х
W(34) = (FX(0, 1) + FX(2, 1) + FX(3, 1) + FX(4, 1) + FX(5, 1) + FX(6, 1) + FX(7, 1) + FX(8, 1) + FX(9, 1) + FX(10, 1)) / m(1) ' ускорение Меркурия по оси Х
W(35) = (FX(0, 2) + FX(1, 2) + FX(3, 2) + FX(4, 2) + FX(5, 2) + FX(6, 2) + FX(7, 2) + FX(8, 2) + FX(9, 2) + FX(10, 2)) / m(2)
W(36) = (FX(0, 3) + FX(1, 3) + FX(2, 3) + FX(4, 3) + FX(5, 3) + FX(6, 3) + FX(7, 3) + FX(8, 3) + FX(9, 3) + FX(10, 3)) / m(3)
W(37) = (FX(0, 4) + FX(1, 4) + FX(2, 4) + FX(3, 4) + FX(5, 4) + FX(6, 4) + FX(7, 4) + FX(8, 4) + FX(9, 4) + FX(10, 4)) / m(4)
W(38) = (FX(0, 5) + FX(1, 5) + FX(2, 5) + FX(3, 5) + FX(4, 5) + FX(6, 5) + FX(7, 5) + FX(8, 5) + FX(9, 5) + FX(10, 5)) / m(5)
W(39) = (FX(0, 6) + FX(1, 6) + FX(2, 6) + FX(3, 6) + FX(4, 6) + FX(5, 6) + FX(7, 6) + FX(8, 6) + FX(9, 6) + FX(10, 6)) / m(6)
W(40) = (FX(0, 7) + FX(1, 7) + FX(2, 7) + FX(3, 7) + FX(4, 7) + FX(5, 7) + FX(6, 7) + FX(8, 7) + FX(9, 7) + FX(10, 7)) / m(7)
W(41) = (FX(0, 8) + FX(1, 8) + FX(2, 8) + FX(3, 8) + FX(4, 8) + FX(5, 8) + FX(6, 8) + FX(7, 8) + FX(9, 8) + FX(10, 8)) / m(8)
W(42) = (FX(0, 9) + FX(1, 9) + FX(2, 9) + FX(3, 9) + FX(4, 9) + FX(5, 9) + FX(6, 9) + FX(7, 9) + FX(8, 9) + FX(10, 9)) / m(9)
W(43) = (FX(0, 10) + FX(1, 10) + FX(2, 10) + FX(3, 10) + FX(4, 10) + FX(5, 10) + FX(6, 10) + FX(7, 10) + FX(8, 10) + FX(9, 10)) / m(10)

W(44) = (FY(1, 0) + FY(2, 0) + FY(3, 0) + FY(4, 0) + FY(5, 0) + FY(6, 0) + FY(7, 0) + FY(8, 0) + FY(9, 0) + FY(10, 0)) / m(0)  ' ускорение Солнца по оси Y
W(45) = (FY(0, 1) + FY(2, 1) + FY(3, 1) + FY(4, 1) + FY(5, 1) + FY(6, 1) + FY(7, 1) + FY(8, 1) + FY(9, 1) + FY(10, 1)) / m(1) ' ускорение Меркурия по оси Y
W(46) = (FY(0, 2) + FY(1, 2) + FY(3, 2) + FY(4, 2) + FY(5, 2) + FY(6, 2) + FY(7, 2) + FY(8, 2) + FY(9, 2) + FY(10, 2)) / m(2)
W(47) = (FY(0, 3) + FY(1, 3) + FY(2, 3) + FY(4, 3) + FY(5, 3) + FY(6, 3) + FY(7, 3) + FY(8, 3) + FY(9, 3) + FY(10, 3)) / m(3)
W(48) = (FY(0, 4) + FY(1, 4) + FY(2, 4) + FY(3, 4) + FY(5, 4) + FY(6, 4) + FY(7, 4) + FY(8, 4) + FY(9, 4) + FY(10, 4)) / m(4)
W(49) = (FY(0, 5) + FY(1, 5) + FY(2, 5) + FY(3, 5) + FY(4, 5) + FY(6, 5) + FY(7, 5) + FY(8, 5) + FY(9, 5) + FY(10, 5)) / m(5)
W(50) = (FY(0, 6) + FY(1, 6) + FY(2, 6) + FY(3, 6) + FY(4, 6) + FY(5, 6) + FY(7, 6) + FY(8, 6) + FY(9, 6) + FY(10, 6)) / m(6)
W(51) = (FY(0, 7) + FY(1, 7) + FY(2, 7) + FY(3, 7) + FY(4, 7) + FY(5, 7) + FY(6, 7) + FY(8, 7) + FY(9, 7) + FY(10, 7)) / m(7)
W(52) = (FY(0, 8) + FY(1, 8) + FY(2, 8) + FY(3, 8) + FY(4, 8) + FY(5, 8) + FY(6, 8) + FY(7, 8) + FY(9, 8) + FY(10, 8)) / m(8)
W(53) = (FY(0, 9) + FY(1, 9) + FY(2, 9) + FY(3, 9) + FY(4, 9) + FY(5, 9) + FY(6, 9) + FY(7, 9) + FY(8, 9) + FY(10, 9)) / m(9)
W(54) = (FY(0, 10) + FY(1, 10) + FY(2, 10) + FY(3, 10) + FY(4, 10) + FY(5, 10) + FY(6, 10) + FY(7, 10) + FY(8, 10) + FY(9, 10)) / m(10)

W(55) = (FZ(1, 0) + FZ(2, 0) + FZ(3, 0) + FZ(4, 0) + FZ(5, 0) + FZ(6, 0) + FZ(7, 0) + FZ(8, 0) + FZ(9, 0) + FZ(10, 0)) / m(0) ' ускорение Солнца по оси Z
W(56) = (FZ(0, 1) + FZ(2, 1) + FZ(3, 1) + FZ(4, 1) + FZ(5, 1) + FZ(6, 1) + FZ(7, 1) + FZ(8, 1) + FZ(9, 1) + FZ(10, 1)) / m(1) ' ускорение Меркурия по оси Z
W(57) = (FZ(0, 2) + FZ(1, 2) + FZ(3, 2) + FZ(4, 2) + FZ(5, 2) + FZ(6, 2) + FZ(7, 2) + FZ(8, 2) + FZ(9, 2) + FZ(10, 2)) / m(2)
W(58) = (FZ(0, 3) + FZ(1, 3) + FZ(2, 3) + FZ(4, 3) + FZ(5, 3) + FZ(6, 3) + FZ(7, 3) + FZ(8, 3) + FZ(9, 3) + FZ(10, 3)) / m(3)
W(59) = (FZ(0, 4) + FZ(1, 4) + FZ(2, 4) + FZ(3, 4) + FZ(5, 4) + FZ(6, 4) + FZ(7, 4) + FZ(8, 4) + FZ(9, 4) + FZ(10, 4)) / m(4)
W(60) = (FZ(0, 5) + FZ(1, 5) + FZ(2, 5) + FZ(3, 5) + FZ(4, 5) + FZ(6, 5) + FZ(7, 5) + FZ(8, 5) + FZ(9, 5) + FZ(10, 5)) / m(5)
W(61) = (FZ(0, 6) + FZ(1, 6) + FZ(2, 6) + FZ(3, 6) + FZ(4, 6) + FZ(5, 6) + FZ(7, 6) + FZ(8, 6) + FZ(9, 6) + FZ(10, 6)) / m(6)
W(62) = (FZ(0, 7) + FZ(1, 7) + FZ(2, 7) + FZ(3, 7) + FZ(4, 7) + FZ(5, 7) + FZ(6, 7) + FZ(8, 7) + FZ(9, 7) + FZ(10, 7)) / m(7)
W(63) = (FZ(0, 8) + FZ(1, 8) + FZ(2, 8) + FZ(3, 8) + FZ(4, 8) + FZ(5, 8) + FZ(6, 8) + FZ(7, 8) + FZ(9, 8) + FZ(10, 8)) / m(8)
W(64) = (FZ(0, 9) + FZ(1, 9) + FZ(2, 9) + FZ(3, 9) + FZ(4, 9) + FZ(5, 9) + FZ(6, 9) + FZ(7, 9) + FZ(8, 9) + FZ(10, 9)) / m(9)
W(65) = (FZ(0, 10) + FZ(1, 10) + FZ(2, 10) + FZ(3, 10) + FZ(4, 10) + FZ(5, 10) + FZ(6, 10) + FZ(7, 10) + FZ(8, 10) + FZ(9, 10)) / m(10)

If Check23.Value = 1 Then '              если вывести текст   dR   dVj-dVi
InfoM(6 + NVgr) = "NVgr=   " & NVgr & "    dR=      " & dR(N1, N2) & "    d2R=   " & d2R(N1, N2)
InfoM(16 + NVgr) = "NVgr=   " & NVgr & "  dVj-dVi= " & dVj(N1, N2) - dVi(N1, N2) & "    d2Vij=   " & d2Vij(N1, N2)
End If

End Sub




'Private Sub Diffur22() 'ОТОшные дифференциальные уравнения для второй формы, т.е. для всей Солнечной системы
'For I% = 0 To NE '    Солнечной системы при движение всех планет по трем координатам и с учетом скорости гравитации
'X(I%) = Q(I%): Y(I%) = Q(I% + NE + 1): Z(I) = Q(I + 2 * (NE + 1)) ' для удобства работы производим замену обозначений
'VX(I%) = Q(I% + 3 * (NE + 1)): VY(I%) = Q(I% + 4 * (NE + 1)): VZ(I) = Q(I + 5 * (NE + 1))
'Next I%

'For I = 0 To NE '  NE=10
'For J = 0 To NE
'If I% = J% Then GoTo 200 '           первое приближение в покоящейся системе
'DX(I, J) = X(J) - X(I)
'DY(I, J) = Y(J) - Y(I)
'DZ(I, J) = Z(J) - Z(I)
'Rxy(I, J) = Sqr(DX(I, J) * DX(I, J) + DY(I, J) * DY(I, J)) '  находим расстояния между объектами
'R(I, J) = Sqr(Rxy(I, J) * Rxy(I, J) + DZ(I, J) * DZ(I, J))
'T(I, J) = R(I, J) / Vgr '            время необходимое для распространения гравитации
'200: Next J
'Next I
'If kodVgr = 0 Then GoTo 310

'For K = 0 To NmaxVgr '               количество итераций равно NmaxVgr
'For I = 0 To NE '  NE=10
'For J = 0 To NE
'If I% = J% Then GoTo 300 '           расчет координат следов планет
'DX(I, J) = X(J) - X(I) + VX(I) * T(I, J) 'расстояние между следом 1-го тела воздействующим на 2-е тело по оси X
'DY(I, J) = Y(J) - Y(I) + VY(I) * T(I, J)
'DZ(I, J) = Z(J) - Z(I) + VZ(I) * T(I, J)
'Rxy(I, J) = Sqr(DX(I, J) * DX(I, J) + DY(I, J) * DY(I, J)) ' расстояния в плоскости XY между следом 1-го тела воздействующим на 2-е тело
'R(I, J) = Sqr(Rxy(I, J) * Rxy(I, J) + DZ(I, J) * DZ(I, J)) ' общее расстояния между следом 1-го тела воздействующим на 2-е тело
'T(I, J) = R(I, J) / Vgr '                                    время необходимое для распространения гравитации
'300: Next J
'Next I
'Next K

'310: For I = 0 To NE '  NE=10
'For J = 0 To NE
'If I% = J% Then GoTo 400 '                                      расчет сил
'If DX(I, J) = 0 Then GoTo 400
'F(I, J) = kodF(I, J) * m(I) * m(J) * gamma / R(I, J) / R(I, J) 'находим силу притяжения между следом 1-го тела воздействующим на 2-е тело
'Betta(I, J) = Atn(Abs(DZ(I, J)) / Rxy(I, J))
'Alfa12(I, J) = Abs(Atn(DY(I, J) / DX(I, J)))
'FX(I, J) = F(I, J) * Cos(Alfa12(I, J)) * Cos(Betta(I, J)) 'находим проекции силы притяжения на оси X и Y
'If DX(I, J) > 0 Then FX(I, J) = -FX(I, J)
'FY(I, J) = F(I, J) * Sin(Alfa12(I, J)) * Cos(Betta(I, J)) 'сила действующая со стороны Jтого тела на Iтое тело
'If DY(I, J) > 0 Then FY(I, J) = -FY(I, J)
'FZ(I, J) = F(I, J) * Sin(Betta(I, J))
'If DZ(I, J) > 0 Then FZ(I, J) = -FZ(I, J)
'400: Next J
'Next I

'500:
'W(0) = Q(33): W(1) = Q(34): W(2) = Q(35): W(3) = Q(36): W(4) = Q(37): W(5) = Q(38) ' скорости по оси Х
'W(6) = Q(39): W(7) = Q(40): W(8) = Q(41): W(9) = Q(42): W(10) = Q(43)

'W(11) = Q(44): W(12) = Q(45): W(13) = Q(46): W(14) = Q(47): W(15) = Q(48): W(16) = Q(49) ' скорости по оси Y
'W(17) = Q(50): W(18) = Q(51): W(19) = Q(52): W(20) = Q(53): W(21) = Q(54)

'W(22) = Q(55): W(23) = Q(56): W(24) = Q(57): W(25) = Q(58): W(26) = Q(59): W(27) = Q(60) ' скорости по оси Y
'W(28) = Q(61): W(29) = Q(62): W(30) = Q(63): W(31) = Q(64): W(32) = Q(65)

'W(33) = (FX(1, 0) + FX(2, 0) + FX(3, 0) + FX(4, 0) + FX(5, 0) + FX(6, 0) + FX(7, 0) + FX(8, 0) + FX(9, 0) + FX(10, 0)) / m(0)  ' ускорение Солнца по оси Х
'W(34) = (FX(0, 1) + FX(2, 1) + FX(3, 1) + FX(4, 1) + FX(5, 1) + FX(6, 1) + FX(7, 1) + FX(8, 1) + FX(9, 1) + FX(10, 1)) / m(1) ' ускорение Меркурия по оси Х
'W(35) = (FX(0, 2) + FX(1, 2) + FX(3, 2) + FX(4, 2) + FX(5, 2) + FX(6, 2) + FX(7, 2) + FX(8, 2) + FX(9, 2) + FX(10, 2)) / m(2)
'W(36) = (FX(0, 3) + FX(1, 3) + FX(2, 3) + FX(4, 3) + FX(5, 3) + FX(6, 3) + FX(7, 3) + FX(8, 3) + FX(9, 3) + FX(10, 3)) / m(3)
'W(37) = (FX(0, 4) + FX(1, 4) + FX(2, 4) + FX(3, 4) + FX(5, 4) + FX(6, 4) + FX(7, 4) + FX(8, 4) + FX(9, 4) + FX(10, 4)) / m(4)
'W(38) = (FX(0, 5) + FX(1, 5) + FX(2, 5) + FX(3, 5) + FX(4, 5) + FX(6, 5) + FX(7, 5) + FX(8, 5) + FX(9, 5) + FX(10, 5)) / m(5)
'W(39) = (FX(0, 6) + FX(1, 6) + FX(2, 6) + FX(3, 6) + FX(4, 6) + FX(5, 6) + FX(7, 6) + FX(8, 6) + FX(9, 6) + FX(10, 6)) / m(6)
'W(40) = (FX(0, 7) + FX(1, 7) + FX(2, 7) + FX(3, 7) + FX(4, 7) + FX(5, 7) + FX(6, 7) + FX(8, 7) + FX(9, 7) + FX(10, 7)) / m(7)
'W(41) = (FX(0, 8) + FX(1, 8) + FX(2, 8) + FX(3, 8) + FX(4, 8) + FX(5, 8) + FX(6, 8) + FX(7, 8) + FX(9, 8) + FX(10, 8)) / m(8)
'W(42) = (FX(0, 9) + FX(1, 9) + FX(2, 9) + FX(3, 9) + FX(4, 9) + FX(5, 9) + FX(6, 9) + FX(7, 9) + FX(8, 9) + FX(10, 9)) / m(9)
'W(43) = (FX(0, 10) + FX(1, 10) + FX(2, 10) + FX(3, 10) + FX(4, 10) + FX(5, 10) + FX(6, 10) + FX(7, 10) + FX(8, 10) + FX(9, 10)) / m(10)

'W(44) = (FY(1, 0) + FY(2, 0) + FY(3, 0) + FY(4, 0) + FY(5, 0) + FY(6, 0) + FY(7, 0) + FY(8, 0) + FY(9, 0) + FY(10, 0)) / m(0)  ' ускорение Солнца по оси Y
'W(45) = (FY(0, 1) + FY(2, 1) + FY(3, 1) + FY(4, 1) + FY(5, 1) + FY(6, 1) + FY(7, 1) + FY(8, 1) + FY(9, 1) + FY(10, 1)) / m(1) ' ускорение Меркурия по оси Y
'W(46) = (FY(0, 2) + FY(1, 2) + FY(3, 2) + FY(4, 2) + FY(5, 2) + FY(6, 2) + FY(7, 2) + FY(8, 2) + FY(9, 2) + FY(10, 2)) / m(2)
'W(47) = (FY(0, 3) + FY(1, 3) + FY(2, 3) + FY(4, 3) + FY(5, 3) + FY(6, 3) + FY(7, 3) + FY(8, 3) + FY(9, 3) + FY(10, 3)) / m(3)
'W(48) = (FY(0, 4) + FY(1, 4) + FY(2, 4) + FY(3, 4) + FY(5, 4) + FY(6, 4) + FY(7, 4) + FY(8, 4) + FY(9, 4) + FY(10, 4)) / m(4)
'W(49) = (FY(0, 5) + FY(1, 5) + FY(2, 5) + FY(3, 5) + FY(4, 5) + FY(6, 5) + FY(7, 5) + FY(8, 5) + FY(9, 5) + FY(10, 5)) / m(5)
'W(50) = (FY(0, 6) + FY(1, 6) + FY(2, 6) + FY(3, 6) + FY(4, 6) + FY(5, 6) + FY(7, 6) + FY(8, 6) + FY(9, 6) + FY(10, 6)) / m(6)
'W(51) = (FY(0, 7) + FY(1, 7) + FY(2, 7) + FY(3, 7) + FY(4, 7) + FY(5, 7) + FY(6, 7) + FY(8, 7) + FY(9, 7) + FY(10, 7)) / m(7)
'W(52) = (FY(0, 8) + FY(1, 8) + FY(2, 8) + FY(3, 8) + FY(4, 8) + FY(5, 8) + FY(6, 8) + FY(7, 8) + FY(9, 8) + FY(10, 8)) / m(8)
'W(53) = (FY(0, 9) + FY(1, 9) + FY(2, 9) + FY(3, 9) + FY(4, 9) + FY(5, 9) + FY(6, 9) + FY(7, 9) + FY(8, 9) + FY(10, 9)) / m(9)
'W(54) = (FY(0, 10) + FY(1, 10) + FY(2, 10) + FY(3, 10) + FY(4, 10) + FY(5, 10) + FY(6, 10) + FY(7, 10) + FY(8, 10) + FY(9, 10)) / m(10)

'W(55) = (FZ(1, 0) + FZ(2, 0) + FZ(3, 0) + FZ(4, 0) + FZ(5, 0) + FZ(6, 0) + FZ(7, 0) + FZ(8, 0) + FZ(9, 0) + FZ(10, 0)) / m(0) ' ускорение Солнца по оси Z
'W(56) = (FZ(0, 1) + FZ(2, 1) + FZ(3, 1) + FZ(4, 1) + FZ(5, 1) + FZ(6, 1) + FZ(7, 1) + FZ(8, 1) + FZ(9, 1) + FZ(10, 1)) / m(1) ' ускорение Меркурия по оси Z
'W(57) = (FZ(0, 2) + FZ(1, 2) + FZ(3, 2) + FZ(4, 2) + FZ(5, 2) + FZ(6, 2) + FZ(7, 2) + FZ(8, 2) + FZ(9, 2) + FZ(10, 2)) / m(2)
'W(58) = (FZ(0, 3) + FZ(1, 3) + FZ(2, 3) + FZ(4, 3) + FZ(5, 3) + FZ(6, 3) + FZ(7, 3) + FZ(8, 3) + FZ(9, 3) + FZ(10, 3)) / m(3)
'W(59) = (FZ(0, 4) + FZ(1, 4) + FZ(2, 4) + FZ(3, 4) + FZ(5, 4) + FZ(6, 4) + FZ(7, 4) + FZ(8, 4) + FZ(9, 4) + FZ(10, 4)) / m(4)
'W(60) = (FZ(0, 5) + FZ(1, 5) + FZ(2, 5) + FZ(3, 5) + FZ(4, 5) + FZ(6, 5) + FZ(7, 5) + FZ(8, 5) + FZ(9, 5) + FZ(10, 5)) / m(5)
'W(61) = (FZ(0, 6) + FZ(1, 6) + FZ(2, 6) + FZ(3, 6) + FZ(4, 6) + FZ(5, 6) + FZ(7, 6) + FZ(8, 6) + FZ(9, 6) + FZ(10, 6)) / m(6)
'W(62) = (FZ(0, 7) + FZ(1, 7) + FZ(2, 7) + FZ(3, 7) + FZ(4, 7) + FZ(5, 7) + FZ(6, 7) + FZ(8, 7) + FZ(9, 7) + FZ(10, 7)) / m(7)
'W(63) = (FZ(0, 8) + FZ(1, 8) + FZ(2, 8) + FZ(3, 8) + FZ(4, 8) + FZ(5, 8) + FZ(6, 8) + FZ(7, 8) + FZ(9, 8) + FZ(10, 8)) / m(8)
'W(64) = (FZ(0, 9) + FZ(1, 9) + FZ(2, 9) + FZ(3, 9) + FZ(4, 9) + FZ(5, 9) + FZ(6, 9) + FZ(7, 9) + FZ(8, 9) + FZ(10, 9)) / m(9)
'W(65) = (FZ(0, 10) + FZ(1, 10) + FZ(2, 10) + FZ(3, 10) + FZ(4, 10) + FZ(5, 10) + FZ(6, 10) + FZ(7, 10) + FZ(8, 10) + FZ(9, 10)) / m(10)
'End Sub





Private Sub File_AlfaP() '             запись данных расчета в файл AlfaP
JD = JD0 + Fix(100 * T3 / sekS) / 100 'для перигелия и узла по JPL0 для JDet с округлением
'можно использовать и JD=TJD, но его надо будет тоже округлить до точности записи DX(i)
i = Ipl: If kodMoon = 1 And Ipl = 10 Then i = 11
RasParamJPL0 ' на дату записи параметров, т.е. пролета перигелия (афелия)
'                  fix(5.2)=5        fix(-5.2)=-5  int(-5.2)=-6
AlfaPsum = Fix(AlfaP0(i) * kUgol / 360) 'к-во полных оборотов параметра
AlfaUsum = Fix(AlfaU0(i) * kUgol / 360)

If AlfaP0(i) > 0 Then
If AlfaP0(i) * kUgol - AlfaPsum * 360 > 300 And AlfaP < 60 Then AlfaPsum = AlfaPsum + 1
If AlfaP0(i) * kUgol - AlfaPsum * 360 < 60 And AlfaP > 300 Then AlfaPsum = AlfaPsum - 1
AlfaPsum = AlfaPsum * 360 + AlfaP
End If
If AlfaU0(i) > 0 Then
If AlfaU0(i) * kUgol - AlfaUsum * 360 > 300 And AlfaU < 60 Then AlfaUsum = AlfaUsum + 1
If AlfaU0(i) * kUgol - AlfaUsum * 360 < 60 And AlfaU > 300 Then AlfaUsum = AlfaUsum - 1
AlfaUsum = AlfaUsum * 360 + AlfaU
End If

If AlfaP0(i) < 0 Then
If AlfaP0(i) * kUgol - AlfaPsum * 360 < -300 And AlfaP > 300 Then AlfaPsum = AlfaPsum - 1
If AlfaP0(i) * kUgol - AlfaPsum * 360 > -60 And AlfaP < 60 Then AlfaPsum = AlfaPsum + 1
AlfaPsum = (AlfaPsum - 1) * 360 + AlfaP
End If
If AlfaU0(i) < 0 Then
If AlfaU0(i) * kUgol - AlfaUsum * 360 < -300 And AlfaU > 300 Then AlfaUsum = AlfaUsum - 1
If AlfaU0(i) * kUgol - AlfaUsum * 360 > -60 And AlfaU < 60 Then AlfaUsum = AlfaUsum + 1
AlfaUsum = (AlfaUsum - 1) * 360 + AlfaU
End If

'AlfaL0 здесь НЕ надо

'Rsol = Rsr(Ipl): If kodMoon = 1 And Ipl = 10 Then Rsol = 1000 * Rsr(Ipl)

NFileP = FreeFile
Open App.Path + FilePlan For Append As #NFileP
Print #NFileP, NP1, Fix(100 * T3 / sekS), Fix(AlfaPsum * 100 * MK), Fix(AlfaUsum * 100 * MK), Fix(kUgol * BettaMax * 100 * MK), Fix(Eks * 100 * MK), Fix(Rsr(Ipl) / 10), Fix(100 * MK * 2 * pi * 36525 / Period)
'все увеличил в 100*Mk, а радиус теперь в десятках метров
'все умножаем на MK и Eks еще умножаем на 1000, а Rsr(Ipl) делим на 1000, но Луна остается реальной (в Solsys6 тоже было /1000)
Close NFileP
End Sub

Private Sub File_dat() '                   запись данных расчета в файл XYZ
'kzn = 1: If kodMoon = 1 Then kzn = mk' можно будет для Луны увеличить
Info = Format(IDnew, "00.000") & "." & Format(IMnew, "00") & "." & Format(IGnew, "0000") _
& "____T= " & Format(T3 / MK, "000000000.000") _
& "____X= " & Format(X(Ipl) / MK, "000000000.000") _
& "____Y= " & Format(Y(Ipl) / MK, "000000000.000") _
& "____Z= " & Format(Z(Ipl) / MK, "000000000.000") _
& "____Rpl= " & Format(Rpl / MK, "000000000.000") _
& "____VX= " & Format(VX(Ipl), "0000000.000") _
& "____VY= " & Format(VY(Ipl), "0000000.000") _
& "____VZ= " & Format(VZ(Ipl), "0000000.000") _
& "____Vpl= " & Format(Vpl, "0000000.000") _
& "____Esys= " & Format(Esys / MK / MK / MK, "000000000.000") '____Msys=
NFileDat = FreeFile
Open App.Path + InfoFile For Append As #NFileDat 'InfoFile="\Dat\Solsys" & Ipl & "XYZ" & Nfile & ".txt"
Print #NFileDat, Info
Close NFileDat
End Sub

Private Sub Command4_Click() '                          Остановить
kodEnd = 1 'это для остановки выполнения единичного эксперимента
NPlan = 0: u = Umax 'это для остановки автоматического выполнения экспериментов по плану
End Sub

Private Sub Command5_Click() '                          Продолжить
kodEnd = 0: NC = 0
If kodTime = 1 And Option1.Value = True Then GoTo 10
If kodTime = -1 And Option2.Value = True Then GoTo 10
For k = 33 To 65: Q(k) = -Q(k): Next k ' если течение времени поменяло направление
MKD0sys = -MKD0sys
10: If Option1.Value = True Then kodTime = 1
If Option2.Value = True Then kodTime = -1
param
End Sub

Private Sub Command2_Click() '       Зделать NC0 шагов решения и остановиться
kodEnd = 1: NC = 0: NC0 = Text37.Text: param
End Sub

Private Sub Command7_Click() '     тест скорости и давления гравитации
If kodFile = 10 Then Close NfileJPL '   если остался открытым файл JPL
Form7.Show
End Sub

Private Sub Command3_Click() 'к статистической обработка данных моделирования
If kodFile = 10 Then Close NfileJPL '   если остался открытым файл JPL
Form6.Show
End Sub

Private Sub Command10_Click() '     перейти к плану эксперимента
If kodFile = 10 Then Close NfileJPL '   если остался открытым файл JPL
Form15.Show
End Sub

Private Sub HScroll1_Change()
Info = "Масштаб = " & HScroll1.Value / 10 & " млн.км/см"
Label1.Caption = Info
ML = HScroll1.Value '              Min.Value=1   Max.Value=10000
ML = 100000 * ML '         переводим в километры в в сантиметре
ML = 1000 * ML '           переводим в метры в в сантиметре
GR = Text36.Text
Risunok21
Risunok22
End Sub

Private Sub Risunok21()
Picture1.Cls '                              очищаем рисунок
Picture1.Scale (-10 + Xmouse * MLmax / ML, 12 + Ymouse * MLmax / ML)-(10 + Xmouse * MLmax / ML, -10 + Ymouse * MLmax / ML)
For i = -10 To 10 Step 1
For j = -10 To 12 Step 1
Picture1.Line (-i, j)-(i, j), RGB(70, 70, 70)
Picture1.Line (i, -j)-(i, j), RGB(70, 70, 70)
Next j
Next i
Picture1.Line (-10, 0 + GR)-(10, 0 + GR), RGB(100, 100, 100)
Picture1.Line (-10, 0)-(10, 0), RGB(100, 100, 100)
Picture1.Line (-10, 11)-(10, 11), RGB(100, 100, 100)
Picture1.Line (0, 12)-(0, -10), RGB(100, 100, 100)
For i = 0 To Ne
Picture1.Circle ((X(i) - X(Ipl0)) / ML, (Y(i) - Y(Ipl0)) / ML), 0.06, RGB(ColorR(i), ColorG(i), ColorB(i))
Picture1.Circle ((X(i) - X(Ipl0)) / ML, (Z(i) - Z(Ipl0)) / ML + 11), 0.06, RGB(ColorR(i), ColorG(i), ColorB(i))
Next i
End Sub

Private Sub Risunok22()
Picture2.Cls '                              очищаем рисунок
Picture2.Line (-5, 0)-(5, 0), RGB(100, 100, 100)
Picture2.Line (0, 4)-(0, -4), RGB(100, 100, 100)
For i = 0 To Ne
Picture2.Circle (X(i) / MLmax, Y(i) / MLmax), 0.05, RGB(ColorR(i), ColorG(i), ColorB(i))  'планеты
Next i
Picture2.Line (Xmouse - 10 * ML / MLmax, Ymouse - 10 * ML / MLmax)-Step(20 * ML / MLmax, 22 * ML / MLmax), RGB(100, 100, 100), B
End Sub

Private Sub Picture2_MouseDown(Button As Integer, Shift As Integer, X As Single, Y As Single)
If Button = 1 Then
Xmouse = X: Ymouse = Y
Risunok21
Risunok22
End If
End Sub
Private Sub Command8_Click() '         Закрыть форму
If kodFile = 10 Then Close NfileJPL
Unload Me
End Sub
Private Sub Command6_Click() '                            Выход
If kodFile = 10 Then Close NfileJPL '   если остался открытым файл JPL
Unload Me
Unload Form1
Unload Form3
Unload Form4
Unload Form5
Unload Form6
Unload Form7
Unload Form8
Unload Form9
Unload Form10
Unload Form11
Unload Form12
Unload Form13
Unload Form14
Unload Form15
Unload Form16
Unload Form17
Unload Form18
Unload Form19
Unload Form20
Unload Form21
Unload Form22
Unload Form23
End Sub


