unit InverseUnit;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, Spin, ComCtrls, ShellCtrls, DBCtrls, FileCtrl,
  Math, HCArrays, HCConduct, HCMaterial, HCModel, UPohMath, UPohStrConv,
  BasicUnit, DirectUnit;

const
  mm=0.001;
  Cel_Kel=273.15;

function FindFTsteps(n: integer; InitTemp: double; L: double; MatID: integer; MatDtb: tMatDB; max_volL:double; dt:double; hloubka:double): integer;
procedure IHCP_Beck(MatDtb: tMatDB; L: double; InitTemp:double; MatID:integer; max_volL:double; n:integer; FTsteps:integer; T: array of double; var q: array of double; var T_comp: array of double; var T_surf: array of double; dt:double; hloubka:double);
procedure Sensitivity(n: integer; L: double; MatID: integer; MatDtb: tMatDB; max_volL:double; dt:double; var del_theta: array of double; hloubka:double);
procedure TempNull(var model: tHCModel; var vTbck: tFloatArray; n: integer; MatDtb: tMatDB; dt:double; var T_qnull: array of double; hloubka: double);
procedure CoefBeck(n: integer; L: double; MatID: integer; MatDtb: tMatDB; max_volL:double; dt:double; var del_theta: array of double; hloubka:double);
procedure IHCP_Full(n:integer; T: array of double; var q: array of double; del_theta: array of double);
procedure IHCP_Full_reg(n:integer; T: array of double; var q: array of double; del_theta: array of double);

implementation

function FindFTsteps(n: integer; InitTemp: double; L: double; MatID: integer; MatDtb: tMatDB; max_volL:double; dt:double; hloubka:double): integer;
var
  q1,q2, T1, T2, diff, T_surf: array of double;
  model: tHCModel;
  vTbck : tFloatArray;
  i, position:integer;
  max:double;
begin
  SetLength(q1,n); SetLength(q2,n); SetLength(T1,n); SetLength(T2,n);
  SetLength(diff,n); SetLength(T_surf,n);
  q1[0]:=-10000000; q1[1]:=0;
  q2[0]:=0; q2[1]:=-10000000;
  for i:=2 to (n-1) do begin
    q1[i]:=0; q2[i]:=0;
  end;
  InitModel(model, L, InitTemp, MatID, max_volL);
  fHC_Init(vTbck);
  SetBoundCond(model);
  fHC_ClearAccuracy(model); // clear accuracy arrays
  fHC_Copy(vTbck, model.vT); //zaloha teplotniho pole
  MainDirectTask(T1, T_surf, MatDtb, model, q1, n, dt, hloubka);
  fHC_Copy(model.vT, vTbck);
  MainDirectTask(T2, T_surf, MatDtb, model, q2, n, dt, hloubka);
  fHC_Free(model);
  for i:=0 to (n-1) do diff[i]:=Abs(T1[i]-T2[i]);
  max:=MaxValue(diff);
  for i:=0 to (n-1) do begin
    if diff[i]=max then position:=i;
  end;
  result:=position+2;
end;

procedure IHCP_Beck(MatDtb: tMatDB; L: double; InitTemp:double; MatID:integer; max_volL:double; n:integer; FTsteps:integer; T : array of double; var q: array of double; var T_comp: array of double; var T_surf: array of double; dt:double; hloubka:double);
var
  q1, T1, T_s, del_theta: array of double;
  model: tHCModel;
  vTbck : tFloatArray;
  sum1, sum2:double;
  i,j:integer;

begin
  SetLength(q1,FTsteps); SetLength(T1,FTsteps);
  SetLength(T_s,FTsteps);
  SetLength(del_theta,FTsteps);

  CoefBeck(FTsteps, L, MatID, MatDtb, max_volL, dt, del_theta, hloubka);
  for i:=0 to (FTsteps-1) do q1[i]:=0;

  // Model
    InitTemp:=T[0];
    InitModel(model, L, InitTemp, MatID, max_volL);
    fHC_Init(vTbck);
    SetBoundCond(model);
    fHC_ClearAccuracy(model);
    fHC_Copy(vTbck, model.vT); //zaloha teplotniho pole

  for i:=1 to (n-FTsteps-1) do begin
    MainDirectTask(T1, T_s, MatDtb, model, q1, FTsteps, dt, hloubka);
    sum1:=0;sum2:=0;
    for j:=0 to (FTsteps-1) do begin
      sum1:=sum1 + ((T[i+j]-T1[j])*del_theta[j]);
      sum2:=sum2 + (del_theta[j]*del_theta[j]);
    end;
    q[i]:=-sum1/sum2;
    fHC_Copy(model.vT, vTbck); //vlozeni zalohovaneho teplotniho pole
    MainDirectTask(T_comp[i], T_surf[i], MatDtb, model, q[i], 1, dt, hloubka);
    fHC_Copy(vTbck, model.vT); //zaloha teplotniho pole
  end;
  fHC_Free(model);
end;

procedure Sensitivity(n: integer; L: double; MatID: integer; MatDtb: tMatDB; max_volL:double; dt:double; var del_theta: array of double; hloubka:double);
var
  q1, T1, T_s: array of double;
  model: tHCModel;
  i:integer;
