Program prog1;
Uses dos,crt,graph;

Type
      Real = Double;
Const
      m    = 98000;    { масса ЛА                                           }
      s    = 201;      { характерная площадь несущей поверхности            }
      fi   = 3*pi/180; { угол установки двигателя                           }
      Vmin = 91.67;    { начальная скорость, м/c                            }
      Vmax = 250.00;   { конечная  скорость, м/c                            }
      Hmin = 500;      { начальная высота, м                                }
      Hmax = 8000;     { конечная высота, м                                 }
      Cyo = -0.25;
Var i,j,maxdiag,dl   : byte; { i,j-переменные цикла; maxdiag- количество    }
                             { точек на максимальной диагонали;             }
                             { dl-количество точек на текущей диагонали     }
    kmaxdiag,k,st,str: byte; { kmaxdiag-количество максимальных диагоналей  }
                             { k- счетчик диагоналей максимальной длины;    }
                             { st- номер                                    }
                             { столбца текущей точки; str- номер строки     }
                             { текущей точки                                }
    stn,strn         : byte; { stn- номер столбца начальной точки текущей   }
                             { диагонали; strn- номер строки начальной      }
                             { точки текущей диагонали                      }
    nH,nV,q          : byte; { nH- количество разбиений по высоты           }
                             { nV- количество разбиений по скорости         }
    SinT,dV,dH,Vn,Hn,min  : double; { dV- элементарное приращение высоты         }
                               { dH- элементарное приращение скорости       }
                               { Vn- скорость в текущей точке               }
                               { Hn- высота в текущей точке                 }
                               { min- минимальный путь в текущей точке      }
    nap              : array[1..80,1..80] of 'a'..'d'; { массив содержащий  }
                             { напраление перехода из текущей точки в       }
                             { следущую                                     }
    znach            : array[1..80,1..80] of double; { массив содержащий    }
                             { время перехода из текущей точки в            }
                             { конечную                                     }
    grDriver,grmode  : integer; { grDriver,grmode- параметры, определяющие  }
                                { графический режим                         }

Procedure vvod; { процедура ввода количества разбиений по высоте и скорости }
 begin
 write('Введите кол-во разбиений по высоте nH= ');
 readln(nH);
 write('Введите кол-во разбиений по скорости nV= ');
 readln(nV);
end;
  {*************************************************************}

 Function g(H:Real):Real;
  begin
  g:=9.8066-0.00000308*H;
  end;

 Function P(H : Double): Double;
  begin
  P:=3*0.8*(66500-3.27*H);
  end;

 Function a(H:Real):Real;
  begin
  a:=340.294-0.00395*H;
  end;

 Function Stream(H,V:Real):Real;
  var
   ro: Double;
  begin
   ro:=1.815-sqrt((H+2131.723)/6125.642);
   Stream := ro * sqr( V ) * S / 2;
  end;

  Function Maxa(H,V:Real): Real;
    begin
   Maxa:=V/a(H)
  end;

 Function Cya(H,V:Real):Real;
  begin
   Cya:= 0.125+1.224*Maxa(H,V)
  end;

 Function ACy(H,V:Real):Real;
  begin
   ACy:=0.0586+0.124*sqr(Maxa(H,V)-0.306)
  end;

 Function Attac_Angle( H, V : Real ) : Real;
  { Расчет угла атаки }
  { Задается H [м], V [м/c], результат в  [рад] }
  begin
   Attac_Angle := (m*g(H)-Cyo*Stream(H,V))/(Cya(H,V)*Stream(H,V));
  end;

 Function Cy(H,V:Real):Real;
  begin
   Cy:=Cyo+Cya(H,V)*Attac_Angle(H,V)
  end;
  Function Cxo(H,V:Real):Real;
  begin
   Cxo:=0.081+0.008*Maxa(H,V)
  end;
 Function Cx(H,V:Real):Real;
  begin
   Cx:=Cxo(H,V)+ACy(H,V)*sqr(Cy(H,V))
  end;

  Function X(H,V:Real) : Real;
   { Вычисление силы лобового сопротивления }
   { Задается H [м], V [м/c], результат в [H] }
   begin
      X := Cx(H,V)*Stream(H,V);
   end;

  Function SInTeta(H,V:Real) : Real;  { Синус угла наклона траектории }
   var
    SinTeta1,SinTeta2,SinTeta3 : Real;
    A,B,C,D : Real;
    P:Double;
   begin
      P:=3*0.8*(66500-3.27*H);
      A:=-ACy(H,V)*sqr(m)*sqr(g(H))/sqr(Stream(H,V));
      B:=-2*m*g(H);
      C:=2*P-(-Stream(H,V)*Cxo(H,V))+(4*ACy(H,V)+P*P-2*P*m*g(H))/(4*sqr(Stream(H,V)));
      D:=sqr(B)-4*A*C;
      SinTeta1:=(-B+sqrt(D))/(2*A);
      SinTeta2:=(-B-sqrt(D))/(2*A);
      IF SinTeta1 > 0 Then SinTeta3:=SinTeta1
      Else SinTeta3:=SinTeta2;
    {  SinTeta3:=-(X(H,V)-P)/(m*g(H));  }
      if (SinTeta3=0) Then begin Writeln('0000000000');Readkey end;
      if (SinTeta3<0) Then begin Writeln('-------------');Readkey end;
     SinTeta:=SinTeta3
   end;


  Function Power(H,V1:Real):Real;
  begin
   Power:=X(H,V1)+m*g(H)*SinTeta(H,V1)
  end;

  Function Horizontal( V1, V2, H : Real ) : Real;
   { Вычисление времени чистого разгона от V1 до V2 }
   { V1-начальная скорость [м/с], V2-конечная скорость [м/с], H-высота [м] }
   var
    P:Double;
   begin
    P:=3*0.8*(66500-3.27*H);
    Horizontal := m * ( V2 - V1 ) / ( P - X(H, V1) );
   end;

  Function Vertical( H1, H2, V1 : Real ) : Real;
   { Вычисление времени чистого подъема с H1 до H2 }
   { H1-начальная высота [м], H2-конечная высота [м], V - скорость [м/с] }
   begin
    Vertical := ( H2 - H1 ) / (V1*SinTeta(H1,V1));
   end;

  Function Diagonal( H1, H2, V1, V2 : Real ) : Real;
   { Вычисление времени подъема с H1 до H2 с разгоном от V1 до V2 }
   { H1, H2 - высоты [м], V1, V2 - скорости [м/с] }
   var
    K:Real;               { Тангенс угла наклона траектории }

   begin
    K := ( H2 - H1 ) / ( V2 - V1 );                   { вычисление тангенса }
    Diagonal := K * ln( V2 / V1 ) /SinTeta(H1,V1);
   end;


