{+doc
+NOMBRE: Formaec4
+CREACION:  16.05.94
+AUTORES: rch
+REGISTRO:
+TIPO: Unidad Pascal.
+PROPOSITO:  formador de ecs. de SiMEEP. M‚todo (AT-0.5I)^(-1)*(AT+0.5I).
+PROYECTO:

+REVISION: 15.11.97
+AUTOR:
+DESCRIPCION:
-doc}

unit Formaec4;

interface

uses
  xMatDefs, variable, FilUSR, horrores, CYTipos, cProcs, MatReal;

procedure Chequeo;
procedure dimencione;
procedure ObtenerEcuaciones(PasoDeTiempo: NReal);

implementation

{$IFDEF testFormaec }
procedure muestra(piv, red: integer);
var
  k, j: word;
begin
  writeln(lst, '------------');
  writeln(lst, 'piv = ', piv: 8, 'redundancia = ', red: 8);
  for k := piv to n do
  begin
    writeln(lst);
    for j := piv to nm do
      if abs(m[k, j]) > 1e-20 then
        Write(lst, 'x')
      else
        Write(lst, 'o');
  end;
  writeln(lst);
end;

{$ENDIF}

procedure limpie; forward;


procedure borre(var x: TVectFila);
var
  k: integer;
begin
  for k := 1 to nm do
    x[k] := 0;
end;


procedure compile;
begin
  chek := False;
  nve := 0;
  nd := 0;
  limpie;
  cirt1;
  nvecirt1 := nve;
  cirt2;
end;

procedure dimencione;
var
  k: integer;
begin
  for k := 1 to qnx do
    indx[k] := k;
  v[0] := 0;
  nm := qnm;
  nx := 1;
  n := qn;
  compile;
  inicio;{---}
  nx := nve;
  ndiodos := nd;
  n := nn + nx;
  nm := n + nx + nindep + 1;
end;

procedure chequeo;
begin
  chek := True;
  nd := 0;
  nve := nvecirt1;
  cirt2;
end;


procedure limpie;
var
  k, j: integer;
begin
  for k := 1 to n do
    borre(m[k]);
  { 20.2.98 cambio inicializando la lista para las ecuaciones }
  for j := 0 to n do
    arbv[j] := -1; {for j:=0 to nn do arbv[j]:=-1;}
  for k := 1 to nx do
  begin
    x[k] := x0[indx[k]];
    indx[k] := k;
  end;
  x0 := x;
end;


procedure CopiarDerivada(var f: TVectEstReal);
var
  j: word;
begin
  for j := 1 to nx do
    f[nn + j] := f[nn + j + nx];
end;



procedure IntercaReal(var x, y: NReal);
var
  z: NReal;
begin
  z := x;
  x := y;
  y := z;
end;


procedure AlteracionDeMatrices(var a: TMatSistema; deltaT: NReal);
var
  Ma, MB: TMatR;
  deltaTdiv2: NReal;
  pf, pc1, pc2: integer;
  k, j: integer;
  m2: NReal;
  res: NReal;
  res_invertible: boolean;
  res_exp10: integer;

begin
  if (nx - redundancia) = 0 then
    exit;
  deltaTdiv2 := deltaT / 2;
  pf := nn;
  pc1 := nn + nx;
  pc2 := nn + nx + nx;

  MB := TMatR.Create_init(nx - redundancia, nm - redundancia - nx - nn);
  Ma := TMatR.Create_init(nx - redundancia, nx - redundancia);
  for k := 1 to Ma.nf do
  begin
    for j := 1 to Ma.nc do
    begin
      m2 := deltaTdiv2 * a[pf + k, pc1 + j];
      Ma.Pon_e(k, j, -m2);
      MB.Pon_e(k, j, m2);
    end;
    m2 := 1;
    Ma.acum_e(k, k, m2);
    MB.acum_e(k, k, m2);
    for j := Ma.nc + 1 to MB.nc do
    begin
      m2 := a[pf + k, pc2 + j - Ma.nc];
      MB.pon_e(k, j, m2);
    end;
  end;

  res := Ma.Escaler(MB, res_invertible, res_exp10);
  if EsCero(res) then
    error('FormaEc4.AlteracionDeMatrices, no puedo escarlerizar cod7');

  for k := 1 to MB.nf do
  begin
    MB.Pon_e(k, k, MB.e(k, k) - 1);

    for j := 1 to Ma.nc do
    begin
      m2 := MB.e(k, j) / deltaT;
      a[pf + k, pc1 + j] := m2;
    end;

    for j := Ma.nc + 1 to MB.nc do
      a[pf + k, pc2 + j - Ma.nc] := MB.e(k, j);
  end;

  MB.Free;
  Ma.Free;
end;



procedure escaler(var a: TMatSistema);

