program field; {$n+,e-,m 65200,0,65200} uses crt,dos;

type  Werte=record
      hf,s,y: integer;
      u0,i,j,k,o,w,z,i1,i2,i3,i4: byte;
      dg0,z0,dr,r,r0,dt,t,su,dg,K0,c1,c2,c3,c4: extended;
      c01,c02,c03,c04,c10,c20,c30,c40,c1s,c2s,c3s,c4s: extended;
      U: array[1..7,0..4,0..120] of extended;
      Ug,Uh: array[1..7,1..120] of single; z1: single;
      ul: array[1..4] of extended;
      V: array[1..7] of extended;
      a,q: array[1..4] of extended;
      g0,h0,m0,g,h,m,p,f,ri: array[1..4,1..4] of extended;
      n0,l,n: array[1..4,1..4,1..4] of extended;
      Cs,Tg: array[0..120] of extended;
      x1: char; a0: double; xi,code: integer;
      bn,bs: Boolean
      end;  {variables being saved at interruption}

var   wr: Werte;
      ca: char; x2:string[3];
      fi1: text;             {parameters and number of steps}
      fi6: file of extended; {all values, for mathematical analyses}
      fi2: file of double;   {values all 15deg, for printing}
      fi3: file of single;   {limiting radii, for graphics}
      fi5: file of Werte;    {for saving variables at
                              interruption}
      bo: Boolean;
      ExitSave: Pointer;

{$F+} procedure MyExit; {$F-}
begin
      close(fi1); if wr.bn then begin
        close(fi2); close(fi3);
        if wr.bs then close (fi6)
      end; ClrScr;
      ExitProc:= ExitSave;
end;

procedure kSchleife;
  {"innerst" loop, across the angle steps  k, contains all formulae}
