Ошибка External: SIGFPE - Pascal

program katinaradost;
uses math;
var z, pk, tk, w, ps,k,fi,ac,bc,acc,alfaT,x,y,fl,fv: array [1..3] of real;
  pp,t,r,am,bm,fipr,fil,fiv,summa,a,b,cc,m,psi,v,zl,zv: real;
  cb: array [1..3,1..3] of real;
  i,j: integer;
procedure kardano(a,b,cc,usl:real; var yy:real);
var alfa,beta,p,q,qb: real; y: array [1..3] of real;
begin
 p:=-a*a/3+b;
 q:=2*(a/3)*(a/3)*(a/3)-a*b/3+cc;
 qb:=(p/3)*(p/3)*(p/3)+(q/2)*(q/2);
 
 if (qb<0) and (p<0) then
  begin
 
   alfa:=arccos(-q/(2*sqrt((-p/3)*(-p/3)*(-p/3))));
   y[1]:=2*sqrt(-p/3)*cos(alfa/3);
   y[2]:=-2*sqrt(-p/3)*cos(alfa/3+3.14/3);
   y[3]:=-2*sqrt(-p/3)*cos(alfa/3-3.14/3);
   if usl=1 then {наибольшее}
    begin
     if y[1]>y[2] then yy:=y[1] else yy:=y[2];
     if yy<y[3] then yy:=y[3];
    end;
   if usl=2 then {наименьшее}
    begin
     if y[1]>y[2] then yy:=y[2] else yy:=y[1];
     if yy>y[3] then yy:=y[3];
    end;
   end;
  if (qb>=0) and (p>0) then
   begin
    beta:=(2/q*sqrt((p/3)*(p/3)*(p/3)));
    alfa:=arctan(exp(1/3*ln(beta)));
    yy:=-2*sqrt(p/3)*cos(2*alfa)/sin(2*alfa);
   end;
  if (qb>=0) and (p<0) then
   begin
   beta:=arcsin(2/q*sqrt((-p/3)*(-p/3)*(-p/3)));
   alfa:=arctan(exp(1/3*ln(sin(beta)/cos(beta))));
   yy:=-2*sqrt(-p/3)*1/sin(2*alfa);
   end;
end;
 
