unit Appr;
{!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		  АППРОКСИМАЦИЯ ФУНКЦИЙ



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!}
interface

uses WinTypes, WinProcs,gltype,glproc,MathObj,Math, Classes, SysUtils, MVTU_TLB;

function TIntPol(at,adt:RealType;var time,dtime:RealType;var AU,AY:TPtrExt;var AX,ADX:array of RealType;
	       var CU,CY : TIntArray;var Prop:TPtrArray;var Vars:Pointer;Action:Integer):Integer;export;
function TMNK(at,adt:RealType;var time,dtime:RealType;var AU,AY:TPtrExt;var AX,ADX:array of RealType;
	       var CU,CY : TIntArray;var Prop:TPtrArray;var Vars:Pointer;Action:Integer):Integer;export;
function TTrap(at,adt:RealType;var time,dtime:RealType;var AU,AY:TPtrExt;var AX,ADX:array of RealType;
	       var CU,CY : TIntArray;var Prop:TPtrArray;var Vars:Pointer;Action:Integer):Integer;export;


{********************************************************************************************************}

                               IMPLEMENTATION

{********************************************************************************************************}
uses LUD,Interpol;

type
PIntPolRec = ^TIntPolRec;
TIntPolRec = record
  Met,
  Order,
  N,M,Nfun      : ^Integer;
 end;

PApprRec = ^TApprRec;
TApprRec = record
  fmType,
  Order,
  N,Nfun  : ^Integer;
 end;

PTrapRec = ^TTrapRec;
TTrapRec = record
  N,Nfun  : ^Integer;
end;

PLAEVar = ^TLAEVar;
TLAEVar = record
 A,B,
 Arg : TExtArray;
 Idn : TIntArray;
end;

PSplineVar = ^TSplineVar;
TSplineVar = record
 Y2,Tmp : Array of RealType;
end;
//**********************************************************************

//     Функции формы для МНК

//**********************************************************************
function TIntPolConvert(Action:  integer;
                          Vars:    Pointer;
                          Prop:    TPtrArray;
                          PropStr: TStringList;
                          BlockId: integer;
                          var MVTU: IMVTU_Server):integer;
  const RecName:PChar = 'Интерполяция';
  var
      Res: integer;
begin
  case Action of
    //Возвращаем ссылку на имя записи в базе МВТУ-4
    cnv_GetRecName: Result:=integer(RecName);
    //Конвертируем и устанавливаем параметры при помощи интерфейса МВТУ
    cnv_Convert: with PIntpolRec(Prop.arr)^ do begin
                   MVTU.SetBlockProp(BlockId,'met',PChar(IntToStr(met^)),Res);
                   MVTU.SetBlockProp(BlockId,'order',PChar(PropStr[1]),Res);
                   MVTU.SetBlockProp(BlockId,'n',PChar(PropStr[2]),Res);
                   MVTU.SetBlockProp(BlockId,'m',PChar(PropStr[3]),Res);
                   MVTU.SetBlockProp(BlockId,'nfun',PChar(PropStr[4]),Res);
                 end;
  end;
end;

function MNKForm(What,m:Integer;x:RealType):RealType;
//
//    What - тип функции формы
//    0 - полином a0+a1*x+...+am*x^m
//    1 - a0*sin(0*pi*x)+ a1*sin(1*pi*x)+...+am*sin(m*pi*x)
//    2 - x^m*(1-x)
//
//    X - аргумент
//    M - текущий порядок функции формы
//
begin
   case What of
    0: Result:=IntPower(x,m);
    1: Result:=sin((m+1)*pi*x);
    2: if m = 0 then Result:=1 else Result:=IntPower(x,m)*(1-x);
    else Result:=0;
   end;
end;
//**********************************************************************

//     Вычисление функции формы в точке для МНК

//**********************************************************************
function MNKFormEval(What,m:Integer;x:RealType;var a : array of RealType):RealType;
//
//    What - тип функции формы
//    X - аргумент
//    M - текущий порядок функции формы
//    a - массив значений коэффициентов функции формы
//
var i : Integer;
begin
  Result:=0;
  for i:=0 to m do
   case What of
    0: Result:=Result+a[i]*IntPower(x,i);
    1: Result:=Result+a[i]*sin((i+1)*pi*x);
    2: if i = 0 then Result:=a[i] else Result:=Result+a[i]*IntPower(x,i)*(1-x);
   end;
end;
//**********************************************************************

//     Коэффициеты матрицы для МНК

//**********************************************************************
function MNKMatrix(What,m,n:Integer;x:RealType):RealType;
//
//    What - тип функции формы
//    X - аргумент
//    M - текущий порядок функции формы
//    N - текущий порядок весовой функции
//
begin
 Result:=MNKForm(What,m,x)*MNKForm(What,n,x);