Begin

clrscr;
vvod;
dV:=(Vmax-Vmin)/nV;
dH:=(Hmax-Hmin)/nh;

Writeln('dV= ',dV:5:2);
Writeln('dH= ',dH:5:2);
readln;
IF nV<nH then maxdiag:=nV+1 else maxdiag:=nH+1;{ если количество разбиений  }
                   { по скорости меньше количества разбиений по высоте, то  }
                   { количество точек на максимальной диагонали = числу     }
                   { разбиений по высоте + 1, иначе числу разбиений по      }
                   { скорости + 1                                           }
znach[nh+1,nv+1]:=0;
dl:=1;
kmaxdiag:=abs(nV-nH)+1;
k:=0;
strn:=nh+1;
stn:=nv+1;

For i:=1 To (nV+nH) Do    { Переход по диагоналям            }
Begin
                          { Считаем длину текущей диагонали }


 If (dl<maxdiag) And (k<kmaxdiag)   Then dl:=dl+1 Else
 If (dl=maxdiag) And (k<kmaxdiag-1) Then k:=k+1   Else
                                         Begin
                                          k:=k+1;
                                          dl:=dl-1;
                                         End;
                        { Если длина меньше максимальной диагонали и     }
                        { счетчик меньше кол-ва мах. диагоналей, то      }
                        { увеличиваем длину на единицу                   }
                        { Если длина равняется максимальной диагонали    }
                        { и счетчик меньше кол-ва мах. диагоналей -1,    }
                        { то счетчик диагоналей увеличивается на 1 иначе }
                        { счетчик увеличиваем на 1, длину уменьшаем на 1 }

 If stn=1 Then strn:=strn-1;
 If stn<>1 Then stn:=stn-1;

                        { Если столбец = 1 то умен. текущую строку на 1 }
                        { Если столбец <> 1 то умен. его на 1 }
 For j:=1 To dl  Do     { переход по текущей диагонали }
  Begin
   st:=stn-1+j;
   str:=strn+1-j;

   Vn:=Vmin+dV*(st-1);
   Hn:=Hmin+dH*(str-1);
   If (Vn>750+0.0214*Hn) Then
                         { Если точка текущая не лежит в области полета ЛА, }
                         { то она исключается из рассмотрения иначе считаем }
                         { время перехода по соответствующим направлениям   }
    Begin
      nap[str,st]:='d';
      znach[str,st]:=0;
    End
   Else

   If str=nH+1 Then
             { Если точка - крайняя правая, то идем вверх }
                    Begin
                     znach[str,st]:=Horizontal( Vn, Vn+dV, Hn )+znach[str,st+1];
                     nap[str,st]:='a';
                    End
               Else If st=nV+1 Then
               { Если точка - крайняя правая, то идем вверх }
                             Begin
                              znach[str,st]:=Vertical( Hn, Hn+dH, Vn )+znach[str+1,st];
                              nap[str,st]:='c';
                             End
               Else Begin
                     If nap[str,st+1]<>'d' Then
     { Если правая точка от текущей входит в заданную полетную область то
          считаем время перехода }
                      Begin
                       znach[str,st]:=Horizontal( Vn, Vn+dV, Hn )+znach[str,st+1];
                       nap[str,st]:='a';
                       If nap[str+1,st+1]<>'d'Then
 { Если правая верхняя точка от текущей входит в заданную полетную область то
          считаем время перехода }
                        Begin
                         min:=Diagonal( Hn, Hn+dH, Vn, Vn+dV )+znach[str+1,st+1];
                         If min<znach[str,st] Then Begin
                       { Берем минимальный путь из двух }
                                                    znach[str,st]:=min;
                                                    nap[str,st]:='b';
                                                   End;
                         If nap[str+1,st]<>'d' Then
  { Если верхняя точка от текущей входит в заданную полетную область то
          считаем время перехода }
                          Begin
                           min:=Vertical( Hn, Hn+dH, Vn )+znach[str+1,st];
                           If min<znach[str,st] Then Begin
                  { Берем минимальный путь из двух }
                                                      znach[str,st]:=min;
                                                      nap[str,st]:='c';
                                                     End;
                          End;
                        End
                       Else
                        Begin
                        End;
                      End
                     Else
                      Begin
                       If nap[str+1,st+1]<>'d'Then
        { Если правая точка от текущей не входит в область то проверяем
          точку на диагонали }
                        Begin
                         znach[str,st]:=Diagonal( Hn, Hn+dH, Vn, Vn+dV )+znach[str+1,st+1];
                         nap[str,st]:='b';
                         min:=Vertical( Hn, Hn+dH, Vn )+znach[str+1,st];
                         If min<znach[str,st] Then
        { Берем минимальный путь из двух }
                          Begin
                           znach[str,st]:=Vertical( Hn, Hn+dH, Vn )+znach[str+1,st];
                           nap[str,st]:='c';
                          End;
                        End
                       Else
   { Если правая и правая верхняя точки не лежат в области - идем вверх }
                        Begin
                         znach[str,st]:=Vertical( Hn, Hn+dH, Vn )+znach[str+1,st];
                         nap[str,st]:='c';
                        End;
                      End;
                    End;
  End;
