Перевод с Pascal на C++ - C++ - Ответ 2500010

program truba_stat;
Const g=9.81; cz=300;
Type matr=array[1..2,1..2] of real;
     vect=array[1..2] of real;
     matr2=array[0..100] of vect; matr3=array[0..100] of matr;
Var a, b, c, e, obr, tmp:matr; f, vec:vect; tau, t, i,j,n:integer;
    h, m, p, ro, v, L, dt, lambda:real;
    y, beta:matr2; alpha:matr3;
    r, r1, r2,ugol,m0,p0:real;
    LR11,LR12,LR21,LR22,R11,R12,R21,R22,a11, a12, b21, b22, f1, f2:real;
Function rho(p:real):real;
begin rho:=p/(cz*cz) end;
procedure sumM(a,b:matr; var c:matr);
var i,j:integer;
begin
 for i:=1 to 2 do for j:=1 to 2 do c[i,j]:=a[i,j]+b[i,j];
end;
procedure sumv(x,y:vect; var z:vect);
var i:integer;
begin for i:=1 to 2 do z[i]:=x[i]+y[i] end;
procedure razM(a,b:matr; var c:matr);
var i,j:integer;
begin
 for i:=1 to 2 do for j:=1 to 2 do c[i,j]:=a[i,j]-b[i,j] end;
procedure razV(x,y:vect; var z:vect);
var i:integer;
begin
 for i:=1 to 2 do z[i]:=x[i]-y[i];
end;
procedure YMM(a,b:matr;var c:matr);
var i,j,k:integer; s:real;
begin
 for i:=1 to 2 do for j:=1 to 2 do
 begin s:=0; for k:=1 to 2 do s:=s+a[i,k]*b[k,j];
  c[i,j]:=s;
 end;
end;
procedure YMV(a:matr;x:vect;var y:vect);
var i,j:integer; s:real;
begin
 for i:=1 to 2 do
 begin s:=0; for j:=1 to 2 do s:=s+a[i,j]*x[j];
   y[i]:=s;
 end;
end;
Procedure obrm(a:matr; var b:matr);
var i,j:integer; d:real;
begin
  d:=a[1,1]*a[2,2]-a[1,2]*a[2,1];
  b[1,1]:=a[2,2]/d;  b[1,2]:=-a[1,2]/d;
  b[2,1]:=-a[2,1]/d; b[2,2]:=a[1,1]/d;
end;
 
begin
 assign(input,'input.txt'); reset(input);
 assign(output,'output.txt'); rewrite(output);
 readln(n,tau,m0,p0,L,Dt,lambda,ugol);
 h:=L/n;
 t:=0;
 for i:=0 to n do
 begin y[i,1]:=0; y[i,2]:=p0 end;
  r1:=tau/h; r2:=tau/(2*h);
 REPEAT  t:=t+tau;
  m:=y[0,1]; p:=y[0,2]; ro:=rho(p); v:=m/ro;
 
  R11:=1; R12:=-(v+cz)/(cz*cz);
  LR11:=v-cz; LR12:=-(v*v-cz*cz)/(cz*cz);
 
  c[1,1]:=R11-r1*LR11; c[1,2]:=R12-r1*LR12;
  c[2,1]:=0; c[2,2]:=1;
 
  b[1,1]:=-r1*LR11;  b[1,2]:=-r1*LR12;
  b[2,1]:=0; b[2,2]:=0;
 
  f1:=-ro*(g*sin(ugol)+lambda*v*abs(v)/(2*Dt));
 
  f[1]:=tau*F1+R11*m+R12*p; f[2]:=p0;
 
  obrm(c,obr); YMM(obr,b,alpha[1]); YMV(obr,f,beta[1]);
  for i:=1 to N-1 do
  begin
    m:=y[i,1]; p:=y[i,2]; ro:=rho(p); v:=m/ro;
 
    R11:=1; R12:=-(v+cz)/(cz*cz);
    R21:=1; R22:=-(v-cz)/(cz*cz);
 
    LR11:=v-cz; LR12:=-(v*v-cz*cz)/(cz*cz);
    LR21:=v+cz; LR22:=-(v*v-cz*cz)/(cz*cz);
 
    f1:=-ro*(g*sin(ugol)+lambda*v*abs(v)/(2*Dt));
    f2:=-ro*(g*sin(ugol)+lambda*v*abs(v)/(2*Dt));
 
    a[1,1]:=r2*LR11; a[1,2]:=r2*LR12;
    a[2,1]:=r2*LR21; a[2,2]:=r2*LR22;
 
    b[1,1]:=-a[1,1];  b[1,2]:=-a[1,2];
    b[2,1]:=-a[2,1];  b[2,2]:=-a[2,2];
 
    c[1,1]:=R11; c[1,2]:=R12;
    c[2,1]:=R21; c[2,2]:=R22;
 
    f1:=-ro*(g*sin(ugol)+lambda*v*abs(v)/(2*Dt));
    f2:=-ro*(g*sin(ugol)+lambda*v*abs(v)/(2*Dt));
 
    f[1]:=tau*f1+R11*m+R12*p; f[2]:=tau*f2+R21*m+R22*p;
 
    YMM(a,alpha[i],tmp); razM(c,tmp,e); obrM(e,obr);
    YMM(obr,b,alpha[i+1]); YMV(a,beta[i],vec);
    sumV(f,vec,f);
    YMV(obr,f,beta[i+1]);
  end;
 
  m:=y[n,1]; p:=y[n,2]; ro:=rho(p); v:=m/ro;
 
  R21:=1; R22:=-(v+cz)/(cz*cz);
 
  LR21:=v+cz; LR22:=-(v*v-cz*cz)/(cz*cz);
 
  c[1,1]:=1; c[1,2]:=0;
  c[2,1]:=R21+r1*LR21; c[2,2]:=R22+r1*LR22;
 
  a[1,1]:=0;         a[1,2]:=0;
  a[2,1]:=r1*LR21;  a[2,2]:=r1*LR22;
 
  f2:=-ro*(g*sin(ugol)+lambda*v*abs(v)/(2*Dt));
 
  f[1]:=m0; f[2]:=tau*F1+R11*m+R12*p;
  YMM(a,alpha[n],tmp); razM(c,tmp,e); obrM(e,obr);
  YMV(a,beta[n],vec);   sumV(f,vec,f);
  YMV(obr,f,y[n]);
 
  for i:=n-1 downto 0 do
  begin YMV(alpha[i+1],y[i+1],vec); sumV(vec,beta[i+1],y[i]) end;
  writeln('t=',t);
  for i:=0 to n do write(y[i,1]:12:2); writeln;
  for i:=0 to n do write(y[i,2]:12:2); writeln;
 UNTIL t>=600;
 close(input); close(output);
end.