{$A+,B-,D+,E+,F-,G-,I+,L+,N-,O-,P-,Q-,R-,S+,T-,V+,X+,Y+}
(*{$М 16384,0,655360}*)
{N+}
 Program traectoriya;
 Uses Crt, Graph;
   type
    deltat=record {тип для массива затрат времени}
    hv:array[1..3] of real; {1-при изменение Н, 2 -при изменении V ,
                             3-при изменении H,V}
   end;
   Trace=record {тип для траекторного массива}
    sum:real; {суммарная затрата времени в точке}
    n,m: byte;   {координаты точки с минимальным временем}
   end;

   const
     massa=47000;  {масса самолета ТУ-134, кг }
     S=127;     {Эффективная площадь несущих поверхностей}
     Stepm=49;
     alfa0=4*pi/180;
     fi=3*pi/180;
   Var {описание глобальных переменных}
    Step:integer;
    Hn,Hk,Vn,Vk:real;
    dV,dH,Vt,Ht,Vnv,Vkv:real;{переменные нужны для изменения скорости и высоты}
    masT:array[1..Stepm,1..Stepm] of deltat ;  {Массив временного расхода}
    i,j,k,num: integer; {счетчики}
    mastr:array[1..Stepm,1..Stepm] of Trace;{массив траекторий}
    vsmas : array [1..3] of real;{вспомогательный массив временного расхода}
    min:real;
    nexti,nextj:integer;
    xn,yn,dx,dy,xn1,xn2,yn1,yn2:integer;
    gdetect,gmode:integer;
    H_string, V_String: String;

