unit D3;
interface
uses graph, dos, crt;
const zero: array[1..3]of longint=(0, 0, 0);
      V0: array[1..3]of real=(0, 0, 0);
type DirVec=1..3;
     Vector = array[DirVec] of real;
     Point  = array[DirVec] of longint;
     sPoint = array[DirVec] of integer;

var O    : Point;   {Pozitia in Oxyz}
    Grvt : Point;   {Vectorul gravitatiei}
    Lght : Point;
    Mv   : Vector;  {Vectorul acceleratiei}
    speed: integer; {Viteza de miscare}
    deeph: real;    {Adancimea proiectiei centrale}
    a, b : real;    {Unghiul orizontal si vertical de privire}
    miu  : real;    {Coeficientul de frecare}
    v: array[DirVec, DirVec] of real; {Matricea operatorului de transformare a coordonatelor}
    sa, sb, ca, cb: real;  {sin si cos}

    TLi: Longint; TR: real; TS: String; TP: Point; {Temp}

{ //--- Calc Coord --- }
procedure CalcSC; {Calcularea Sin Cos}
procedure CalcVec;
function  PCrd(p: point; d: dirvec): longint; Function VCrd(p: vector; d: dirvec): real;
Procedure PRot(pr: point; var pv: point);    Procedure VRot(pr: vector; var pv: vector);
{ --- Calc Coord ---// }
Procedure ProdVec(v1, v2: Point; var res: Point);                         {Produsul vectorial}
Procedure MiniVec(p: Point; var res: vector);                             {Se obtine un vector cu coordonata maxima 1, paralel la p}
Procedure vProdVec(v1, v2: Point; var res: vector);                       {Produsul vectorial, rezultatul de tip Vector, coordonata maxima 1}
Procedure pGetPerp(p1, p2, p3: Point; var res: Point);                    {Se obtine un vector perpendicular la planul determinat de punctele p1, p2, p3}
Procedure vGetPerp(p1, p2, p3: Point; var res: Vector);                   {Ca si pGetPerp, cu coordonata maxima 1}
Procedure pGetPerpL(p1, p2, p3: Point; var res: Point; lng: longint);     {Ca si pGetPerp, cu modulul lng}
Procedure pGetPerpPoint(p1, p2, p3: Point; var res: Point; lng: longint); {Punctul situat la distanta lng de la planul determinat de p1, p2, p3}
Procedure GetVecPar(v1: vector; L: LongInt; var res: point);              {Obtinerea unui vector de tip Point, de modul L, paralel cu v1}
Procedure GetVec(p1: Point; var v1: Vector);                              {Obtinerea unui vector de modul 1 paralel la p1}
Procedure VecUnit(v1: vector; var res: vector);                           {Obtinerea unui vector cu modulul 1 din v1}
Procedure PSum(v1, v2: Point; var res: Point);
Procedure PDif(v1, v2: Point; var res: Point);
Procedure vSum(v1, v2: Vector; var res: Vector);
Procedure vDif(v1, v2: Vector; var res: Vector);
Function  pProd(v1, v2: Point): Longint;
Function  vProd(v1, v2: Vector): real;
Function  vpProd(v1: vector; p1: Point): Extended;
function  r(g: real): real;
function  g(r: real): real;
Function  DphP(v: real): real;
Procedure Move(d: DirVec; st: integer);
Procedure vMove(p: Vector);
function  Gravity(var Spd: Vector; var poz: Point; zer: longint): longint;
Procedure CalcDynam;
Function  Dist(p1, p2: vector): real;
Function  vModul(p: Vector): extended;      Function PModul(p: point): real;
Procedure vWrite(v1: vector);
Procedure pWrite(v1: point);

implementation
var i, j, k: LongInt;
    bt: byte;

Procedure vWrite;
begin
for bt:=1to 3 do write(v1[bt]: 10: 2, ' '); writeln;
end;

Procedure pWrite;
begin
for bt:=1to 3 do write(v1[bt]: 10, ' '); writeln;
end;

Procedure GetVecPar;
var l1, l2: real;
begin
 l2:=vModul(v1);
 if(l2<>0)then begin
  l1:=l/l2;
  for bt:=1 to 3 do begin
   res[bt]:=round(v1[bt]*l1);
  end;
 end;
end;

Procedure GetVec(p1: Point; var v1: Vector);
var l: real;
begin
 l:=pModul(p1);
 if(l<>0)then begin
  for bt:=1 to 3 do begin
   v1[bt]:=(p1[bt]/l);
  end;
 end;
end;

Procedure VecUnit(v1: vector; var res: vector);
var l: real;
begin
 l:=vModul(v1);
 if(l<>0)then begin
  for bt:=1 to 3 do begin
   res[bt]:=(v1[bt]/l);
  end;
 end;
end;

Procedure ProdVec(v1, v2: Point; var res: Point);
begin
 res[1] := v1[2]*v2[3] - v2[2]*v1[3];
 res[2] := v1[3]*v2[1] - v2[3]*v1[1];
 res[3] := v1[1]*v2[2] - v2[1]*v1[2];
end;

Procedure MiniVec;
var i, j: byte;
begin
 res:=v0; j:=1;
 for i:=1 to 3 do if p[i]>p[j] then j:=i;
 if(j<3)and(p[j]<p[3])then j:=3;
 for i:=1 to 3 do res[i]:=(p[i]/p[j]);
end;

Procedure pGetPerp(p1, p2, p3: Point; var res: Point);
var t1, t2: Point;
begin
 PDif(p2, p1, t1); PDif(p3, p1, t2);
 ProdVec(t1, t2, res);
end;

Procedure vProdVec; begin ProdVec(v1, v2, tp); MiniVec(tp, res); end;
Procedure vGetPerp; begin pGetPerp(p1, p2, p3, tp); MiniVec(tp, res); end;

Procedure pGetPerpL(p1, p2, p3: Point; var res: Point; lng: longint);
var v: Vector;
begin
 vGetPerp(p1, p2, p3, v);
 GetVecPar(v, lng, res);
end;

Procedure pGetPerpPoint(p1, p2, p3: Point; var res: Point; lng: longint);
var v: Vector;
begin
 vGetPerp(p1, p2, p3, v);
 GetVecPar(v, lng, res);
 PSum(p1, res, res);
end;

function Gravity(var Spd: Vector; var poz: Point; zer: longint): longint;
var i: integer;
begin
//   if poz[2]>zer then Spd[2] := Spd[2]+grvt[2] else
//   if poz[2]<zer then begin
//     if Spd[2]<zer then Spd[2]:=-Spd[2];
//     Spd[2]:=round(Spd[2]*2/3);
//     inc(poz[2], (500-poz[2])div 2+1);
//   end;
  if poz[2]<zer then begin
    inc(poz[2], (zer-poz[2])div 2+1);
    if Spd[2]<0 then begin
      Spd[2] := abs(Spd[2]);
      Spd[2] := round(Spd[2]*2/3); {Ciocnirea consuma o parte din energia cinetica pentru a o transforma in energie termica}
    end;
  end else
  if (poz[2] = zer) and (Spd[2] < 0) then Spd[2] := -grvt[2]; {La nivelul zero acceleratia grav. se ekilibreaza}

  {Acceleratia gravitationala}
  for i := 1 to 3 do Spd[i] := Spd[i] + grvt[i];

  Result:=poz[2];
end;

Procedure CalcDynam;
begin
  {Inertia}
  vMove(Mv);

  lght := o;

  {Gravitatia}
  Gravity(Mv, o, 500);

  {Forta de frecare}
  for i:=1 to 3 do Mv[i] := Mv[i]*miu;
end;

Procedure Move(d: DirVec; st: integer);
begin
//    TP:=zero; TP[d]:=st;
//    for i:=1 to 3 do mv[i]:=mv[i]+round(v[1, i]*tp[1]+v[2, i]*tp[2]+v[3, i]*tp[3]);
   for i:=1 to 3 do mv[i]:=mv[i]+round(v[d, i]*st);
end;

Procedure vMove(p: Vector); begin for i:=1 to 3 do O[i] := round(O[i]+p[i]); end;

function r; begin r:=g*pi/180; end;
function g; begin g:=r*180/pi; end;

Function  Dist; begin dist:=sqrt(sqr(p2[1]-p1[1])+sqr(p2[2]-p1[2])+sqr(p2[3]-p1[3])); end;
Function  vModul; begin vModul:=sqrt(sqr(p[1])+sqr(p[2])+sqr(p[3])); end;
Function  PModul; begin pModul:=sqrt(abs(sqr(p[1])+sqr(p[2])+sqr(p[3]))); end;
Procedure pSum; begin for bt:=1 to 3 do res[bt]:=v1[bt]+v2[bt]; end;
Procedure pDif; begin for bt:=1 to 3 do res[bt]:=v1[bt]-v2[bt]; end;
Procedure vSum; begin for bt:=1 to 3 do res[bt]:=v1[bt]+v2[bt]; end;
Procedure vDif; begin for bt:=1 to 3 do res[bt]:=v1[bt]-v2[bt]; end;
Function  pProd; var p: longint; begin p:=0; for bt:=1 to 3 do p:=v1[bt]*v2[bt]+p; result:=p; end;
Function  vProd; var p: real; begin p:=0; for bt:=1 to 3 do p:=v1[bt]*v2[bt]+p; result:=p; end;
Function  vpProd; var p: Extended; begin p:=0; for bt:=1 to 3 do p:=v1[bt]*p1[bt]+p; result:=p; end;
Function  DphP; begin if v+deeph<>0 then DphP:=deeph/(v+deeph) else Dphp:=0; end;

{//------------------------- Calcule de transformare a coordonatelor -------------------------}
procedure CalcSC;
begin
    sa := sin(a);      sb := sin(b);
    ca := cos(a);      cb := cos(b);
end;

procedure CalcVec;
begin CalcSC;
      v[1, 1] := ca*cb;  v[2, 1] := ca*sb;  v[3, 1] := -sa;
      v[1, 2] := -sb;    v[2, 2] := cb;     v[3, 2] := 0;
      v[1, 3] := sa*cb;  v[2, 3] := sa*sb;  v[3, 3] := ca;
end;

Function  PCrd; begin PCrd:=round((p[1]-o[1])*v[d, 1]+(p[2]-o[2])*v[d, 2]+(p[3]-o[3])*v[d, 3]); end;
Function  VCrd; begin VCrd:=p[1]*v[d, 1]+p[2]*v[d, 2]+p[3]*v[d, 3]; end;
Procedure pRot; begin for i:=1 to 3 do pv[i]:=PCrd(pr, i); end;
Procedure vRot; begin for i:=1 to 3 do pv[i]:=VCrd(pr, i); end;
{------------------------- Calcule de transformare a coordonatelor -------------------------//}

Begin
     a      := 0;
     b      := 0;
     o      := zero;
     grvt   := o;
     Lght   := o;
     Lght[2]:= round(Deeph);
//     CalcVec;
     speed  := 200;
     miu    := (1-1/100);
     grvt[2]:=-speed div 2;
End.