procedure MEDFIT(X, Y:array of real; NDATA:integer; var A, B, ABDEV:real);
Label 1,2,3;
var
    J:integer;
    SX,SY,SXY,SXX,DEL,AA,BB,CHISQ,SIGB,B1,F1,B2,F2,ZZ,F:real;
begin
    SX:=0; 
    SY:=0;
    SXY:=0; 
    SXX:=0; 
    For J:=1 To NDATA do
    begin
        XT[J]:=X[J];
        YT[J]:=Y[J];
        SX:=SX + X[J];
        SY:=SY + Y[J];
        SXY:=SXY + X[J] * Y[J];
        SXX:=SXX + Sqr(X[J]);
    end; 
    NDATAT:=NDATA;
    DEL:=NDATA * SXX - Sqr(SX);
    AA:=(SXX * SY - SX * SXY) / DEL;
    BB:=(NDATA * SXY - SX * SY) / DEL;
    CHISQ:=0;
    For J:=1 To NDATA do
        CHISQ:=CHISQ + Sqr(Y[J] - (AA + BB * X[J]));
    SIGB:=Sqrt(CHISQ / DEL);
    B1:=BB;
    F1:=ROFUNC(B1);
    If f1 >= 0 then
        ZZ:=1
    Else
        ZZ:=-1;
    B2:=BB + Abs(3 * SIGB) * ZZ;
    F2:=ROFUNC(B2);
1:  If F1 * F2 > 0  Then
    begin
        BB:=2 * B2 - B1;
        B1:=B2;
        F1:=F2;
        B2:=BB;
        F2:=ROFUNC(B2);
        GoTo 1;
    end;
    SIGB:=0.01 * SIGB;
2:  If Abs(B2 - B1) > SIGB Then
    begin
        BB:=0.5 * (B1 + B2);
        If (BB = B1) Or (BB = B2) Then GoTo 3;
        F:=ROFUNC(BB);
        If F * F1 >= 0  Then
        begin
            F1:=F;
            B1:=BB;
        end
        Else
        begin
            F2:=F;
            B2:=BB;
        end;
        GoTo 2;
    end;
3:  A:=AA;
    B:=BB;
    ABDEV:=ABDEVT / NDATA;
end;