const
  Cero = 1.0e-20; (* Estudiar un poco donde poner la definicion *)

var
  piv, pm: word;
  temp: NReal;
  pc1, pf1: word;


  procedure ElimineRedundantes;
  var
    pfr: integer;  (* indice a la primer fila redundancia *)
    pcr: integer; (* indice a primer columna de redundancia *)
    k, j: integer;

  begin
    pfr := nn + nx - redundancia + 1;
    pcr := nn + nx + nx - redundancia + 1;

    { Limpiamos los "ceros sucios" dejados por escaler 1}
    for k := nn + 2 to pfr - 1 do
      for j := nn + 1 to k - 1 do
        a[k, j] := 0;

    { Derivamos la relación entre las redundantes }
    for k := pfr to pf1 do
    begin
      for j := pf1 to pf1 + nx do
      begin
        a[k, j - nx] := a[k, j];
        a[k, j] := 0;
      end;

      { la derivada de las entradas constantes es CERO }
      a[k, pf1 + nx + 1] := 0;

      { Detectamos si es necesario derivar las entradas }
      for j := nn + nx + nx + 2 to pc1 do
        if not EsCero(a[k, j]) then
          error('FormaEc4.escaler, ES NECESARIO DERIVAR UNA ENTRADA');
    end;

    redundancia := 0;

    { retroceso del pivote }
    piv := nn + 1;
  end;


  function PruebaPivote: boolean;
  var
    Estado: (Probando, Falso, Bueno);
  begin
    Estado := Probando;
    repeat
      if piv > pf1 - redundancia then
      begin
        PruebaPivote := False;
        exit;
      end;

      pm := PointMaxCol(a, piv, pf1, piv);
      if piv <> pm then
        IntercaFila(a[piv], a[pm], piv, pc1);

      temp := a[piv, piv];
      if abs(temp) < Cero then
        if piv > nn then
        begin
          Inc(redundancia);
          if piv <= pf1 - redundancia then
          begin
            IntercaCol(a, piv, nn + nx - redundancia + 1, 1, pf1);
            IntercaCol(a, piv + nx, nn + nx + nx - redundancia + 1, 1, pf1);
            pm := indx[piv - nn];
            indx[piv - nn] := indx[nx - redundancia + 1];
            indx[nx - redundancia + 1] := pm;
          end
          else
            Estado := falso;
        end
        else
          error('FromaEc4.PruebaPivote, CIRCUITO INCONEXO')  {cir. inconexo}
      else
        Estado := Bueno;
    until Estado <> Probando;
    PruebaPivote := Estado = Bueno;
  end;

var
  salida: boolean;

begin
  {----------Escaler1------------------}
  piv := 1;
  redundancia := 0;

  pc1 := nn + 2 * nx + nindep + 1;
  pf1 := nn + nx;

  salida := False;

  { Mientras se pueda realizar la escalerización directa }
  while PruebaPivote do
  begin
    Eliminacion1(a, piv, pf1, pc1);
    Inc(piv);
  end;

  { Si se produjo redundancia durante la escalerización directa
  tratar la matriz para eliminar la redundancia }
  if redundancia > 0 then
  begin
    ElimineRedundantes;
    { Mientras se pueda realizar la escalerización directa }
    while PruebaPivote do
    begin
      Eliminacion1(a, piv, pf1, pc1);
      Inc(piv);
    end;
  end;

  if redundancia > 0 then
    error('FromaEc4.escaler, NO PUEDO RESOLVER LA REDUNDANCIA');
  { no puedo resolver la redundancia }

  { inicio de ESCA2 }
  for piv := pf1 downto 2 do
  begin

    if abs(a[piv, piv]) > Cero then
      DivFila(a[piv], nn + nx + 1, pc1, a[piv, piv])
    else
      error('FormaEc4.escaler, PIVOTE NULO en estapa ESCA2 !IMPOSILBE!'); { imposible }

    for pm := 1 to piv - 1 do
      if abs(a[pm, piv]) > Cero then
        SumFilaFac(a[pm], a[piv], nn + nx + 1, pc1, -a[pm, piv]);
  end;

  DivFila(a[1], nn + nx + 1, pc1, a[1, 1]);

  { Fin de ESCA2 }

end;{fin escaler}


procedure ObtenerEcuaciones(PasoDeTiempo: NReal);
begin
  Compile;
  Escaler(m);

{ Ahora modificamos el sistema de la ecuaci¢n de estado
para tener en cuenta el m‚todo de integraci¢n del trapesio
para el paso DeltaT y luego arreglamos todo para que
se pueda seguir utilizando el algoritmo de integraci¢n  de Simul
usado para Euler Xk+1 = (AXk+BRk)*dt +Xk }
  AlteracionDeMatrices(m, PasoDeTiempo);

end;

end.(* FormaEc *)
