{ Калибровка шкалы масс }

{---------------------------------------------------------------------------
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-47-15                       тел. 75-47-15

   E-mail: aleks@dpt.ustu.ru
----------------------------------------------------------------------------}
unit MassClbr;
{$X+}

interface
   Uses MItypes;

type
  tCounter=longint;
  tSqrCounter=Comp;

  tNodeNumber=1..3;

  tMassBaseCalibration=record
    { основные данные калибровки M=K*(C+С0)^2+M0}
    dC,M0,K:tMass;
  end;

  tPrivateData=record
    BaseCalibration:tMassBaseCalibration;
//    CounterShift:tCounter;
  end;

  tMassValid=(vMassOK, vMassTooSmall, vMassTooLarge);

  tPMassCalibration=^tMassCalibration;
  tMassCalibration=object
    procedure Init;
    procedure Validate;

    function Valid:boolean;
    function Equal(AMassCalibration:tMassCalibration):boolean;

//    function  BaseCalibrationGetCoeff:tMass;
//    function  BaseCalibrationGetOrigin:tMass;
//    function  BaseCalibrationGetCounterOrigin:tMass;

    function  BaseCalibrationCalculate(M1,M2:tMass; C1,C2:tCounter):boolean;
    function  BaseCalibrationCalculateM(NewM1,NewM2,OldM1,OldM2:tMass):boolean;

    function  Get_K:tMass;
    function  Get_dC:tMass;
    function  Get_M0:tMass;

    procedure Set_K(K:tMass);
    procedure Set_M0(M0:tMass);
    procedure Set_dC(dC:tMass);


//    function  CalculateAndSet(M1,M2:tMass; C1,C2:tCounter):boolean;
    procedure CalculateAndSet_dC(C:tCounter; M:tMass);
    function  Calculate_dC(C:tCounter; M:tMass):tMass;
    function  CalculateAndSetM(OldM1,NewM1, OldM2,NewM2:tMass):boolean;
    function  Calculate(OldM1,NewM1, OldM2,NewM2:tMass;
                        var dC:tMass; var K:tMass):boolean;

//    {  Вычисление новой калибровки }
//    function  CalculateNew(M1,M2,         {текущие значения масс}
//                           M01,M02:tMass; {точные значения масс}):boolean;

    function  MassHigh:tMass;
    function  MassLow:tMass;
    function  Mass2Counter(m:tMass):tCounter;
    function  ValidMass(m:tMass):boolean;
    function  ValidMassX(m:tMass):tMassValid;
    function  Counter2Mass(c:tCounter):tMass;

    function  Convert(m:tMass; AOldCalibration:tMassCalibration):tMass;

{    function Refine(M:tMass; Cnt:tCounter):boolean;
    procedure SetCalibration(M:tMass; Cnt:tCounter);
    function GetMass(i:tNodeNumber):tMass;
    function GetCount(i:tNodeNumber):longint;}
  private
    prData:tPrivateData;
    function Mass2CounterX(m:tMass):extended;
  end;



implementation
{$IfNDef Seg16}
    Uses sysutils;
{$EndIf}

function tMassCalibration.Valid:boolean;
begin
  Valid:=prData.BaseCalibration.K>1.0E-30;
end;

function tMassCalibration.Equal(AMassCalibration:tMassCalibration):boolean;
begin
  Equal:=(prData.BaseCalibration.K=AMassCalibration.prData.BaseCalibration.K)
     and (prData.BaseCalibration.M0=AMassCalibration.prData.BaseCalibration.M0)
     and (prData.BaseCalibration.dC=AMassCalibration.prData.BaseCalibration.dC);
end;

function  tMassCalibration.Convert(m:tMass; AOldCalibration:tMassCalibration):tMass;
begin
  Convert:=Counter2Mass(AOldCalibration.Mass2Counter(m));
end;

type
  tCompareResult=(crInvalid, crEqual, crBelow, crAbove);