end;
//**********************************************************************

//     Коэффициенты вектор-столбца правых частей для МНК

//**********************************************************************
function MNKRight(What,m:Integer;x,y:RealType):RealType;
//
//    What - тип функции формы
//    X - аргумент
//    Y - значение функции
//    M - текущий порядок функции формы
//
begin
 Result:=y*MNKForm(What,m,x);
end;
{---------------------------------------------------------------------------

---------------------------------------------------------------------------}
function TIntPol;
var i,j   : Integer;
    px,py : PExtArr;
    pxx   : TExtArray;
begin
 Result:=0;
 with PIntpolRec(Prop.arr)^ do
 case Action of
 f_GetConvertFuncAdr: Result:=integer(@TIntPolConvert);
 f_GetDeriCount,
 f_GetStateCount,
 f_GetInit    :  Result:=0;
 f_GetCount   :  begin
                  CU.arr^[0]:=N^;
                  CU.arr^[1]:=N^*Nfun^;
                  CY.arr^[0]:=CU.arr^[2]*Nfun^;
                 end;
  f_Create   : begin
                GetMem(Vars,SizeOf(TSplineVar));
                with PSplineVar(Vars)^ do begin
                 SetLength(Y2,0);
                 SetLength(Tmp,0);
                end
               end;
  f_Free     : begin
                with PSplineVar(Vars)^ do begin
                 SetLength(Y2,0);
                 SetLength(Tmp,0);
                end;
                FreeMem(Vars,SizeOf(TLAEVar));
               end;
  f_InitObjects:
               with PSplineVar(Vars)^ do begin
                SetLength(Y2,N^*SOfR);
                SetLength(Tmp,N^*SOfR);
               end;
 f_RestoreOuts,
 f_UpdateJacoby,
 f_GoodStep,
 f_InitState  :  with PSplineVar(Vars)^ do begin
                  px:=AU.Ptr(0).arr;
                  py:=AU.Ptr(1).arr;
                  pxx:=AU.Ptr(2);
                  for i:=0 to Nfun^-1 do begin
                   if i > 0 then py:=@py^[N^];
                   case Met^ of
                   0: for j:=0 to pxx.Count-1 do
                       AY.Ptr(0).arr^[i*pxx.Count+j]:=Lagrange(px^,py^,pxx.arr^[j],Order^,M^);
                   1: begin
                       spline(px^,py^,y2,Tmp,n^,0,0);
                       for j:=0 to pxx.Count-1 do
                        Result:=splint(px^,py^,y2,n^,pxx.arr^[j],AY.Ptr(0).arr^[i*pxx.Count+j]);
                      end;
                   end
		  end
                 end;
 end
end;
{---------------------------------------------------------------------------

---------------------------------------------------------------------------}
function TMNKConvert(Action:  integer;
                          Vars:    Pointer;
                          Prop:    TPtrArray;
                          PropStr: TStringList;
                          BlockId: integer;
                          var MVTU: IMVTU_Server):integer;
  const RecName:PChar = 'МНК-аппроксимация';
  var
      Res: integer;
begin
  case Action of
    //Возвращаем ссылку на имя записи в базе МВТУ-4
    cnv_GetRecName: Result:=integer(RecName);
    //Конвертируем и устанавливаем параметры при помощи интерфейса МВТУ
    cnv_Convert: with PApprRec(Prop.arr)^ do begin
                   MVTU.SetBlockProp(BlockId,'fmType',PChar(IntToStr(fmType^)),Res);
                   MVTU.SetBlockProp(BlockId,'order',PChar(PropStr[1]),Res);
                   MVTU.SetBlockProp(BlockId,'n',PChar(PropStr[2]),Res);
                   MVTU.SetBlockProp(BlockId,'nfun',PChar(PropStr[3]),Res);
                 end;
  end;
end;

function TMNK;
var i,j,k : Integer;
    sum   : RealType;
    fLU   : Boolean;