begin
  pk[1]:= 54.4*1e5;  {метан}
  tk[1]:= 190.6;
  w[1]:= 0.008;
  z[1]:= 0.6;
  pk[2]:= 48.2*1e5; {этан}
  tk[2]:= 305.4;
  w[2]:= 0.098;
  z[2]:= 0.15;
  pk[3]:= 37.5*1e5;  {бутан}
  tk[3]:= 425.2;
  w[3]:= 0.193;
  z[3]:= 0.25;
  pp:=2*1e7;
  t:=343;
  r:=8.31;
 
  cb[2,1]:=1;
  cb[1,2]:=cb[2,1];
  cb[1,3]:=3.1;
  cb[3,1]:= cb[1,3];
  cb[2,3]:=0.45;
  cb[3,2]:=cb[2,3];
  cb[1,1]:=0;
  cb[2,2]:=0;
  cb[3,3]:=0;

  for i:=1 to 3 do
   begin
    ps[i]:=exp(5.373*(1+w[i])*(1-tk[i]/t))*pk[i];
    k[i]:=ps[i]/pp;
   end;
 
  repeat
  a:=(k[2]*k[3]*z[1]+k[1]*k[3]*z[2]+k[1]*k[2]*z[3]-k[2]*z[1]-k[3]*z[2]-k[1]*z[3]-k[3]*z[1]-k[1]*z[2]-k[2]*z[3]+z[1]+z[2]+z[3]-k[1]*k[2]-k[1]*k[3]-k[2]*k[3]+2*(k[1]+k[2]+k[3]-1))/(-k[1]*k[2]*k[3]+k[1]*k[2]+k[1]*k[3]+k[2]*k[3]-k[1]-k[2]-k[3]+1);
  b:=(k[3]*z[1]+k[1]*z[2]+k[2]*z[3]+k[2]*z[1]+k[3]*z[2]+k[1]*z[3]-2*(z[1]+z[2]+z[3]-1)-k[1]-k[2]-k[3])/(-k[1]*k[2]*k[3]+k[1]*k[2]+k[1]*k[3]+k[2]*k[3]-k[1]-k[2]-k[3]+1);
  cc:=(z[1]+z[2]+z[3]-1)/(-k[1]*k[2]*k[3]+k[1]*k[2]+k[1]*k[3]+k[2]*k[3]-k[1]-k[2]-k[3]+1);
 
  kardano(a,b,cc,0,v); {мольная доля пар. фазы V}
 
  for i:=1 to 3 do
   begin
    x[i]:=z[i]/(1-v+k[i]*v);
    y[i]:=z[i]/(1-v+k[i]*v)*k[i];
   end;
 
  for i:=1 to 3 do
   begin
    acc[i]:=0.46*r*r*tk[i]*tk[i]/pk[i];
    if w[i]<0.49 then m:=0.37+1.54*w[i]-0.27*w[i]*w[i] else m:=0.38+1.41*w[i]-0.16*w[i]*w[i]+0.017*w[i]*w[i]*w[i];
    alfaT[i]:=(1+m*(1-sqrt(t/tk[i])))*(1+m*(1-sqrt(t/tk[i])));
    ac[i]:=acc[i]*alfaT[i];
    bc[i]:=0.078*r*tk[i]/pk[i];
   end;

  for i:=1 to 3 do
   for j:=1 to 3 do
    begin
     am:=am+(1-cb[i,j])*x[i]*x[j]*sqrt(ac[i]*ac[j]);
    end;
  for i:=1 to 3 do
   begin
    bm:=bm+x[i]*bc[i];
   end;
 
  a:=bm-r*t/pp;
  b:=(r*t*(am-2*b))/pp-3*bm*bm;
  cc:=bm*bm*bm+r*t*(bm*bm-am*bm)/pp;
 
  kardano(a,b,cc,2,zl); {коэф. сферхсжим жля жидк. фазы, из ур. PR}
 
  for i:=1 to 3 do
   for j:=1 to 3 do
    begin
     am:=am+(1-cb[i,j])*y[i]*y[j]*sqrt(ac[i]*ac[j]);
    end;
  for i:=1 to 3 do
   begin
    bm:=bm+y[i]*bc[i];
   end;

  a:=bm-r*t/pp;
  b:=(r*t*(am-2*b))/pp-3*bm*bm;
  cc:=bm*bm*bm+r*t*(bm*bm-am*bm)/pp;
 
  kardano(a,b,cc,1,zv); {коэф. сферхсжим жля пар. фазы, из ур. PR}

  for i:=1 to 3 do
   begin
  psi:=0.37+1.57*w[i]-0.27*sqr(w[i]);
  fipr:=sqr(1+psi*(1-sqrt(t/tk[i])));
  fil:=exp(zl-1-ln(zl-0.078*(pp/pk[i])/(t/tk[i]))-2.078*fipr/(t/tk[i])*ln((zl+0.19*(pp/pk[i])/(t/tk[i]))/(zl-0.032*(pp/pk[i])/(t/tk[i]))));
  fiv:=exp(zv-1-ln(zv-0.078*(pp/pk[i])/(t/tk[i]))-2.078*fipr/(t/tk[i])*ln((zv+0.19*(pp/pk[i])/(t/tk[i]))/(zv-0.032*(pp/pk[i])/(t/tk[i]))));
  fl[i]:=fil*pp*z[i];
  fv[i]:=fiv*pp*z[i];
   end;

  summa:=0;
 
  for i:=1 to 3 do
   begin
  k[i]:=fl[i]/fv[i];
  summa:=summa+(1-k[i])*(1-k[i]);
   end;
 
  until summa<1e-8;
 
  writeln(x[1],x[2],x[3]);
  writeln(y[1],y[2],y[3]);

end.