{function Compare(const x1,x2:tMassCounterRelay):tCompareResult;
begin
  if (x1.Mass<x2.Mass) and (x1.Counter<x2.Counter) then begin
    Compare:=crBelow;
  end else if  (x1.Mass>x2.Mass) and (x1.Counter>x2.Counter)  then begin
    Compare:=crAbove;
  end else if  (x1.Mass=x2.Mass) or (x1.Counter=x2.Counter)  then begin
    Compare:=crEqual;
  end else begin
    Compare:=crInvalid;
  end;
end;}

procedure tMassCalibration.Init;
begin
  Set_K(1);
  Set_M0(0);
  Set_dC(0);
end;

procedure tMassCalibration.Validate;
begin
  if prData.BaseCalibration.K<=0 then
    prData.BaseCalibration.K:=1;
end;

procedure tMassCalibration.Set_K;
begin
  if K<0 then K:=ABS(K);
  if K=0 then K:=1.0E-30;
  prData.BaseCalibration.K:=K;
end;

procedure tMassCalibration.Set_M0;
begin
  prData.BaseCalibration.M0:=M0;
end;

procedure tMassCalibration.Set_dC;
begin
  prData.BaseCalibration.dC:=dC;
end;

procedure tMassCalibration.CalculateAndSet_dC(C:tCounter; M:tMass);
begin
  Set_dC(Calculate_dC(C,M));
end;

function tMassCalibration.Calculate_dC(C:tCounter; M:tMass):tMass;
  var K:tMass;
begin
  K:=Get_K;
  if K>0 then begin
    Calculate_dC:=SQRT(M/K)-C;
  end else begin
    Calculate_dC:=0;
  end;
end;

function  tMassCalibration.Get_K;
begin
  if prData.BaseCalibration.K<=0 then
{$IfDef Seg16}
    RunError(201);
{$Else IfDef Seg16}
    raise ERangeError.Create('tMassCalibration.GetK: Неверное значение К.');
{$EndIf Def Seg16}
  Get_K:=prData.BaseCalibration.K;
end;

function  tMassCalibration.Get_M0;
begin
  Get_M0:=prData.BaseCalibration.M0;
end;

function  tMassCalibration.Get_dC;
begin
  Get_dC:=prData.BaseCalibration.dC;
end;

function tMassCalibration.Calculate(OldM1,NewM1,OldM2,NewM2:tMass;
                                                   var dC:tMass;
                                                   var K:tMass):boolean;
  var
    _OldC1,_c1:extended;
    _dC:extended;
begin
  if (NewM2=NewM1) or (OldM2=OldM1)
    or ((NewM1<NewM2) and (OldM1>OldM2)) or ((NewM1>NewM2) and (OldM1<OldM2))
    or (NewM1<1) or (NewM2<1) then begin
    Calculate:=FALSE;
  end else begin
    _OldC1:=Mass2CounterX(OldM1);
    _dC:=Mass2CounterX(OldM2)-_OldC1;
    if Abs(_dC)<1 then begin
      Calculate:=FALSE;
    end else begin
      K:=SQR( (SQRT(NewM2)-SQRT(NewM1))/_dC );
      _c1:=SQRT(NewM1/K);
      dC:=_c1-_OldC1;
      Calculate:=TRUE;
    end;
  end;
end;

function tMassCalibration.CalculateAndSetM(OldM1,NewM1, OldM2,NewM2:tMass):boolean;
  var dC:tMass; K:tMass;
begin
  if Calculate(OldM1,NewM1, OldM2,NewM2, dC,K) then begin
    CalculateAndSetM:=TRUE;
    Set_K(K);
    Set_M0(0);
    Set_dC(dC);
  end else begin
    CalculateAndSetM:=FALSE;
  end;
end;

function  tMassCalibration.BaseCalibrationCalculate(M1,M2:tMass; C1,C2:tCounter):boolean;
  var
    dm:tMass;
    cc1,cc2, dcc:tMass;
    K,M0:tMass;
