comment programma R 543c, M.L.Potters, codenr. MLP 201161 / 1510,		MCP: SUM, FIXT, FLOT,	ALD ;
begin	integer i; real array fct[0: 159];
	real procedure theta (x); value x; integer x;
	theta:= if x <= -160 then -1 else if x < 0 then -fct[-x] else if x < 160 then fct[x] else 1;
	procedure TRIDIAG (T,r,n); value n; integer n; array T,r;
	begin	integer k,h; real w; n:= n - 1; for k:= 1 step 1 until n do
		begin	h:= 3 * k; w:= T[h - 2] / T[0]; if abs (w) < abs (T[h] / T[2]) then
			begin	T[h - 2]:= T[0]; T[0]:= T[h] - w * T[2]; T[h]:= T[2];
				if k != n then begin T[2]:= T[h + 2]; T[h + 2]:= 0 end;
				w:= r[k] - w * r[0]; r[k]:= r[0]; r[0]:= w
			end	else
			begin	T[0]:= T[2] - T[h] / w; if k != n then T[2]:= -T[h + 2] / w; r[0]:= r[0] - r[k] / w end
		end;
		r[0]:= r[0] / T[0]; for k:= n step -1 until 1 do
		begin	h:= 3 * k; w:= if k != n
			then (r[k] - T[h] * r[0] - T[h + 2] * r[k + 1]) / T[h - 2]
			else (r[k] - T[h] * r[0]) / T[h - 2]; r[k]:= r[0]; r[0]:= w
		end
	end;
	for i:= 0 step 1 until 159 do begin real u; u:= exp (i / 10); fct[i]:= (u - 1) / (u + 1) end;