End;

Writeln('-----------------');

q:=0;
i:=1;
j:=1;

While (i<>nh+1) Or (j<>nV+1) Do
Begin
   If nap[i,j]='a' Then Begin
                         Write('a[',i,',',j,']=',znach[i,j]:5:2);
                         Write(nap[i,j],'  ');
                         i:=i;

                         If j<>nv+1 Then j:=j+1;

                        q:=q+1;
                        If q/4=trunc(q/4) Then Writeln;

                        End;
   If nap[i,j]='b' Then Begin
                         Write('a[',i,',',j,']=',znach[i,j]:5:2);
                         Write(nap[i,j],'  ');
                         If i<>nh+1 Then i:=i+1;
                         If j<>nv+1 Then j:=j+1;
                          q:=q+1;
                        If q/4=trunc(q/4) Then Writeln;

                        End;
   If nap[i,j]='c' Then Begin
                         Write('a[',i,',',j,']=',znach[i,j]:5:2);
                         Write(nap[i,j],'  ');
                         If i<>nh+1 Then i:=i+1;
                         j:=j;
                          q:=q+1;
                        If q/4=trunc(q/4) Then Writeln;

                        End;
End;

Readln;


grDriver := Detect;       { Вывод графиков }
InitGraph(grDriver, grMode,'f:\work\bp7\bgi');
Setbkcolor(black);
Cleardevice;
SETCOLOR(white);

dv:=590/nv;
dh:=400/nh;

LINE(20,20,20,420);
LINE(20,20,610,20);

LINE(20,420,610,420);
LINE(610,20,610,420);

Setlinestyle(0,1,3);

i:=1;
j:=1;

While (i<>nh+1) or (j<>nV+1) do
Begin
   If nap[i,j]='a' Then Begin
                         line(20+round((j-1)*dv),420-round((i-1)*dh)
                         ,20+round(j*dv),420-round((i-1)*dh));
                         i:=i;
                         If j<>nv+1 Then j:=j+1;
                        End;
   If nap[i,j]='b' Then Begin
                         line(20+round((j-1)*dv),420-round((i-1)*dh)
                         ,20+round(j*dv),420-round(i*dh));
                         If i<>nh+1 Then i:=i+1;
                         If j<>nv+1 Then j:=j+1;
                        End;
   If nap[i,j]='c' Then Begin
                         line(20+round((j-1)*dv),420-round((i-1)*dh)
                         ,20+round((j-1)*dv),420-round(i*dh));
                         If i<>nh+1 Then i:=i+1;
                         j:=j;
                        End;
End;


REPEAT
UNTIL KEYPRESSED;
End.