begin
  if (C1<0) or (C2<0) or (C2=C1) or (M2=M1) then begin
    BaseCalibrationCalculate:=FALSE;
  end else begin
    cc1:=C1; cc1:=cc1*C1;
    cc2:=C2; cc2:=cc2*C2;
    dcc:=cc2-cc1;
    dm:=M2-M1;

    K:=dm/dcc;
    M0:=M1-K*cc1;
    Set_M0(M0);
    Set_K(K);
    BaseCalibrationCalculate:=TRUE;
  end;
end;

function tMassCalibration.BaseCalibrationCalculateM(NewM1,NewM2,OldM1,OldM2:tMass):boolean;
begin
  BaseCalibrationCalculateM:=BaseCalibrationCalculate(
    NewM1,NewM2,
    Mass2Counter(OldM1),Mass2Counter(OldM2)
  );
end;

(*function tMassCalibration.Refine;

  function TryCalc:boolean;
  begin
    If Valid then begin
      CalculateCalibration;
      TryCalc:=True;
      Refine:=True;
    end else begin
      TryCalc:=False;
    end;
  end;

  var i:integer; mcr,tmp:tMassCounterRelay;
begin
  mcr.Mass:=M;
  mcr.Counter:=Cnt*Cnt;
  i:=1;
  Refine:=False;
  while (i<=prNC) and (Compare(prMC[i],mcr)=crBelow)  do Inc(i);
  if (i>prNC) then begin
    tmp:=prMC[prNC];
    prMC[prNC]:=mcr;
    If not TryCalc then begin
      prMC[prNC]:=tmp;
    end;
  end else begin
    case Compare(prMC[i],mcr) of
      crAbove: begin
        if prNC<High(prNC) then begin
          Inc(prNC);
          prMC[prNC]:=mcr;
          If not TryCalc then begin
            Dec(prNC);
          end;
        end else begin
          if (i>1) and ((prMC[i].Mass-mcr.Mass)>(mcr.Mass-prMC[i-1].Mass)) then Dec(i);
          tmp:=prMC[i];
          prMC[i]:=mcr;
          If not TryCalc then begin
            prMC[i]:=tmp;
          end;
        end;
      end;

      crEqual: begin
        tmp:=prMC[i];
        prMC[i]:=mcr;
        If not TryCalc then begin
          prMC[i]:=tmp;
        end;
      end;

      else

    end;
  end;
  CalculateCalibration;
end;

function tMassCalibration.GetMass;
begin
  if i<=prNC then GetMass:=prMC[i].Mass else GetMass:=-1;
end;

function tMassCalibration.GetCount;
begin
  if i<=prNC then GetCount:=Round(SQRT(prMC[i].Counter)) else GetCount:=-1;
end;

procedure tMassCalibration.SetCalibration;
begin
  prMC[1].Mass:=M;
  prMC[1].Counter:=Cnt*Cnt;
  prNC:=1;
  CalculateCalibration;
end;
*)

function tMassCalibration.Mass2CounterX;
  var x:extended;
begin
  x:=(m-prData.BaseCalibration.M0)/prData.BaseCalibration.K;
  if x<0 then x:=0;
  Mass2CounterX:=SQRT(x)-prData.BaseCalibration.dC;
end;

function tMassCalibration.Mass2Counter;
begin
  if ValidMass(m) then begin
//  Mass2Counter:=Round(SQRT((m-prData.BaseCalibration.M0)/prData.BaseCalibration.K));
    Mass2Counter:=Round(Mass2CounterX(m));
  end else begin
{$IfDef Seg16}
    RunError(201);
{$Else IfDef Seg16}
    raise ERangeError.Create('tMassCalibration.Mass2Counter: Неверное значение массы.');
{$EndIf Def Seg16}
  end;
end;

function  tMassCalibration.ValidMass(m:tMass):boolean;
begin
  ValidMass:=ValidMassX(m)=vMassOK;
end;

function tMassCalibration.ValidMassX;
  var x:extended;
begin
  x:=Mass2CounterX(m);
  if x<Low(tCounter) then
    ValidMassX:=vMassTooSmall
  else if x>High(tCounter) then
    ValidMassX:=vMassTooLarge
  else
    ValidMassX:=vMassOK;