begin with wr do
  begin
          for i:=1 to 4 do
            if (i=1) or (i=4) then
              begin g0[i,i]:=U[i,2,k]; l[i,i,1]:=(U[i,1,k]-U[i,3,k])/2/dr;
              l[i,i,2]:=(U[i,2,k+1]-U[i,2,k-1])/2/dt;
              m0[i,i]:=(U[i,2,k+1]-2*U[i,2,k]+U[i,2,k-1])/sqr(dt);
              g[i,i]:=g0[i,i]+1; m[i,i]:=m0[i,i] end;
          g0[2,2]:=r*r*U[2,2,k]; l[2,2,1]:=2*r*U[2,2,k]+r*r*(U[2,1,k]-U[2,3,k])
          /2/dr; l[2,2,2]:=r*r*(U[2,2,k+1]-U[2,2,k-1])/2/dt;
          m0[2,2]:=r*r*(U[2,2,k+1]-2*U[2,2,k]+U[2,2,k-1])/sqr(dt);
          g0[3,3]:=sqr(r*Cs[k])*U[3,2,k]; l[3,3,1]:=sqr(r*Cs[k])*(2*U[3,2,k]/r
          +(U[3,1,k]-U[3,3,k])/2/dr); l[3,3,2]:=sqr(r*Cs[k])*((U[3,2,k+1]-
          U[3,2,k-1])/2/dt-2*U[3,2,k]*Tg[k]);
          g0[3,4]:=U[5,2,k]*sqr(Cs[k])/r; g0[4,3]:=g0[3,4];
          a[3]:=U[6,2,k]*sqr(Cs[k])/r; a[4]:=U[7,2,k]/r;
          l[3,4,1]:=sqr(Cs[k])*(-U[5,2,k]/r+(U[5,1,k]-U[5,3,k])/2/dr)/r;
          l[3,4,2]:=sqr(Cs[k])*((U[5,2,k+1]-U[5,2,k-1])/2/dt-2*U[5,2,k]*
          Tg[k])/r; l[4,3,1]:=l[3,4,1]; l[4,3,2]:=l[3,4,2];
          p[3,1]:=sqr(Cs[k])*(-U[6,2,k]/r+(U[6,1,k]-U[6,3,k])/2/dr)/r;
          p[3,2]:=sqr(Cs[k])*((U[6,2,k+1]-U[6,2,k-1])/2/dt-2*U[6,2,k]*Tg[k])
          /r; p[4,1]:=-U[7,2,k]/r/r+(U[7,1,k]-U[7,3,k])/2/dr/r;
          p[4,2]:=(U[7,2,k+1]-U[7,2,k-1])/2/dt/r;
          m0[3,4]:=sqr(Cs[k])*((U[5,2,k+1]-2*U[5,2,k]+U[5,2,k-1])/dt/dt
          -2*(U[5,2,k+1]-U[5,2,k-1])/dt*Tg[k]+2*U[5,2,k]*(sqr(Tg[k])-1))/r;
          m0[4,3]:=m0[3,4];
          q[3]:=sqr(Cs[k])*((U[6,2,k+1]-2*U[6,2,k]+U[6,2,k-1])/sqr(dt)
          -2*(U[6,2,k+1]-U[6,2,k-1])/dt*Tg[k]+2*U[6,2,k]*(sqr(Tg[k])-1))/r;
          q[4]:=(U[7,2,k+1]-2*U[7,2,k]+U[7,2,k-1])/sqr(dt)/r;
          g[2,2]:=g0[2,2]+r*r; m[2,2]:=m0[2,2];
          g[3,3]:=g0[3,3]+sqr(r*Cs[k]);
          g[3,4]:=g0[3,4]; g[4,3]:=g0[3,4];
          m[3,4]:=m0[3,4]; m[4,3]:=m0[3,4];
          {Potentials and derivations, here are used quantities avoiding
           the formulae with terms of the proportion of 1 .
           These special quantities are marked with '0'.
           All components are formed from a "storing field" U .}

          if (g[1,1]<-100) or (g[2,2]<-100*r*r) or (g[3,3]<-100*sqr(r*Cs[k]))
          or (g[4,4]<-100) then hf:=y;
          if (abs(g[3,4])>100*r*Cs[k]) or (abs(a[3])>100*r*Cs[k]) or
          (abs(a[4])>100) then hf:=y;
          {for controlled drop out}

          dg0:=g0[3,3]+sqr(r*Cs[k])*g0[4,4]+g0[3,3]*g0[4,4]+sqr(g0[3,4]);
          dg:=dg0+sqr(r*Cs[k]);
          h0[1,1]:=(r*r*dg0+g0[2,2]*dg)/sqr(r*r*Cs[k]);
          h[1,1]:=h0[1,1]+1;
          h0[2,2]:=(dg0+g0[1,1]*dg)/sqr(r*r*Cs[k]);
          h[2,2]:=h0[2,2]+1/r/r;
          h0[3,3]:=(g[1,1]*g[2,2]*g0[4,4]+g0[1,1]*g[2,2]+g0[2,2])
          /sqr(r*r*Cs[k]); h[3,3]:=h0[3,3]+1/sqr(r*Cs[k]);
          h0[4,4]:=(g0[1,1]*g[2,2]*g[3,3]+g0[2,2]*g[3,3]+r*r*g0[3,3])
          /sqr(r*r*Cs[k]); h[4,4]:=h0[4,4]+1;
          h0[3,4]:=-g[1,1]*g[2,2]*g[3,4]/sqr(r*r*Cs[k]);
          h0[4,3]:=h0[3,4]; h[3,4]:=h0[3,4]; h[4,3]:=h[3,4];
          {contravariant components of the metric tensor}

          for i:=1 to 4 do
            for j:=1 to 4 do f[i,j]:=p[i,j]-p[j,i];
          {field tensor}

          su:=0;
          for i:=1 to 4 do
            for j:=1 to 4 do
              for o:=1 to 4 do
                for w:=1 to 4 do begin
                  if (i=4) or (j=4) or (o=4) or (w=4) then s:=-1 else s:=1;
                  {s=-1 for products from imaginary components}
                  su:=su+s*h[i,j]*h[o,w]*f[i,o]*f[j,w] end;

          for i:=1 to 4 do
            for j:=1 to 4 do
              begin ri[i,j]:=g[i,j]*(su/4-3*K0); for o:=1 to 4 do
                for w:=1 to 4 do begin
                  if (i=4) and (j=4) or (o=4) or (w=4) then s:=-1 else s:=1;
                  ri[i,j]:=ri[i,j]-s*h[o,w]*f[i,o]*f[j,w] end
              end;
          {Ricci tensor}

          for i:=1 to 4 do
            for j:=1 to 4 do
              for o:=1 to 4 do
                begin n0[i,j,o]:=0; for w:=1 to 4 do begin
                  if (o=3) and (w=4) and ((i=3) or (j=3)) or (o=4) and (w=3)
                  and ((i=4) or (j=4)) then s:=-1 else s:=1;
                  n0[i,j,o]:=n0[i,j,o]+s*
                  h[o,w]*(l[j,w,i]+l[i,w,j]-l[i,j,w])/2 end;
                n[i,j,o]:=n0[i,j,o]
                end;

          n[2,2,1]:=n0[2,2,1]-r*h[1,1];
          n[1,2,2]:=n0[1,2,2]+r*h[2,2]; n[2,1,2]:=n[1,2,2];
          n[3,3,1]:=n0[3,3,1]-r*sqr(Cs[k])*h[1,1];
          n[1,3,3]:=n0[1,3,3]+r*sqr(Cs[k])*h[3,3]; n[3,1,3]:=n[1,3,3];
          n[3,3,2]:=n0[3,3,2]+Tg[k]*sqr(r*Cs[k])*h[2,2];
          n[2,3,3]:=n0[2,3,3]-Tg[k]*sqr(r*Cs[k])*h[3,3]; n[3,2,3]:=n[2,3,3];
          {Christoffel symbols}

          V[1]:=2*n[1,1,1]*(n[1,1,1]-1/r)+n[1,1,2]*(Tg[k]+
               2*n[2,2,2])+m[1,1]*h[2,2]/2-ri[1,1];
          V[1]:=V[1]+r*h0[2,2]*n[1,2,2]+n0[1,2,2]/r+h0[2,2];
          V[1]:=V[1]+r*sqr(Cs[k])*h0[3,3]*n[1,3,3]
               +n0[1,3,3]/r+sqr(Cs[k])*h0[3,3];
          for i:= 1 to 4 do
            for j:= 1 to 4 do begin
              if (i=3) and (j=4) or (i=4) and (j=3) then s:=-1 else s:=1;
              V[1]:=V[1]+s*n0[1,i,j]*n[1,j,i] end;

          V[1]:=2*g[1,1]*V[1];

          V[2]:=2*(-n[2,2,1]*n[1,1,1]+n0[2,2,1]/r
               -h0[1,1])+ri[2,2];
          V[2]:=V[2]-n[2,2,2]*(2*n[2,2,2]+Tg[k])+h[2,2]*m[2,2]/2
               -r*h0[2,2]*n[2,2,1]+r*h0[1,1]*n[2,1,2];
          V[2]:=V[2]+Tg[k]*sqr(r*Cs[k])*h0[3,3]*n[2,3,3]
               -n0[2,2,1]/r+h0[1,1];
          V[2]:=V[2]+r*(n0[2,1,2]+r*h0[2,2])+Tg[k]*(n0[2,3,3]
               -Tg[k]*sqr(r*Cs[k])*h0[3,3]);
          for i:=1 to 4 do
            for j:=1 to 4 do begin
              if (i=3) and (j=4) or (i=4) and (j=3) then s:=-1 else s:=1;
              V[2]:=V[2]-s*n0[2,i,j]*n[2,j,i] end;

          V[2]:=2*g0[1,1]+2*g[1,1]*V[2];

          V[4]:=2*g[1,1]*(2*n[4,4,1]*(1/r-n[1,1,1])-n[4,4,2]*(Tg[k]+
          2*n[2,2,2])-m[4,4]*h[2,2]/2+ri[4,4]);
          for i:=1 to 4 do
            for j:= 1 to 4 do begin
              if (i=3) or (j=3)
              then s:=-1 else s:=1;
              V[4]:=V[4]-s*2*g[1,1]*n[4,i,j]*n[4,j,i] end;

          V[5]:=2*g[1,1]*(2*n[3,4,1]*(1/r-n[1,1,1])-n[3,4,2]*(Tg[k]+2*n[2,2,2]
          )-m[3,4]*h[2,2]/2+ri[3,4]);
          for i:=1 to 4 do
            for j:=1 to 4 do V[5]:=V[5]-2*g[1,1]*n[3,i,j]*n[4,j,i];

          for w:=6 to 7 do
            begin V[w]:=g[1,1]*(2*n[1,1,1]*h[1,1]*f[w-3,1]+2*n[2,2,2]*h[2,2]*
            f[w-3,2]-q[w-3]*h[2,2]);
            for i:=1 to 4 do
              for j:=1 to 4 do
                for o:=1 to 4 do begin
                 if (w=6) and ((i=4) or (j=4) or (o=4)) then s:=-1 else s:=1;
                 V[w]:=V[w]+g[1,1]*(s*n[w-3,i,j]*f[j,o]
                 -n[i,j,j]*f[w-3,o])*h[o,i] end
            end;
          {2. derivations to  r}

          for i:=1 to 4 do
            if (i=1) or (i=4) then
            U[i,4,k]:=2*U[i,2,k]-U[i,0,k]+V[i]*4*dr*dr;
          U[2,4,k]:=2*U[2,2,k]-U[2,0,k]+sqr(2*dr/r)*(V[2]-
          2*U[2,2,k]-2*r*(U[2,1,k]-U[2,3,k])/dr);
          U[7,4,k]:=2*U[7,2,k]-U[7,0,k]+sqr(2*dr)*V[7]*r-
          8*U[7,2,k]*sqr(dr/r)+4*(U[7,1,k]-U[7,3,k])/dr/r;
          for i:=5 to 6 do U[i,4,k]:=2*U[i,2,k]-U[i,0,k]+(U[i,1,k]
          -U[i,3,k])*4*dr/r-U[i,2,k]*8*sqr(dr/r)+V[i]*sqr(2*dr/Cs[k])*r;
          U[3,4,k]:=-((1+U[1,4,k])*(1+U[2,4,k])*(U[4,4,k]+sqr(U[5,4,k]
          *Cs[k]/r/r))+U[2,4,k]*(1+U[1,4,k])+U[1,4,k]);
          U[3,4,k]:=U[3,4,k]/((1+U[4,4,k])*(1+U[1,4,k])*(1+U[2,4,k]))
          {the next values for the storing field from the last,
          they refer to the next radius step}
  end;