begin
 Result:=0;
 with PApprRec(Prop.arr)^ do
 case Action of
 f_GetConvertFuncAdr: Result:=integer(@TMNKConvert);
 f_GetDeriCount,
 f_GetStateCount,
 f_GetInit    :  Result:=0;
 f_GetCount   :  begin
                  CU.arr^[0]:=N^;
                  CU.arr^[1]:=N^*Nfun^;
                  CY.arr^[0]:=Nfun^*CU.arr^[2];
                  CY.arr^[1]:=(Order^+1)*Nfun^;
                 end;
  f_Create   : begin
                GetMem(Vars,SizeOf(TLAEVar));
                with PLAEVar(Vars)^ do begin
                 A:=TExtArray.Create(1);
                 B:=TExtArray.Create(1);
                 Arg:=TExtArray.Create(1);
                 Idn:=TIntArray.Create(1);
                end
               end;
  f_Free     : begin
                with PLAEVar(Vars)^ do begin
                 A.Free;
                 B.Free;
                 Idn.Free;
                 Arg.Free;
                end;
                FreeMem(Vars,SizeOf(TLAEVar));
               end;
  f_InitObjects:
               with PLAEVar(Vars)^ do begin
                A.ChangeCount(Sqr(Order^+1));
                B.ChangeCount(Order^+1);
                Idn.ChangeCount(Order^+1);
                Arg.ChangeCount(N^);
                Arg.FillArray(0);
               end;
 f_RestoreOuts,
 f_UpdateJacoby,
 f_GoodStep,
 f_InitState  :  with PLAEVar(Vars)^ do begin

                  fLU:=false;
                  for j:=0 to N^-1 do begin
                   fLU:=Arg.arr^[j] <> AU.Ptr(0).arr^[j];
                   if fLU then break;
                  end;
                  Move(AU.Ptr(0).arr^,Arg.arr^,N^*SOfR);

                  if fLU then begin
                   for i:=0 to Order^ do
                    for k:=0 to Order^ do begin
                     sum:=0;
                     for j:=0 to N^-1 do
                      sum:=sum+MNKMatrix(fmType^,i,k,Arg.arr^[j]);
                      A.arr^[i*(Order^+1)+k]:=sum;
                    end;

                   Result:=ludcmp(a.arr^,b.arr^,idn.arr^,Order^+1);
                   if Result <> 0 then exit
                  end;

                  for k:=0 to Nfun^-1 do begin

                   for i:=0 to Order^ do begin
                    sum:=0;
                    for j:=0 to N^-1 do
                     sum:=sum+MNKRight(fmType^,i,Arg.arr^[j],AU.Ptr(1).arr^[k*N^+j]);
                    B.arr^[i]:=sum;
                   end;

                   lubksb(a.arr^,b.arr^,idn.arr^,Order^+1);

                   for j:=0 to AU.Ptr(2).Count-1 do
                    AY.Ptr(0).arr^[k*AU.Ptr(2).Count+j]:=MNKFormEval(fmType^,Order^,AU.Ptr(2).arr^[j],B.arr^);

                   for j:=0 to Order^ do AY.Ptr(1).arr^[k*(Order^+1)+j]:=B.arr^[j];
                  end
		 end;
 end
end;

{---------------------------------------------------------------------------

---------------------------------------------------------------------------}
function TTrapConvert(Action:  integer;
                          Vars:    Pointer;
                          Prop:    TPtrArray;
                          PropStr: TStringList;
                          BlockId: integer;
                          var MVTU: IMVTU_Server):integer;
  const RecName:PChar = 'Интеграл от табличной функции';
  var
      Res: integer;
begin
  case Action of
    //Возвращаем ссылку на имя записи в базе МВТУ-4
    cnv_GetRecName: Result:=integer(RecName);
    //Конвертируем и устанавливаем параметры при помощи интерфейса МВТУ
    cnv_Convert: with PIntpolRec(Prop.arr)^ do begin
                   MVTU.SetBlockProp(BlockId,'npoint',PChar(PropStr[0]),Res);
                   MVTU.SetBlockProp(BlockId,'nfun',PChar(PropStr[1]),Res);
                 end;
  end;
end;

function TTrap;
var i,j   : Integer;
    sum   : RealType;
    px,py : PExtArr;
begin
 Result:=0;
 with PTrapRec(Prop.arr)^ do
 case Action of
 f_GetConvertFuncAdr: Result:=integer(@TTrapConvert);
 f_GetDeriCount,
 f_GetStateCount,
 f_GetInit    :  Result:=0;
 f_GetCount   :  begin
                  CU.arr^[0]:=N^;
                  CU.arr^[1]:=N^*Nfun^;
                  CY.arr^[0]:=Nfun^;
                 end;
 f_RestoreOuts,
 f_UpdateJacoby,
 f_GoodStep,
 f_InitState  :  begin
                  px:=AU.Ptr(0).arr;
                  py:=AU.Ptr(1).arr;
                  for i:=0 to Nfun^-1 do begin
                   sum:=0;
                   if i > 0 then py:=@py^[N^];
                   for j:=1 to N^-1 do sum:=sum+(px^[j]-px^[j-1])*(py^[j]+py^[j-1])/2;
                   AY.Ptr(0).arr^[i]:=sum;
                  end
		 end;
 end
end;

end.