|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 16128 (0x3f00)
Types: TextFile
Names: »pftx1«
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
└─⟦4334b4c0b⟧
└─⟦this⟧ »pftx1«
Fit to Planck function of temperature
begin
boolean cel,dow,nup,ter;
integer i,j,l,m,mm,mp,n,p,q,r,s,u,w,we,wf,wt;
long dT,D,F,S,T,T0,T1,T2,U;
real a,b,C,d,e,eng,f,fix,g,h,K,K1,K2,R,t0,v,x,y,y0,z;
array c(1:3,0:7),cx(1:3,1:6);
real procedure joke(o);
value o; integer o;
begin
f:=h:=0;
z:=(T-T0)/T;
v:=(T-T0)/(T+T0);
y:=v*v;
if y<.554 then begin
b:=T0/(T+T0);
for i:=15 step -1 until 3 do h:=i*i*y/(i+i+1-h);
g:=y*4/(5-h); e:=y/(3-g);
for i:=mp step -1 until o do case 8-i of begin
f:=((((v*15-h*12+48)*v+h*23-87)*v-h*4)*v+h)*c(m,i)/b/80;
f:=((((v*15-h*8+27)*v+h*15-60)*v-h)*c(m,i)/48 + f)*g/b;
f:=((((v*4+g*3-8)*v-g*2)*v+g)*c(m,i)/6 + f)/b;
f:=(((v*2+g*3-6)*v-g)*c(m,i)/4 + f)*e/b;
f:=((y-e)*c(m,i) - f)/b;
f:=((v+e)*c(m,i) + f)/(1-e);
f:=(c(m,i) + f)*z + c(m,0)
end;
joke:=case o of (f,f*z);
end else begin
for i:=mp-1 step -1 until 1 do f:=f*x + cx(m,i);
joke:=case o of ((c(m,1) + f)*z + ln(T0/T)*cx(m,1) + c(m,0),
f*z + ln(T0/T)*cx(m,1))
end
end;
real procedure rent(o);
value o; integer o;
begin
f:=0;
for i:=mp step -1 until o do f:=f*x + c(m,i);
rent:=f
end;
real procedure mhcp(o);
value o; integer o;
begin
f:=0;
for i:=mp step -1 until o do f:=f*x + c(m,i)*(i-1);
mhcp:=f
end;
procedure devi;
begin
if g=0 then g:='-12;
l:=entier(ln(g)/ln10 + 13);
q:=if m=mm then 12 else 10;
l:=if l>q then q else if l<0 then 0 else l;
if f<0 then begin
write(out,"sp",q-l);
for i:=1 step 1 until l do outchar(out,35)
end else write(out,"sp",q);
outchar(out,54);
if f>0 then begin
for i:=1 step 1 until l do outchar(out,35);
write(out,"sp",q-l)
end else write(out,"sp",q)
end;
procedure tbhd(o);
value o; integer o;
begin
write(out,<:
temp./:>,if cel then <:C:> else <:K:>,"sp",wt-7);
if o//100>0 then begin
we:=wf//100+2;
we:=we+we//3+(if o//100=1 then 5 else 1);
write(out,<: !:>,"sp",we-we//2,case o//100 of (<:K:>,
<:pK:>,<:Yø:>),"sp",if o mod 100=0 then 0 else we//2)
end;
if o mod 100>9 then begin
we:=wf mod 100//10+2;
we:=we+we//3-2;
write(out,<: !:>,"sp",we-we//2,<:Hø/'3:>,
"sp",if o mod 10=0 then 0 else we//2)
end;
if o mod 10>0 then begin
we:=wf mod 10+2;
we:=we+we//3;
write(out,<: !:>,"sp",we-we//2,<:Cpø:>)
end;
write(out,"nl",1)
end;
procedure tbln(o);
value o; integer o;
begin
firm(w);
if n=1 then write(out,"sp",wt) else write(out,"nl",1,string fix,
T/F-C,"sp",if w=0 then 2 else 0);
x:=(T-T0)/T0;
if o//100>0 then begin
if o//100=1 then flow(wf//100) else firm(wf//100);
write(out,<: !:>,if o//100=1 then string eng else string fix,
case o//100 of (exp(joke(1))*y0,-joke(1)/ln10,joke(1)*R))
end;
if o mod 100>9 then begin
firm(wf mod 100//10);
write(out,<: !:>,string fix,rent(1)*T0*R/(F*1000))
end;
if o mod 10>0 then begin
firm(wf mod 10);
write(out,<: !:>,string fix,mhcp(2)*R)
end
end;
procedure flow(o);
value o; integer o;
eng:=real(case o+1 of (<<-d00'-dd>,<<-dd0.0'-dd>,<<-ddd.00'-dd>,
<<-ddd.d00'-dd>,<<-ddd.dd0_0'-dd>,<<-ddd.ddd_00'-dd>,
<<-ddd.ddd_d00'-dd>,<<-ddd.ddd_dd0_0'-dd>,
<<-ddd.ddd_ddd_00'-dd>,<<-ddd.ddd_ddd_d00'-dd>));
procedure firm(o);
value o; integer o;
fix:=real(case o+1 of (<<-dddd>,<<-dddd.d>,<<-dddd.dd>,
<<-dddd.ddd>,<<-dddd.ddd_d>,<<-dddd.ddd_dd>,<<-dddd.ddd_ddd>,
<<-dddd.ddd_ddd_d>,<<-dddd.ddd_ddd_dd>,<<-dddd.ddd_ddd_ddd>));
for r:=readi(<:
input mode:>) while r>0 and r<5 do begin
nup:=r<4;
if nup then begin
n:=readi(case r of
(<: number of free energy values to be fitted:>,
<: number of enthalpy values to be fitted:>,
<: number of heat capacity values to be fitted:>));
if n<3-r then goto L0;
R:=if readi(<: SI ?:>) = 0 then readr(<: R:>)
else gasconstant;
cel:=readi(<:
Celsius ?:>) = 1;
C:=if cel then 273.15 else 0;
w:=readi(<: temp.decimals:>);
if w<0 or w>4 then goto L0;
firm(w); flow(2);
F:=10**w;
wt:=case w+1 of (7,7,8,9,11);
case r//2+1 of begin
begin
long array N(1:n);
array U,U2,UU,KU(1:3),k(1:n);
mp:=if n>3 then 3 else n-1;
mm:=if n>3 then 1 else 5-n;
u:=(n+1)//2;
dow:=readi(<: equidistant temp.s ?:>)=1;
if dow then begin
T1:=(readr(<: first temp.:>) + C)*F;
dT:=readr(<: temp.interval:>)*F;
if T1<=0 or dT<=0 then goto L0
end;
p:=readi(<:
free energy type:>);
if p<1 or p>3 then goto L0;
if p=1 then y0:=readr(<: K/Kø:>);
d:=case p of (ln(readr(<: K0:>)/y0),
-readr(<: pK0:>)*ln10,readr(<: Y0:>)/R);
ter:=readi(<: term.input?:>)=1;
if -,ter then connectcuri(<:pfind:>);
for j:=1 step 1 until n do begin
if ter then write(out,"nl",1,<<ddddd>,j,"sp",2);
if dow then begin
N(j):=T:=(j-1)*dT + T1;
if ter then write(out,string fix,T/F-C,
if cel then <:C :> else <:K :>);
if ter then setposition(out,0,0)
end else begin
if ter then setposition(out,0,0);
read(in,a);
N(j):=T:=(a+C)*F
end;
if j<=u then T0:=T;
if T+T0>11863283 then goto L0;
read(in,h);
if case p of (h<=0,abs h>'10,abs h>'10) then goto L0;
k(j):=case p of (ln(h),-h*ln10,h/R)
end;
if -,ter then unstackcuri;
for l:=1,2,3 do U(l):=U2(l):=UU(l):=KU(l):=0;
K:=k(u);
m:=1;
for j:=u+m while j>0 and j<=n do begin
m:=if m>0 then -m else 1-m;
T:=N(j);
D:=T-T0;
S:=T+T0;
f:=D/T;
y:=(D*D)/(S*S);
if y<.445 then begin
e:=0;
for i:=13 step -1 until 1 do e:=i*i*y/(i+i+1-e);
b:=1-e;
g:=(D/S+e)/b*f;
h:=(y-e)/(T0/S)/b*f
end else begin
g:=ln(T/T0) - f;
h:=(D*S)/(T*T0) - ln((T*T)/(T0*T0))
end;
e:=k(j);
U(1):=U(1) + f;
U(2):=U(2) + g;
U(3):=U(3) + h;
U2(1):=U2(1) + f*f;
U2(2):=U2(2) + g*g;
U2(3):=U2(3) + h*h;
KU(1):=KU(1) + e*f;
KU(2):=KU(2) + e*g;
KU(3):=KU(3) + e*h;
UU(1):=UU(1) + g*h;
UU(2):=UU(2) + f*h;
UU(3):=UU(3) + f*g;
K:=e+K
end;
b:=U(2); e:=U(3);
for l:=1,2,3 do begin
a:=U(l);
U2(l):=U2(l)*n - a*a;
KU(l):=KU(l)*n - a*K;
UU(l):=UU(l)*n - b*e;
b:=e; e:=a
end;
f:=U2(2)*U2(1) - UU(3)*UU(3);
g:=UU(1)*U2(1) - UU(2)*UU(3);
h:=KU(2)*U2(1) - KU(1)*UU(3);
for m:=3 step -1 until mm do begin
c(m,3):=cx(m,2):=if m<mm+2 then 0 else
((KU(3)*U2(1) - KU(1)*UU(2))*f - g*h)/
((U2(3)*U2(1) - UU(2)*UU(2))*f - g*g);
c(m,2):=if m<mm+1 then 0 else (h - c(m,3)*g)/f;
cx(m,1):=c(m,3)*2 - c(m,2);
c(m,1):=(KU(1) - c(m,2)*UU(3) - c(m,3)*UU(2))/U2(1);
c(m,0):=(K - c(m,1)*U(1) - c(m,2)*U(2) - c(m,3)*U(3))/n
end;
write(out,<:
number of parameters in H:
:>);
for m:=mm step 1 until 3 do begin
KU(m):=0;
write(out,<<d>,m-mm+1,"sp",case m of (mm*2+20,mm*2+18,0))
end;
write(out,"nl",1);
for j:=1 step 1 until n do begin
write(out,"nl",1,<<ddd>,j);
T:=N(j);
x:=(T-T0)/T0;
for m:=mm step 1 until 3 do begin
f:=k(j)-joke(1);
g:=f*f;
KU(m):=KU(m) + g;
devi
end
end;
write(out,"nl",2,"sp",3);
for m:=mm step 1 until 3 do c(m,0):=c(m,0) + d;
for m:=mm step 1 until if n>4 then 3 else 2 do
write(out,"sp",7,string eng,sqrt(KU(m)/(n-mp-m+2)),"sp",5)
end;
begin
long array M,N(0:6);
array KU,KV,F0(0:6),k(1:if n>0 then n else 1);
mp:=if n>6 then 7 else n+r-2;
s:=mp-r+1;
if s*9+n>79 then goto L0;
mm:=if n>2 then 1 else 6-r-n;
T0:=(readr(if cel then <: ref.temp./C:>
else <: ref.temp./K:>) + C)*F;
p:=readi(<:
free energy type:>);
if p<1 or p>3 then goto L0;
if p=1 then y0:=readr(<: K/Kø:>);
dow:=if r=3 then false
else readi(<: Hø(T0) known exactly ?:>)=1;
K1:=case p of (ln(readr(<: K(T0):>)/y0),
-readr(<: pK(T0):>)*ln10,readr(<: Yø(T0):>)/R);
K2:=if dow or r=3 then readr(<: Hø(T0):>) else 0;
T1:=if n<1 then T0 else (readr(<:
first temp.:>) + C)*F;
dT:=if n<2 then T0 else readr(<: temp.interval:>)*F;
if T1<=0 or dT<=0 then goto L0;
for i:=0 step 1 until s do begin
KU(i):=0;
N(i):=n*n-i*i
end;
ter:=readi(<: term.input?:>)=1;
if -,ter then connectcuri(<:pfind:>);
for j:=1 step 1 until n do begin
if ter then
write(out,"nl",1,<<ddddd>,j,"sp",2,string fix,((j-1)*dT+T1)/F - C,
if cel then <:C :> else <:K :>);
if ter then setposition(out,0,0);
read(in,b);
if abs b>'10 then goto L0;
a:=k(j):=if r=3 then b/R else (b-K2)*F/T0/R;
m:=j*2-n-1;
S:=0; T:=1;
for i:=0 step 1 until s do begin
KU(i):=KU(i) + a*T;
D:=S; S:=T;
u:=i+i+1;
T:=(u*m*S - N(i)*D)//(i*i+u)
end
end;
if -,ter then unstackcuri;
for i:=0 step 1 until s do begin
D:=if i=0 then n else D*N(i)//(i*i);
M(i):=D//(i+i+1);
KU(i):=KU(i)/M(i)
end;
U:=(T0-T1)*2 - dT*(n-1);
g:=U/2/dT;
h:=(U*U)/4/(dT*dT);
D:=N(0); S:=D*3;
T:=(S*5-230)*D+407;
U:=(D*5-110)*D+329;
for i:=0 step 1 until s do case i+1 of begin
F0(i):=1;
F0(i):=g*2;
F0(i):=h*3 - N(1)/4;
F0(i):=(h*10 - (S-7)/2)*g/3;
F0(i):=(h*7 - (S-13)/2)*h*5/12 + N(1)*N(3)/64;
F0(i):=((h*.3 - (D-7)/12)*h*7 + T/480)*g;
F0(i):=((h*11/60-(S-31)/48)*h + U/960)*h*7 - N(1)*N(3)*N(5)/2304
end;
f:=T0/dT;
for m:=3 step -1 until mm do begin
if dow then begin
a:=b:=0;
for i:=s step -1 until 0 do begin
d:=F0(i);
a:=KU(i)*d + a;
b:=d*d/M(i) + b
end;
d:=a/b
end;
for i:=0 step 1 until 6 do
KV(i):=if i>s then 0
else if dow then KU(i) - F0(i)/M(i)*d else KU(i);
for i:=0 step 1 until s do case i+1 of begin
begin
a:=0;
for j:=s step -1 until 0 do a:=F0(j)*KV(j) + a;
KV(i):=a
end;
KV(i):=f*(((h*1.1 - (S-31)/12)*h + U/480)*g*7*KV(6)
+ ((h*3 - (D-7)/2)*h*3.5 + T/480)*KV(5)
+ (h*7 - (S-13)/4)*g*5/3*KV(4)
+ (h*10 - (S-7)/6)*KV(3)
+ g*6*KV(2)
+ 2*KV(1));
KV(i):=f*f*(((h*11 - (S-31)/2)*h + U/240)*1.75*KV(6)
+ (h*3 - (D-7)/4)*g*7*KV(5)
+ (h*7 - (S-13)/12)*2.5*KV(4)
+ g*10*KV(3)
+ 3*KV(2));
KV(i):=f**3*((h*11 - (S-31)/4)*g*7/3*KV(6)
+ (h*3 - (D-7)/12)*7*KV(5)
+ g*35/3*KV(4)
+ 10/3*KV(3));
KV(i):=f**4*((h*11 - (S-31)/12)*1.75*KV(6)
+ g*10.5*KV(5)
+ 35/12*KV(4));
KV(i):=f**5*(g*7.7*KV(6) + 2.1*KV(5));
KV(i):=f**6*77/60*KV(6)
end;
c(m,0):=K1;
c(m,1):=if dow or r=3 then K2*F/T0/R else KV(0);
for i:=2 step 1 until mp do
c(m,i):=if r=3 then KV(i-2)/(i-1) else KV(i-1);
s:=s-1
end;
write(out,<:
number of parameters in H:
:>);
for m:=mm step 1 until 3 do begin
KU(m):=0;
write(out,<<d>,m+mp-3,"sp",case m of (mm*2+20,mm*2+18,0))
end;
write(out,"nl",1);
for j:=1 step 1 until n do begin
write(out,"nl",1,<<ddd>,j);
x:=((j-1)*dT + T1 - T0)/T0;
for m:=mm step 1 until 3 do begin
f:=if r=3 then k(j)-mhcp(2) else
if dow then k(j) - rent(2)*x else k(j) - rent(1);
g:=f*f;
KU(m):=KU(m) + g;
devi
end
end;
write(out,"nl",2,"sp",3);
s:=if r=3 or dow then n-mp+4 else n-mp+3;
for m:=mm step 1 until if s<4 then s-1 else 3 do
write(out,"sp",7,string eng,sqrt(KU(m)/(s-m)),"sp",5)
end end end;
if -,nup then begin
mp:=7; mm:=-3
end;
for q:=readi(<:
number of parameters in H:>) while q>mp+mm-4 and q<=mp do begin
if r<4 then flow(9) else begin
mp:=q; mm:=4-q
end;
m:=q-mp+3;
ter:=readi(<: term.coef.?:>)=1;
if -,ter and r=4 then connectcuri(<:pfcoef:>);
if -,ter and r<4 then connectcuro(<:pfcoef:>);
if r<4 or ter then write(out,<:
ref.temp./K = :>);
if nup then begin
t0:=T0/F;
firm(w);
write(out,string fix,t0,"nl",1)
end else begin
if ter then setposition(out,0,0);
read(in,t0)
end;
for l:=0 step 1 until q do begin
if r<4 or ter then write(out,"nl",1);
if ter then write(out,<: c(:>,<<d>,l,<:) = :>);
if r<4 then write(out,string eng,c(m,l)) else begin
if ter then setposition(out,0,0);
read(in,c(m,l))
end
end;
if -,ter and r=4 then unstackcuri;
if -,ter and r<4 then closeout;
if r>1 then begin
h:=0;
for i:=mp-1 step -1 until 1 do begin
h:=c(m,i+1)*i - h;
cx(m,i):=if i>1 then h/i/(i-1) else -h
end
end;
write(out,"nl",1);
for l:=readi(<:
output type:>) while l>0 and l<312 do begin
if nup then nup:=l//100<>1 or p=1;
p:=l//100;
wf:=readi(<: decimals:>);
if wf<0 or wf>999 then goto L1;
if nup then nup:=readi(<: new K/Kø,R,C!K,t.dec.?:>)<>1;
if-,nup then begin
if p=1 then y0:=readr(<: K/Kø:>);
R:=if readi(<: SI ?:>) = 0 then readr(<: R:>)
else gasconstant;
cel:=readi(<:
Celsius ?:>) = 1;
C:=if cel then 273.15 else 0;
w:=readi(<: temp.decimals:>);
if w<0 or w>9 then goto L1;
F:=10**w;
T0:=t0*F;
wt:=case w+1 of (7,7,8,9,11,12,13,15,16,17);
nup:=true
end;
if wf//100 + wf mod 100//10 + wf mod 10 + w
> (if p=1 then 32 else 34) then goto L1;
for n:=readi(<:
output mode:>) while n>0 and n<4 do case n of begin
begin
tbhd(l);
write(out,"nl",1);
setposition(out,0,0);
read(in,a);
for T:=(a+C)*F while T>0 do begin
tbln(l);
write(out,"nl",2);
setposition(out,0,0);
read(in,a)
end
end;
begin
T1:=(readr(<: first temp.:>) + C)*F;
dT:=readr(<: temp.interval:>)*F;
if T1<=0 or dT<=0 then goto L2;
T2:=(readr(<: last temp.:>) + C)*F;
if T2<T1 then goto L2;
ter:=readi(<: term.output?:>)=1;
if -,ter then connectcuro(<:pfud:>);
tbhd(l);
for T:=T1 step dT until T2 do tbln(l);
closeout
end;
begin
if l<100 then goto L1;
T1:=(readr(<: lower temp.:>) + C)*F;
if T1<=0 then goto L2;
D:=S:=T2:=(readr(<: higher temp.:>) + C)*F;
if T2<=T1 then goto L2;
x:=(T1-T0)/T0;
j:=sgn(rent(1));
s:=sgn(mhcp(2));
dow:=j<>s;
x:=(T2-T0)/T0;
u:=if j<>sgn(rent(1)) then 3
else if s<>sgn(mhcp(2)) and dow then 4 else 2;
d:=maxlong/2;
for j:=1 step 1 until u-2 do begin
T:=if j=1 and dow then T1 else T2;
for s:=0,s+1 while dT<>T do begin
x:=(T-T0)/T0;
a:=mhcp(3)*x;
b:=((a-rent(3)*x)*x-c(m,1))/(c(m,2)+a)*T0 + T0;
if b<.5 or b>d or s>99 then begin
u:=2;
goto L5
end;
dT:=T;
T:=b
end;
case j of begin
D:=T; S:=T
end
end;
L5: K1:=maxreal; K2:=-K1;
tbhd(p*100);
for j:=1 step 1 until u do begin
T:=case j of (T1,D,S,T2);
tbln(p*100);
if f>K2 then K2:=f;
if f<K1 then K1:=f
end;
write(out,"nl",3);
tbhd(l);
L4: write(out,"nl",1,"sp",wt+2);
setposition(out,0,0);
read(in,h);
if p=1 and h<=0 then goto L2;
K:=case p of (ln(h/y0),-h*ln10,h/R);
if K<K1 or K>K2 then goto L2;
K:=c(m,0) + (c(m,1) - K);
L3: write(out,"sp",4);
setposition(out,0,0);
read(in,a);
T:=(a+C)*F;
if T<T1 or T>T2 then goto L4;
for s:=0,s+1 while dT<>T do begin
x:=(T-T0)/T0;
a:=rent(2);
b:=(a*x + c(m,1))/(joke(2) + a*z + K)*T0;
if b<.5 or b>d or s>99 then goto L3;
dT:=T;
T:=b
end;
tbln(l);
write(out,"nl",2);
goto L3
end;
L2: end;
L1: end
end;
L0: end end
IN1: en▶EOF◀