end;

begin { main program }
      ExitSave := ExitProc;
      ExitProc := @MyExit;
      {for the own exit from the program,
      in order to close all open data files}

with wr do begin

      write(#13#10'File No. (begin): '); ReadLn(x2);
      write('Search (U), Interrupted on (W), Store (S): ');
      ca:=ReadKey; WriteLn(ca);
      if upcase(ca)='W' then bo:=false else bo:=true;
      if upcase(ca)='U' then bn:=false else bn:=true;
      if upcase(ca)='S' then bs:=true else bs:=false;
  {In the mode "Search" (not bn) is documented only the number of steps
   until end of program (in fi1). An other character means "normal" run
   that is all data files except fi6 are opened.}

  if not bo then begin  {when interrupted program is continued}
    assign(fi5,'c:\data\fsp'+x2+'.dat'); reset(fi5);
    read(fi5,wr); close(fi5); erase(fi5);
  {reads all "spooled" values and continues with them}
    assign(fi1,'c:\data\fet'+x2+'.dat'); append(fi1);
    if bn then begin
    if bs then begin
      assign(fi6,'c:\data\fs'+x2+'.dat'); reset(fi6); seek(fi6,FileSize(fi6))
    end;
    assign(fi2,'c:\data\fe'+x2+'.dat'); reset(fi2); seek(fi2,FileSize(fi2));
    assign(fi3,'c:\data\fgr'+x2+'.dat'); reset(fi3); seek(fi3,FileSize(fi3))
    end
  end else begin
      write('Number of angle steps (z<=120)  z='); ReadLn(z); z0:=z/1;
      write('Initial radius  r0='); ReadLn(r0);
      write('Step size  dr='); ReadLn(dr);
      {difference of radius}
      write('Constant curvature  K0='); ReadLn(K0);
      if bn then begin
        write('Automatic (a):  '); x1:=ReadKey; WriteLn(x1)
      end else x1:='A';
      if upcase(x1)='A' then begin
      {automatic runs across the following given ranges of values
       and intervals of the integration constants
       ("Search" always automatically)}
        write('Smallest mass  m min = '); ReadLn(c01);
        write('Largest mass   m max = '); ReadLn(c10);
        write('step  '); ReadLn(c1s);
        write('Smallest spin  s min = '); ReadLn(c02);
        write('Largest spin   s max = '); ReadLn(c20);
        write('step  '); ReadLn(c2s);
        write('Smallest charge  Q min = '); ReadLn(c03);
        write('Largest charge   Q max = '); ReadLn(c30);
        write('step  '); ReadLn(c3s);
        write('Smallest magn.momentum  M min = '); ReadLn(c04);
        write('Largest magn.momentum   M max = '); ReadLn(c40);
        write('step  '); ReadLn(c4s);
        c1:=c01; c2:=c02; c3:=c03; c4:=c04
      end else begin  {1 run only}
        write('Mass  m='); ReadLn(c1);
        write('Spin  s='); ReadLn(c2);
        write('Charge  Q='); ReadLn(c3);
        write('Magn.momentum  M='); ReadLn(c4)
      end;
      write('Interrupt with ^U.'); delay(2000);

      k:=0; dt:=pi/(2*z-3);  {angle difference}
      repeat
        t:=(k-1)*dt;
        Cs[k]:=cos(t); Tg[k]:=sin(t)/Cs[k];
        inc(k)
      until k=z+1;
      {list for Cosinus and Tangens across the angle steps}
      for i:=1 to 4 do
        begin a[i]:=0; q[i]:=0; for j:=1 to 4 do
          begin g[i,j]:=0; h[i,j]:=0; p[i,j]:=0; m[i,j]:=0;
                g0[i,j]:=0; h0[i,j]:=0; m0[i,j]:=0;
          for o:=1 to 4 do l[i,j,o]:=0 end
        end;
      {all components set to zero, the not interesting keep zero always}
      z1:=(z-1)/1
  end;

repeat clrscr;
  writeln('z=',z:2,',r0=',r0:3:1,',dr=',dr:5:3,',K0=',K0:8,',');
  writeln('m=',c1:11,',s=',c2:11,',Q=',c3:11,',M=',c4:11);
  window(1,3,80,25);
  if bo then begin hf:=0;
      assign(fi1,'c:\data\fet'+x2+'.dat');
    if bn then begin
      if bs then begin
        assign(fi6,'c:\data\fs'+x2+'.dat'); rewrite(fi6) end;
      rewrite(fi1);
      assign(fi2,'c:\data\fe'+x2+'.dat'); rewrite(fi2);
      assign(fi3,'c:\data\fgr'+x2+'.dat'); rewrite(fi3);
      write(fi3,z1)
    end else
      if (c1=c01) and (c2=c02) and (c3=c03) and (c4=c04)
        then rewrite(fi1) else append(fi1);
      writeln(fi1,'z=',z:2,',r0=',r0:3:1,',dr=',dr:5:3,',K0=',K0:8,',');
      write(fi1,'m=',c1:11,',s=',c2:11,',Q=',c3:11,',M=',c4:11,',');
      {registration of the parameters}

      y:=0;
      repeat
        r:=r0-y*dr; k:=0;
        repeat
 {initial conditions}
          U[1,y,k]:=K0*r*r+c1/r-sqr(c3/r)/2+sqr(c4/r/r)*(1+sqr(Cs[k]))/10;
          U[2,y,k]:=sqr(c4/r/r)*(sqr(Cs[k])/3-0.3);
          U[3,y,k]:=sqr(c4/r/r)*(sqr(Cs[k])/15-0.3);
          U[4,y,k]:=-K0*r*r-c1/r+sqr(c3/r)/2+sqr(c4/r/r)*(1-sqr(Cs[k]))/2;
          U[5,y,k]:=c2-c3*c4/r/2;
          U[6,y,k]:=c4; U[7,y,k]:=c3;
          inc(k)
        until k=z+1; inc(y)
      until y=4 {beginning values set};

      for y:=0 to 3 do
        begin r:=r0-y*dr; for i:=1 to 4 do
          begin write(r:6:3,' G',i,i); for u0:=0 to 6 do
            begin k:=round(u0*(z-1.5)/6)+1;    {G means \gamma}
            a0:=U[i,y,k]; write(a0:8);
            if bn then write(fi2,a0) end
          end;
        write(r:6:3,' G34'); for u0:=0 to 6 do
          begin k:=round(u0*(z-1.5)/6)+1;
          a0:=U[5,y,k]*Cs[k]/r/r; write(a0:8);
          if bn then write(fi2,a0) end;
        write(r:6:3,'  A3'); for u0:=0 to 6 do
          begin k:=round(u0*(z-1.5)/6)+1;
          a0:=U[6,y,k]*Cs[k]/r/r; write(a0:8);
          if bn then write(fi2,a0) end;
        write(r:6:3,'  A4'); for u0:=0 to 6 do
          begin k:=round(u0*(z-1.5)/6)+1;
          a0:=U[7,y,k]/r; write(a0:8);
          if bn then write(fi2,a0) end
        end;

      for y:=1 to 7 do
        for k:=1 to 120 do begin Uh[y,k]:=0; Ug[y,k]:=0 end;

      y:=2
    end;

      repeat if bo then begin r:=r0-y*dr; k:=1 end;
        repeat
          bo:=true; {At this place starts the program, if an
          interrupted computation goes on}
          if KeyPressed then
            if ReadKey=#21 then begin
              assign(fi5,'c:\data\fsp'+x2+'.dat');
              rewrite(fi5); write(fi5,wr); close(fi5); halt
            end;
          {computation interrupted, all variables stored in fi5}

          kSchleife;

          inc(k)
        until k=z;

        for i:=1 to 7 do
          begin U[i,4,0]:=U[i,4,2];
            U[i,4,z]:=U[i,4,z-1]
          end   {"marginal values" loaded};
        for k:=0 to z do
          for i:=1 to 7 do begin
            if bs and bn then write(fi6,U[i,4,k]);
            for j:=0 to 3 do U[i,j,k]:=U[i,j+1,k]
          end   {shifted};

        r:=r-2*dr;
        for k:=1 to z-1 do begin
          ul[1]:=K0*r*r+c1/r-sqr(c3/r)/2+sqr(c4/r/r)*(1+sqr(Cs[k]))/10;
          ul[2]:=sqr(c4/r/r)*(sqr(Cs[k])/3-0.3);
          ul[3]:=sqr(c4/r/r)*(sqr(Cs[k])/15-0.3);
          ul[4]:=-K0*r*r-c1/r+sqr(c3/r)/2+sqr(c4/r/r)*(1-sqr(Cs[k]))/2;
          for i:=1 to 4 do
            begin if abs(U[i,4,k])>1e12 then hf:=y;
            if (Uh[i,k]=0) and (U[i,4,k]<-0.9) then Uh[i,k]:=r;
            if (Ug[i,k]=0) and (abs(U[i,4,k])>2*abs(ul[i]))
            and (ul[i]<>0) then Ug[i,k]:=r end;
          for i:=5 to 6 do
            begin if abs(U[i,4,k])*Cs[k]>1e12*r*r then hf:=y;
            if (Uh[i,k]=0) and (abs(U[i,4,k])*Cs[k]>0.9*r*r) then
            Uh[i,k]:=r end;
          if (Ug[5,k]=0) and (abs(U[5,4,k])>2*abs(c2-c3*c4/r/2))
          and (c2-c3*c4/r/2<>0) then Ug[5,k]:=r;
          if (Ug[6,k]=0) and (abs(U[6,4,k])>2*abs(c4))
          and (c4<>0) then Ug[6,k]:=r;
          if (Ug[7,k]=0) and (abs(U[7,4,k])>2*abs(c3))
          and (c3<>0) then Ug[7,k]:=r;
          if (Uh[7,k]=0) and (abs(U[7,4,k])>0.9*r) then Uh[7,k]:=r;
          if abs(U[7,4,k])>1e12*r then hf:=y end;
        {registration of the limiting radii for graphics (Uh: geometr. limit,
         Ug: entry into the non-linearity)  and setting the conditions
         for exit}

        for i:=1 to 4 do
          begin write(r:6:3,' G',i,i); for u0:=0 to 6 do
            begin k:=round(u0*(z-1.5)/6)+1;
            a0:=U[i,4,k]; write(a0:8);
            if bn then write(fi2,a0) end
          end;
        write(r:6:3,' G34'); for u0:=0 to 6 do
          begin k:=round(u0*(z-1.5)/6)+1;
          if abs(r)<dr/2 then a0:=0 else
          a0:=U[5,4,k]*Cs[k]/r/r; write(a0:8);
          if bn then write(fi2,a0) end;
        write(r:6:3,'  A3'); for u0:=0 to 6 do
          begin k:=round(u0*(z-1.5)/6)+1;
          if abs(r)<dr/2 then a0:=0 else
          a0:=U[6,4,k]*Cs[k]/r/r; write(a0:8);
          if bn then write(fi2,a0) end;
        write(r:6:3,'  A4'); for u0:=0 to 6 do
          begin k:=round(u0*(z-1.5)/6)+1;
          if abs(r)<dr/2 then a0:=0 else
          a0:=U[7,4,k]/r; write(a0:8);
          if bn then write(fi2,a0) end;
        {output and registration of the values all 15deg}

        inc(y);   {radius step}
        if hf>0 then y:=round(r0/dr)
      until y=round(r0/dr);
      if hf=0 then writeln(fi1,y+1:4) else writeln(fi1,hf+2:4);
      {registration of the number of steps}
      close(fi1);

    if bn then begin
      for y:=1 to 7 do begin
        for k:=1 to z-1 do write(fi3,Ug[y,k]);
        for k:=1 to z-1 do write(fi3,Uh[y,k]) end;
      if bs then begin
        write(fi6,z0,z0,z0,z0,z0,z0,z0); close(fi6) end;
      close(fi2); close(fi3)
    end;
    {all data saved}

      if upcase(x1)='A' then begin  {switching up the
           integration constants during "automatic"}
        c1:=c1+c1s; if c1>c10 then begin c1:=c01; c2:=c2+c2s end;
        if c2>c20 then begin c2:=c02;c3:=c3+c3s end;
        if c3>c30 then begin c3:=c03;c4:=c4+c4s end;
        if bn then begin
          val(x2,xi,code); inc(xi); str(xi,x2)
        end else code:=0
      end;
      window(1,1,80,25)
    until (c4>c40) or (upcase(x1)<>'A') or (code<>0);

    reset(fi1); if bn then begin
      reset(fi2); reset(fi3);
      if bs then reset(fi6) end
    {for the last exit}
  end;
end.