{ Определение аппроксимирующей функции (F) и ее частных производных (GradF)
для процедуры обобщенного поиска минимума функции многих переменных 
(см. модуль uGenFit). 

Аппроксимирующая функция: 
  F(x):=A*(erf(s*((x-x0)+h))-erf(s*((x-x0)-h)))+)+B
Параметры и их порядок в структуре данных:
  P=(A, B, x0, h, s);
}

{---------------------------------------------------------------------------
The control units for mass-spectrometer MI1201-AGM
(c) Copyright	 Aleksandrov O.E., 1998

Модуль управления масс-спектрометром МИ1201-АГМ
(c) Собственность     Александрова О.Е., 1998

   Molecular Physics department	 	620002, Екатеринбург, К-2
   USTU, Ekaterinsburg, K-2, 620002     УГТУ,
   RUSSIA     	   			Кафедра молекулярной физики
   phone 75-48-39                       тел. 75-48-39

   E-mail: aleks@dpt.ustu.ru
----------------------------------------------------------------------------}

unit uFuncDef;

interface
   USES uGenFit;

type
  tFuncParams=record
    Size:tParameterNumber;
    A,B,x0,h,s:Double;
  end;
  tPtrFuncParams=^tFuncParams;

const
  cParSize=SizeOf(tFuncParams) div SizeOf(Double);

function F(const P:tParameters; X:Double):Double;
procedure GradF(const P:tParameters; X:Double; var Grad:tParameters);

procedure SelectInitials(var P:tParameters; const X,Y:tData);

implementation
  USES uErf;

function F;
begin
  with tPtrFuncParams(@P)^ do begin
{$IfOpt R+}
    If Size<>cParSize then RunError(201);
{$EndIf}
    x:=x-x0;
    F:=A*(Erf(s*(x+h))-Erf(s*(x-h)))+B;
  end;
end;

procedure GradF;
  const
    c2divSqrtPI=2/1.7724538509055160273; { = 2/Sqrt(PI) }
  var
    x1,x2,sx1,sx2,expSx1,expSx2:Extended;
    A_2divSqrtPI,AS_2divSqrtPI:Extended;
begin
  with tPtrFuncParams(@P)^ do begin
{$IfOpt R+}
    If Size<>cParSize then RunError(201);
{$EndIf}
    x:=x-x0;
    x1:=(x+h);
    x2:=(x-h);
    sx1:=s*x1;
    sx2:=s*x2;
    expSx1:=exp(-sx1*sx1);
    expSx2:=exp(-sx2*sx2);
    A_2divSqrtPI:=A*c2divSqrtPI;
    AS_2divSqrtPI:=S*A_2divSqrtPI;
  end;
  with tPtrFuncParams(@Grad)^ do begin
    A :=Erf(sx1)-Erf(sx2); 		    { dF/dA }
    B :=1; 				    { dF/dB }
    X0:=AS_2divSqrtPI*(-expSx1+expSx2);     { dF/dX0 }
    H :=AS_2divSqrtPI*(expSx1+expSx2);      { dF/dH }
    S :=A_2divSqrtPI*(x1*expSx1-x2*expSx2); { dF/dS }
  end;
end;

{ Выбор начального приближения для параметров }
procedure SelectInitials;
  var
    dX,maxY,minY,minX,maxX,iX0:Double;
    i:word;
begin
  with tPtrFuncParams(@P)^ do begin
{$IfOpt R+}
    If Size<>cParSize then RunError(201);
    If (X.Size<>Y.Size) then RunError(201);
    If (X.Size<=Size+2) then RunError(201);
    If (X.Size>cMaxDataSize) then RunError(201);
{$EndIf}
    maxY:=Y.Data[1];
    minY:=Y.Data[1];
    minX:=X.Data[1];
    maxX:=X.Data[1];
    iX0:=minX;
    for i:=2 to X.Size do begin
      if maxY<Y.Data[i] then begin
	maxY:=Y.Data[i];
	iX0:=X.Data[i];
      end;
      if Y.Data[i]<minY then minY:=Y.Data[i];
      if maxX<X.Data[i] then maxX:=X.Data[i];
      if X.Data[i]<minX then minX:=X.Data[i];
    end;
{$IfOpt R+}
    If maxX=minX then RunError(201);
    If maxY=minY then RunError(201);
{$EndIf}
    A:=maxY;
    B:=minY;
    X0:=iX0;
    dX:=(maxX-minX);
    h:=(0.05/2)*dX;
    s:=(1/(Sqrt(2)*0.25))/dX;
  end;
end;

end.