begin
  SetLength(q1,n); SetLength(T1,n); SetLength(T_s,n);
  for i:=0 to (n-1) do begin
    q1[i]:=-10000000;
  end;
  InitModel(model, L, 0, MatID, max_volL);
  SetBoundCond(model);
  fHC_ClearAccuracy(model);
  MainDirectTask(T1, T_s, MatDtb, model, q1, n, dt, hloubka);
  fHC_Free(model);
  del_theta[0]:=(T1[0])/(-q1[0]);
  for i:=1 to (n-1) do del_theta[i]:=(T1[i]-T1[i-1])/(-q1[0]);
end;

procedure TempNull(var model: tHCModel; var vTbck: tFloatArray; n: integer; MatDtb: tMatDB; dt:double; var T_qnull: array of double; hloubka: double);
var
  q1, T_s: array of double;
  i:integer;
begin
  SetLength(q1,n); SetLength(T_s,n);
  for i:=0 to (n-1) do
    q1[i]:=0;
  fHC_Copy(model.vT,vTbck);
  MainDirectTask(T_qnull, T_s, MatDtb, model, q1, n, dt, hloubka);
end;

procedure CoefBeck(n: integer; L: double; MatID: integer; MatDtb: tMatDB; max_volL:double; dt:double; var del_theta: array of double; hloubka:double);
var
  q1, T1, T_s: array of double;
  model: tHCModel;
  i:integer;
begin
  SetLength(q1,n); SetLength(T1,n); SetLength(T_s,n);
  for i:=0 to (n-1) do begin
    q1[i]:=-10000000;
  end;
  InitModel(model, L, 0, MatID, max_volL);
  SetBoundCond(model);
  fHC_ClearAccuracy(model);
  MainDirectTask(T1, T_s, MatDtb, model, q1, n, dt, hloubka);
  fHC_Free(model);
  for i:=0 to (n-1) do del_theta[i]:=(T1[i])/(-q1[0]);
end;


procedure IHCP_Full(n:integer; T: array of double; var q: array of double; del_theta: array of double);
var
  X, inv_X: array of array of double;
  i, j, k:integer;
  pom:double;

begin
  SetLength(X,n,n);   SetLength(inv_X,n,n);
  for i:=0 to (n-1) do begin
    q[i]:=0;
    for j:=0 to (n-1) do begin
      if i<j then X[i,j]:=0;
      if i>=j then X[i,j]:= del_theta[i-j];
      if i=j then inv_X[i,j]:=1
      else inv_X[i,j]:=0;
    end;
  end;

  for i:=0 to (n-1) do
    for j:=0 to i do begin
      if i=j then inv_X[i,j]:=inv_X[i,j]/X[i,i];
      X[i,j]:=X[i,j]/X[i,i];
    end;

  for i:=(n-1) downto 1 do
    for j:=(i-1) downto 0 do begin
      pom:=X[i,j];
      for k:=j downto 0 do begin
        X[i,k]:=X[i,k]-X[j,k]*pom;
        inv_X[i,k]:=inv_X[i,k]-inv_X[j,k]*pom;
      end;
    end;

  for i:=0 to (n-1) do begin
    for j:= 0 to (n-1) do begin
      q[i]:=q[i]+inv_X[i,j]*T[j];
    end;
  end;

end;

procedure IHCP_Full_reg(n:integer; T: array of double; var q: array of double; del_theta: array of double);
var
  X, Xt, M, inv_M, F: array of array of double;
  i, j, k:integer;
  alpha, pom:double;

begin
  alpha:=del_theta[0]*del_theta[0];
  SetLength(X,n,n); SetLength(Xt,n,n);
  SetLength(M,n,n); SetLength(inv_M,n,n); SetLength(F,n,n);
  for i:=0 to (n-1) do begin
    q[i]:=0;
    for j:=0 to (n-1) do begin
      if i<j then X[i,j]:=0;                       //naplneni matice X
      if i>=j then X[i,j]:= del_theta[i-j];
      Xt[j,i]:=X[i,j];                             //naplneni matice Xt
      M[i,j]:=0; F[i,j]:=0;
      if i=j then inv_M[i,j]:=1                    //naplneni inverzni matice M (zatim jednotkova)
      else inv_M[i,j]:=0;
    end;
  end;

  for i:=0 to (n-1) do
    for j:=0 to (n-1) do begin
      for k:=0 to (n-1) do
        M[i,j]:= M[i,j]+Xt[i,k]*X[k,j];
      if i=j then M[i,j]:=M[i,j]+alpha;
    end;

  for j:=0 to (n-1) do
    for i:=0 to (n-1) do
      if i<>j then begin
        pom:=M[i,j]/M[j,j];
        for k:=0 to (n-1) do begin
          M[i,k]:=M[i,k]-M[j,k]*pom;
          inv_M[i,k]:=inv_M[i,k]-inv_M[j,k]*pom;
        end;
      end;

  for i:=0 to (n-1) do begin
    for j:=0 to (n-1) do begin
      inv_M[i,j]:=inv_M[i,j]/M[i,i];
      M[i,j]:=M[i,j]/M[i,i];
    end;
  end;

  for i:=0 to (n-1) do
    for j:=0 to (n-1) do begin
      if j>i then inv_M[i,j]:=inv_M[j,i];
    end;

  for i:=0 to (n-1) do begin
    for j:=0 to (n-1) do begin
      for k:=0 to (n-1) do
        F[i,j]:= F[i,j]+inv_M[i,k]*Xt[k,j];
      q[i]:=q[i]+F[i,j]*T[j];
    end;
  end;

end;
end.

