Error Function Calculator

Source Code:

program erfo;

{$G+,N+}

uses crt;

const
a1: real = 0.254829592;
a2: real = -0.284496736;
a3: real = 1.421413741;
a4: real = -1.453152027;
a5: real = 1.061405429;
p:  real = 0.3275911;
Si: real = 10e-4;
Se: real = 10e-2;
H:  real = 10e5;
Hw: real = 1.3e5;
C0: real = 10e3;
T0:real = 10e3;
dTemp: real = 10;
Mu: real = 1.26e-6;
{-------------------------}

var
tx: double;
i: integer;
F: Text;
T: integer;
Z: Real;
V0,Alfa,Eta,S0,S,Delta,Sp,Sm,Em,Es,B,B2: real;
Tau,Beta:Real;

procedure calcparams;
begin
V0 := (dTemp/T0);
Alfa := sqrt(Se/Si);
Eta := (1+Alfa)/(Hw*Alfa);
S0 := Mu*Si;
S := S0*C0;
Delta := H-Hw;
Sp := S/(Eta+S);
Sm := S/(Eta-S);
Em := Eta/(Eta-S);
Es := Eta/S;
B := Delta{-Alfa*z};
B2 := B/2;
end;


function erf(x:double):double;
var t:real;
begin
if x<=0 then erf:=0 else begin
   t:=1/(1+p*x);
   erf:=1-t*(a1+t*(a2+t*(a3+t*(a4+t*a5))))*exp(-x*x);
   end;
end;

function Teta(X:Real):Byte;
begin
if x<0 then Teta:=0 else Teta:=1;
end;


function X0(T:Real):Real;
begin
X0:=1-erf(sqrt(S0/T)*(Tau-B2));
end;

function X1(T:Real):Real;
begin
X1:=1-erf(sqrt(S0/T)*(Tau+B2));
end;

function X2(T:Real):Real;
begin
X2:=1-erf(sqrt(S0/T)*(Tau*Es+B2));
end;

function X3(Beta:Real):Real;
begin
X3:=1-erf(sqrt(S*Beta));
end;

function X4(Beta:Real):Real;
begin
X4:=1-erf(sqrt(Beta/S)*Eta);
end;

function Hx(T,Z:Real):double;
begin
Tau:= C0*T;
Beta:=Tau-Delta;
Hx:=V0/Hw*Teta(T)*Sp*Sm*(Exp(S*(Tau+B)*X1(T)-Exp(Eta*(Tau*Es+B)*X2(T))-
V0/(2*Hw)*Teta(Beta)*Sp*(Exp(S*Beta)*X3(Beta)-Exp(Eta*Beta*Es)*X4(Beta))));
end;

begin
clrscr;
  Assign(F, '');
  Rewrite(F);
for T:=1 to 100 do begin
    writeln(F,'значение Hх :',Hx(T,0));
    end;
  Close(F);
end.