Procedure graf(st:integer);
var x,y:integer;
 begin
      xn:=70; yn:=400;
      dx:=round((620-xn)/Step);
      dy:=round((400-50)/Step);
      settextstyle(0,0,0);
      setcolor(3);
      x:=xn+step*dx;
      y:=yn-step*dy;
      for i:=1 to Step+1 do {Сетка}
       begin
         line(xn,yn,x,yn);
         yn:=yn-dy ;
       end;
      for j:=1 to Step+1 do
       begin
        line(xn,y,xn,400);
        xn:=xn+dx ;
       end;

    {Вывод текста}
   setcolor(2);
   outtextxy(100,25,'Для выхода из программы нажмите любую клавишу');
   setcolor(7);
   xn1:=50; xn2:=10;
   yn1:=410; yn2:=400;
   Str(Hn:4:1, H_string);
   Str(Vnv:4:1, V_string);
   outtextxy(xn1,yn1,V_string);
   outtextxy(xn2,yn2,H_string);
        xn1:=xn1+Step*dx-20;
        yn2:=yn2-Step*dy;
        Str(Hk:4:1, H_string);
        Str(Vkv:4:1, V_string);
        outtextxy(xn1,yn1,V_string);
        outtextxy(xn2,yn2,H_string);
  end;


    Function g(H:real ):real;  {Ускорение свободного падения}{OK}
      begin
        g:=9.8066-0.000003*H;
      end;

    Function Sound_speed(H:real):real;{скорость звука}
      begin
        Sound_speed:=340.294-(27.988/7000)*H;
      end;

    Function density(H:Real): Real; {Расчет плотности} {OK}
      begin
       density:=1.225*exp(-0.00009654630994*H);
      end;

    Function Power(H,V:Real):Real;{Расчет силы тяги двигателя}
     var
       M:Real;    { Число Маха  }
       P:Real;    { Вспогательная переменная }
     begin
       M:=V/Sound_speed(H); {Вычисления числа маха}
       If (H>0)and(H<=4000) then P:=77500-20000*M; {по графику для Н=О }
       If (H>4000) and (H<=6000) then P:=71666.7-21666.7*M;    {для H=4000 }
       If (H>6000) and (H<=8000) then P:=55833-11666.6*M; {для Н=6000 }
       if (H>8000) and (H<=10000) then P:=45000-4444.533*M; {для H=8000}
       Power:=P*1.05;     {По условию Max P 105% от номинальной}
     end;
   Function kapa(H,V,dH,dV:real):real;
     begin
       if dH=0 then kapa:=0 else  kapa:=1/(1+V*dV/g(H)/dH);
     end;

   Function Cx0(H,V:real):real;
     var
       M:real;
     begin
       M:=V/Sound_speed(H);
       Cx0:=0.016*M+0.007;
     end;

   Function A(H,V:real):real; {коэффициент развала поляры}
     var
       M:real;
     begin
       M:=V/Sound_speed(H);
       A:=-0.07*M+0.096;
     end;

   Procedure itera(var alfa1, teta1,Cy,Cx,X,H,V,dH1,dV1:real);
   {расчет альфа, тэта, и т.д. методом итераций}
    var
     alpha,teta:real;
     Cy0,Xo,Yo,Y:real;
     k:integer;
      begin
       alpha:=5*pi/180;
       teta:=5*pi/180;
       Cy0:=0.1;
       Cy:=Cy0*(alpha-alfa0);
       Cx:=Cx0(H,V)+A(H,V)*Cy*Cy;
       Xo:=Cx*density(H)*V*V*S*0.5;
       Yo:=Cy*density(H)*V*V*S*0.5;
        for k:=1 to 4 do
         begin
           alpha:=(Cy0*alfa0*S*density(H)*V*V*0.5+massa*g(H)*cos(teta)-
           Power(H,V)*sin(alpha+fi))/(Cy0*S*density(H)*V*V*0.5);
           teta:=arctan((kapa(H,V,dH,dV)*Power(H,V)*cos(alpha+fi)-Xo*kapa(H,V,dH,dV))/(Power(H,V)*sin(alpha+fi)+Yo));
         end;
          alfa1:=alpha;
          teta1:=teta;
          Cy:=Cy0*(alfa1-alfa0);
          Cx:=Cx0(H,V)+A(H,V)*Cy*Cy;
          X:=Cx*S*density(H)*V*V*0.5;
        end;

   Function dt(H,V,dH,dV:real):real; {временной расход}
     var
      alfa1,teta1,Cy,Cx,X: real;
     begin
      itera(alfa1,teta1,Cy,Cx,X,H,V,dH,dV);
      if dV=0 then dt:=dH/(V*sin(teta1));
      if dH=0 then dt:=(massa*dV)/(Power(H,V)*cos(alfa1+fi)-X-massa*g(H)*sin(teta1));
      if (dV<>0)and(dH<>0) then dt:=dH/(V*sin(teta1));
     end;

   Procedure tracebuilding;
    begin
     moveto(70,400);
     xn:=70; yn:=400;
     j:=1;i:=1;
     nexti:=mastr[i,j].n;
     nextj:=mastr[i,j].m;  {получаем ссылку на cледующую точку}
      repeat
      {По полученным значениям ссылок на последукщию точку,направляем линию дальше}
       if ((nexti-i)<>0) and ((nextj-j)=0) then xn:=xn+dx ;
       if ((nexti-i)=0) and ((nextj-j)<>0) then yn:=yn-dy ;
       if ((nexti-i)<>0) and ((nextj-j)<>0) then  begin
                                                 xn:=xn+dx;
                                                 yn:=yn-dy;
                                               end;
       setcolor(5);
       lineto(xn,yn) ;
       j:=nextj; i:=nexti; {переходим в точку на которую получили ссылку}
       nextj:=mastr[i,j].m;    {получаем ссылку на cледующую точку}
       nexti:=mastr[i,j].n;
      until (j=Step+1) and (i=Step+1);
    end;



  BEGIN     {MAIN 3аполнение массива временного расхода}
  clrscr;
  write('Введите начальную высоту в метрах:');readln(Hn);
  write('Введите конечную высоту в метрах не более 11000 метров:');readln(Hk);
  write('Введите начальную скорость в км/ч:');readln(Vnv);
  write('Введите конечную скорость в км/ч не более 900 км/ч:');readln(Vkv);
  Vn:=Vnv*1000/3600;{начальная скорость [м/с] }
  Vk:=Vkv*1000/3600;
  write('Введите число разбиений не более 40:'); read(Step);
  dH:=(Hk-Hn)/Step;  {Разбиваем на части всю длину высоты и скорости}
  dV:=(Vk-Vn)/Step;
  Ht:=Hn;Vt:=Vn; {Уходим в начальную точку}


  for j:=1 to Step do  {Подсчет временного расхода при изменении высоты}
    begin
     Ht:=Hn;  {В начальную точку}
      for i:=1 to Step  do
       begin
        if (Ht<36*Vt+5318.18)and(Ht<113.142*Vt-3428.527)
           and(Ht<263.994*Vt-15999.5)
           and (Ht>247.491*Vt-43123.195)
           and (Ht>158.4*Vt-32200.94)
          then masT[i,j].hv[1]:=dt(Ht,Vt,dH,0)
          else
            masT[i,j].hv[1]:=100000 ;
             Ht:=Ht+dH;
       end;
             Vt:=Vt+dV;
     end;

   {Подсчет временного расхода при изменении скорости}
      Ht:=Hn;  Vt:=Vn; {в начальную точку}
       for i:=1 to Step  do
        begin
         Vt:=Vn;
         for j:=1 to Step do
           begin
            if (Vt<918.1818)
                and (Ht>247.491*Vt-43123.195) and (Ht>158.4*Vt-32200.94)
                and (Ht<36*Vt+5318.18)and(Ht<113.142*Vt-3428.527)
                and (Ht<263.994*Vt-15999.5)
               then masT[i,j].hv[2]:=dt(Ht,Vt,0,dV)
               else masT[i,j].hv[2]:=100000;
            Vt:=Vt+dV
            end;
               Ht:=Ht+dH;
        end;

  {Подсчет временного расхода при изменении скорости и высоты одновременно}

   Vt:=Vn;  Ht:=Hn; {в начальную точку}
   for i:=1 to step do
    begin
     Vt:=Vn;
     Ht:=Ht+dH;
      for j:=1 to Step do
       begin
        Vt:=Vt+dV;
        if (Vt<918.1818) and (Ht<11000)
           and (Ht>247.491*Vt-43123.195) and (Ht>158.4*Vt-32200.94)
           and (Ht<36*Vt+5318.18)and(Ht<113.142*Vt-3428.527)
           and (Ht<263.994*Vt-15999.5)
           then masT[i,j].hv[3]:=dt(Ht,Vt,dH,dV)
           else masT[i,j].hv[3]:=100000 ;
       end;
    end;

                 {Нахождение траектории}
   i:=Step+1; j:=Step+1;{в конечную точку}
   mastr[i,j].sum:=0;{обнуление траекторного временного расхода}
   for j:=Step downto 1 do
    begin
     mastr[i,j].sum:=mastr[i,j+1].sum+masT[i,j].hv[1];{Начинаем суммировать по горизонтали}
     mastr[i,j].n:=i;
     mastr[i,j].m:=j+1; {Запоминаем координаты предыдущей точки}
    end;

   i:=Step+1; j:=Step+1;{в конечную точку}
   mastr[i,j].sum:=0;{обнуление траекторного временного расхода}
   for i:=step downto 1 do
    begin
     mastr[i,j].sum:=mastr[i+1,j].sum+mast[i,j].hv[2] ;{Начинаем суммировать по вертикали}
     mastr[i,j].n:=i+1;
     mastr[i,j].m:=j ;  {Запоминаem координаты предыдущей точки}
    end;


   vsmas[1]:=0; vsmas[2]:=0; vsmas[3]:=0; {обнуляем вспомогательный массив}
    for i:=Step downto 1 do
     begin
      for j:=Step downto 1 do
       begin   {3 временных расхода для данной точки}
        vsmas[1]:=mastr[i,j+1].sum+masT[i,j].hv[2];
        vsmas[2]:=mastr[i+1,j].sum+masT[i,j].hv[1];
        vsmas[3]:=mastr[i+1,j+1].sum+masT[i,j].hv[3];
        min:=vsmas[1] ;
        num:=1;
         for k:=1 to 3 do
          begin
           if vsmas[k]<min then {определение минимального расхода}
            begin
             min:=vsmas[k] ; num:=k;
            end;
          end;
         mastr[i,j].sum:=min;
         if num=1 then begin mastr[i,j].n:=i; mastr[i,j].m:=j+1;end;
         if num=2 then begin mastr[i,j].n:=i+1; mastr[i,j].m:=j ;end;
         if num=3 then begin mastr[i,j].n:=i+1;mastr[i,j].m:=j+1;end;
       end;
     end;{ массив траекторий }


  clrscr;     {Построение траектории}
  gdetect:=Detect;
  InitGraph(gdetect, gmode, 'c:\progra~1\Borland\Pascal\bgi');
  if GraphResult <> grOk then Halt(1);
  graf(Step);  {Построение сетки и высчитывание шага по х и у}
  tracebuilding;

repeat until keypressed;
closegraph ;
End.