end;

(*function tMassCalibration.Mass2Counter;
  var x:double;
begin
  Mass2Counter:=0;
  if m>0 then begin
    if mA>(1e-300) then begin
      x:=mB*mB+4*mA*(m-mC);
      Mass2Counter:=Round(SQRT((0.5/mA)*(SQRT(x)-mB)));
    end else if mB>(1e-300) then begin
      Mass2Counter:=Round(SQRT((m-mC)/mB));
    end;
  end;
end;*)

function tMassCalibration.Counter2Mass;
begin
  with prData.BaseCalibration do begin
//    Counter2Mass:=M0+K*c*c;
    Counter2Mass:=M0+K*SQR(c+dC);
  end;
end;

function  tMassCalibration.MassHigh:tMass;
begin
  MassHigh:=Counter2Mass(High(tCounter));
end;

function  tMassCalibration.MassLow:tMass;
begin
  MassLow:=Counter2Mass(Low(tCounter));
end;

(*function tMassCalibration.Counter2Mass;
  var rc:tSqrCounter;
begin
  rc:=c;
  rc:=rc*rc;
  Counter2Mass:=(mA*rc*rc+mB*rc+mC);
end;*)

function C123(c1,c2,c3:tSqrCounter):double;
begin
  C123:=c1*c2*c2-c1*c3*c3+c2*c3*c3-c2*c1*c1+c3*c1*c1-c3*c2*c2;
end;

function CM123(c1,c2,c3:tSqrCounter; m1,m2,m3:tMass):double;
begin
  CM123:=-(c3*m2+c2*m1+c1*m3)+(c3*m1+c2*m3+c1*m2);
end;

function CCM123(c1,c2,c3:tSqrCounter; m1,m2,m3:tMass):double;
begin
  c1:=c1*c1;
  c2:=c2*c2;
  c3:=c3*c3;
  CCM123:=-(c1*m2+c2*m3+c3*m1)+(c1*m3+c2*m1+c3*m2);
end;

function CCCM123(c1,c2,c3:tSqrCounter; m1,m2,m3:tMass):double;
  var cc1,cc2,cc3:double;
begin
  cc1:=c1*c1;
  cc2:=c2*c2;
  cc3:=c3*c3;
  CCCM123:=c1*cc2*m3+c2*cc3*m1+c3*cc1*m2-(c3*cc2*m1+c2*cc1*m3+c1*cc3*m2);
end;

(*procedure tMassCalibration.CalculateCalibration;
  var dc:double;
begin
  case prNC of
    3: begin
      dc:= C123(prMC[1].Counter, prMC[2].Counter, prMC[3].Counter);
      mA:=CM123(prMC[1].Counter, prMC[2].Counter, prMC[3].Counter,
                prMC[1].Mass, prMC[2].Mass, prMC[3].Mass)/dc;
      mB:=CCM123(prMC[1].Counter, prMC[2].Counter, prMC[3].Counter,
                prMC[1].Mass, prMC[2].Mass, prMC[3].Mass)/dc;
      mC:=CCCM123(prMC[1].Counter, prMC[2].Counter, prMC[3].Counter,
                prMC[1].Mass, prMC[2].Mass, prMC[3].Mass)/dc;
    end;
    2: begin
      dc:=(prMC[2].Counter-prMC[1].Counter);
      mA:=0;
      mB:=(prMC[2].Mass-prMC[1].Mass)/dc;
      mC:=(prMC[2].Counter*prMC[1].Mass-prMC[1].Counter*prMC[2].Mass)/dc;
    end;

    else begin
      mA:=0;
      mB:=prMC[1].Mass/prMC[1].Counter;
      mC:=0;
    end;
  end;
end;

function tMassCalibration.Valid:boolean;
begin
  Valid:=(prNC<3) or
         (C123(prMC[1].Counter, prMC[2].Counter, prMC[3].Counter)>1e-300);
end;
*)
end.