Start:	if read > 0 then
	begin	integer mu; real A,C,D,E,S,n;
		A:= read; C:= read; D:= read; E:= read; S:= read; n:= read;
		for i:= 1 step 1 until 5 do NLCR; FIXT (1,3,A); FIXT (1,3,C); FIXT (1,3,D);
		NLCR; FIXT (1,3,E); FIXT (1,3,S); FIXT (1,3,n); NLCR; NLCR;
		for mu:= 0,2,4,8,16,32,64,128,1000 do
		begin	real m1,m2,p1,p2,q1,q2,t1,t2,u1,u2,v,w,Jc0,Jc1,Jc2,Jh0,Jh1;
			real array cp,hp,cm,hm[0: 1],M1,M2[0: 3];
			Jc0:= Jc1:= Jc2:= Jh0:= Jh1:= 0; if mu = 0 | mu = 1000 then
			begin	if mu = 0 then Jc0:= 80 / (A + (D - E) * S * (S + 1)) else
				begin	v:= A + (D - E) * S ** 2; w:= 2 * C * S * (1 - 2 * n); Jc0:= 40 / (v + w) + 40 / (v - w) end;
				Jc2:= Jc0 * 328.9868134; goto Typ
			end;
			m1:= (S + 0.5) / theta ((2 * S + 1) * mu) - 0.5 / theta (mu); v:= m1 / theta (mu); m2:= S * (S + 1) - v;
			p1:= -D * m1 * n; p2:= -D * m1 * (1 - n); t1:= -E * m1 * n; t2:= -E * m1 * (1 - n);
			u1:= -E * v * n; u2:= -E * v * (1 - n); v:= A + (D - E) *  m2 + D * v; w:= 2 * C * m1 * (1 - 2 * n);
			q1:= v + w; q2:= v - w; v:= t1 + t2 - u1 + u2; w:= t1 + t2 + u1 - u2; if n = 1 then
			begin	integer eta; real array c,h[-160: 160];
				real procedure int (F); real F; int:= SUM (eta,-160,160,(1 - theta (eta) ** 2) * F);
				for eta:= -160 step 1 until 160 do
				begin	real thp,tho,thm,wr1,wr2,wr3;
					thp:= theta (eta + mu); tho:= theta (eta); thm:= theta (eta - mu);
					wr1:= -p1 * tho + q2; wr2:= t1 * thp - u1; wr3:= (p1 * thp + q1) * wr1 - (-t1 * tho - u1) * wr2;
					h[eta]:= wr2:= wr2 / wr3; c[eta]:= wr2 + wr1 / wr3;
					wr1:= p1 * tho + q1; wr2:= -t1 * thm - u1; wr3:= (-p1 * thm + q2) * wr1 - (t1 * tho - u1) * wr2;
					wr2:= wr2 / wr3; h[eta]:= h[eta] - wr2; c[eta]:= c[eta] + wr2 + wr1 / wr3
				end;
				Jc0:= int (c[eta]); Jc1:= int (c[eta] * eta); Jc2:= int (c[eta] * eta ** 2);
				Jh0:= int (h[eta]); Jh1:= int (h[eta] * eta); goto Typ
			end;
			M1[0]:= M2[0]:= p1 - p2 + q1; M1[1]:= M2[1]:= t1 - t2 + u1 + u2;
			M1[2]:= M2[2]:= -t1 + t2 + u1 + u2; M1[3]:= M2[3]:= -p1 + p2 + q2;
			cp[0]:= cp[1]:= 1; TRIDIAG (M1,cp,2);
			hp[0]:= v * cp[1]; hp[1]:= w * cp[0]; TRIDIAG (M2,hp,2);
			M1[0]:= M2[0]:= -p1 + p2 + q1; M1[1]:= M2[1]:= -t1 + t2 + u1 + u2;
			M1[2]:= M2[2]:= t1 - t2 + u1 + u2; M1[3]:= M2[3]:= p1 - p2 + q2;
			cm[0]:= cm[1]:= 1; TRIDIAG (M1,cm,2);
			hm[0]:= -w * cm[1]; hm[1]:= -v * cm[0]; TRIDIAG (M2,hm,2);
			for i:= 0 step 1 until mu - 1 do
			begin	integer J,K; K:= 1 + (320 + i) % mu;
				begin	integer j,r,s; real array M[-2: 3 * K + 2],c,h[0: K];
					integer procedure eta; eta:= -160 - i + j * mu;
					real procedure int (F); real F; int:= SUM (j,0,K,(1 - theta (eta) ** 2) * F);
					procedure matrix;
					begin	s:= r; for j:= 0 step 1 until K do
						begin	real thp,thm; thp:= theta (eta + mu); thm:= theta (eta - mu); J:= 3 * j;
							M[J - 2]:= if s > 0 then t2 * thm + u2 else t1 * thm + u1;
							M[J + 2]:= if s > 0 then -t1 * thp + u1 else -t2 * thp + u2;
							M[J]:= if s > 0 then p1 * thp - p2 * thm + q1
									    else p2 * thp - p1 * thm + q2; s:= -s
						end
					end; for r:= +1,-1 do
					begin	integer z1,z2;
						matrix; z1:= if r > 0 then 1 else 0; z2:= if s > 0 then 0 else 1;
						c[0]:= 1 - M[-2] * cm[z1]; c[K]:= 1 - M[3 * K + 2] * cp[z2];
						for j:= 1 step 1 until K - 1 do c[j]:= 1; TRIDIAG (M,c,K + 1);
						matrix; h[0]:= M[-2] * (cm[z1] - hm[z1]) - M[2] * c[1];
						h[K]:= M[3 * K - 2] * c[K - 1] - M[3 * K + 2] * (cp[z2] + hp[z2]);
						for j:= 1 step 1 until K - 1 do
						begin	J:= 3 * j; h[j]:= M[J - 2] * c[j - 1] - M[J + 2] * c[j + 1] end;
						TRIDIAG (M,h,K + 1); Jc0:= Jc0 + int (c[j]); Jc1:= Jc1 + int (c[j] * eta);
						Jc2:= Jc2 + int (c[j] * eta ** 2); Jh0:= Jh0 + int (h[j]); Jh1:= Jh1 + int (h[j] * eta)
					end
				end
			end;
		Typ:	NLCR; FIXT (2,1,mu / 10); TAB; FLOT (5,Jc0 / 40); TAB; FLOT (5,Jc1 / 400); TAB; v:= Jc1 / Jc0; FLOT (5,v / 10);
			TAB; FLOT (5,mu * Jh0 / 400); TAB; w:= -v * Jc1 + Jc2 + mu * Jh1; FLOT (5,w / 4000); TAB;
			FLOT (5,w / (Jc0 * 328.9868134) - 1)
		end;
		goto Start
	end
end
q