Высотные ракеты SashaMaks

//Функция расчёта Cf по давлению P0 и константе адиабаты k:
function GetCF(const P0, K: Double): Double;
begin
//Простой расчёт Cf по формуле:
	Result := Power(2 * SQR(K) * Power(2 / (K + 1), (K + 1)/(K - 1)) * (1.0 - Power(1.0 / P0, (K - 1) / K)) / (K - 1), 0.5);
end;

//Расчёт площади критического сечения:
Sk := SQR(Di * 0.0005) * Pi;
//Di - текущий диаметр для измерения;
//Pi - просто число ПИ (3,14...);
//Sk - площадь критического сечения;
//Условие выбора между расчётом для сопла Лаваля и без него:
	If Laval then begin
	//Запоминаем отдельное измерение тяги из массива данных в отдельную переменную:
	Force := FData.Y;
	//Условие защиты алгоритма от неадекватных отрицательных значений тяги:
	If Force < 0.0 then Force := 0.0;
	//Считаем верхний предел давления Pmax без учёта Cf:
	Pmax := Force * g_const / (Sk * 101325);
	//Рассчитываем Cf по функции: GetCF
	CF1 := GetCF(Pmax + 1.0, K_Adiabat.GetY(Pmax + 1.0));
	//При этом подставляется абсолютное значение давления Pmax + 1 атмосфера и берётся для этого давления своя константа адиабаты из функции зависимости K_Adiabat (k = f(P0)). Данная зависимость строиться в программе Propellant.
	//Далее делается анализ на предмет полученных значений Cf так, чтобы не было деления на 0.
	If CF1 <> 0.0 then
		//Для ненулевых значений Cf рассчитывается нижний предел давления Pmin
		Pmin := Force * g_const / (Sk * CF1 * 101325) else
		//Для нулевых значений Cf минимальное давление берётся равным максимальному
		Pmin := Pmax;

		//Алгоритм итерации
		//Берётся значение шага давления 1/100 от полученного диапазона (Pmax - Pmin) и далее выполняется последовательный подбор давления к истинному значению
		dP := (Pmax - Pmin) * 0.01;
		//Итерирование продолжается до тех пор, пока минимальное давление меньше максимального
		While Pmin < Pmax do begin
			//Минимальное давление увеличивается на dP
			Pmin := Pmin + dP;
			//Вычисляется новый Cf по формуле в функции GetCF по текущему новому минимальному давлению и новому значению адиабаты для него k
			CF1 := GetCF(Pmin + 1.0, K_Adiabat.GetY(Pmin + 1.0));
			//Вычисляется сила тяги от этого давления и Cf для него
			Force := Sk * Pmin * 101325 * CF1 / g_const;
			//После чего полученное значение силы тяги сравнивается с измеренным с заданной точностью 1%
			If ABS(Force - FData.Y) / FData.Y < 0.01 then
				Break;
			//Если разница между расчётным значением силы тяги и измеренным оказывается меньше 1/100, то итерирование прерывается, а начисленное значение Pmin записывается в функцию давления, как и Cf для вывода в графическом виде
		end;
		PData.Y  := Pmin;
		CFData.Y := CF1;
	end else begin
		//Без сопла Лаваля  расчёт идёт по упрощенной схеме
		//Запоминаем отдельное измерение тяги из массива данных в отдельную переменную:
		Force := FData.Y;
		//Условие защиты алгоритма от неадекватных отрицательных значений тяги:
		//Этот код повторяется, так как они принципиально разделены и чтобы потом не забыть, где были его отдельные части, дабы не было ошибок при его корректировке.
		If Force < 0.0 then Force := 0.0;
		//Просто вычисляется значение давления через площадь критического сечения
		PData.Y  := Force * g_const / (Sk * 101325) + 1.0;
		//Cf принимается равным 1
		CFData.Y := 1.0;
	end;