DataMuseum.dk

begin
real al,c,d,r1,r2,r,gamma,delta,dr,t1,t2,ke,pe;
integer i,j,k,gem,itr;
boolean lay;
real  procedure expint(y);
value y; real y;
begin
real u,v,b,t,a0,a1,b0,b1;
integer i,c,p,a;
if y>1 then
  begin
  a0:=0; b0:=a1:=1; b1:=y;
  c:=p:=0;
  v:=(if abs(y)<10-600 then 10600 else 1/y);
rep:  b:=y⨯⨯p;
   a:=p+c;
     t:=b⨯b0+a⨯a0; a0:=b0; b0:=t;
    t:=b⨯b1+a⨯a1; a1:=b1; b1:=t;
   u:=(if abs(b0)⨯10-154<abs(b1) then b0/b1 else 10154);
   if abs(u/v-1)>510-11 then begin v:=u; goto rep end;
expint:=u⨯exp(-y);
   end else
expint:=((((.00107857⨯y-.00976004)⨯y+.05519968)⨯y
-.24991055)⨯y+.99999193)⨯y-.57721566-ln(y);
end expint;
procedure min2(x,y,dr,delta,f,itr);
value dr,delta; real x,y,dr,delta;
integer itr;
real procedure f;
begin
real p1,p2,p3,p4,p0,dx,dy;
itr:=0;  dx:=dr⨯x; dy:=dr⨯y;
LI:
p0:=f(x,y);
L:
p1:=f(x,y+dy); p2:=f(x-dx,y);
p3:=f(x,y-dy); p4:=f(x+dx,y);
LL:
itr:=itr+1;
if p0<p1∧p0<p2∧p0<p3∧p0<p4 then
  begin
  dx:=dx/2; dy:=dy/2;
  if abs(dx)>delta/2⨯abs(x)|abs(dy)>delta/2⨯abs(y) then goto L
  end else begin
    if p1<p0∧p1<p2∧p1<p3∧p1<p4 then begin y:=y+dy;
    p3:=p0; p0:=p1; p1:=f(x,y+dy);
      p2:=f(x-dx,y); p4:=f(x+dx,y) end else
    if p2<p0∧p2<p1∧p2<p3∧p2<p4 then begin x:=x-dx;
    p4:=p0; p0:=p2; p1:=f(x,y+dy);
      p2:=f(x-dx,y); p3:=f(x,y-dy) end else
    if p3<p0∧p3<p2∧p3<p1∧p3<p4 then begin y:=y-dy;
    p1:=p0; p0:=p3; p2:=f(x-dx,y);
      p3:=f(x,y-dy); p4:=f(x+dx,y) end else begin
      x:=x+dx; p2:=p0; p0:=p4; p1:=f(x,y+dy);
      p3:=f(x,y-dy); p4:=f(x+dx,y) end;
    goto LL
  end
end min2;
real procedure wein(al,c);
real al,c;
begin
real ro1,ro2,ro3,e1ro,e2ro,s,s1,s2,ke0,ke1,pe0,pe1
     ,expi2,expi4,cou,ion,e1,e2,ex1;
ro1:=al⨯r;
  ro2:=ro1⨯⨯2; ro3:=ro2⨯ro1;
  e1ro:=exp(-ro1);
  e2ro:=e1ro⨯⨯2;
expi2:=expint(2⨯ro1);
expi4:=expint(4⨯ro1);
  s:=e1ro⨯(1+ro1 +ro2/3);
s1:=(1-e2ro⨯(1+ro1))/ro1 + 1;
  ke0:=1+s⨯e1ro⨯(1 +ro1-ro2/3);
  ke1:=2⨯e1ro⨯(1+ro1);
  s2:=(1-ro1+ro2/3)/e1ro;
  e1:=((s⨯⨯2⨯(ln(ro1)+gamma)+2⨯s⨯s2⨯expi2)
      -s2⨯⨯2⨯expi4)/ro1;
  e2:=e2ro⨯(.625 -1.15⨯ro1 -.6⨯ro2 -ro3/15);
  ex1:=1.2⨯e1+e2;
cou:=(1-e2ro⨯(1 +1.375⨯ro1 +.75⨯ro2 +ro3/6))/ro1;
ion:=(e1ro⨯(5 +2⨯ro1 +16⨯ro2) -e2ro⨯e1ro⨯(5+2⨯ro1))/8/ro1;
pe0:= -2⨯(s1 + s⨯ke1) + ex1;
pe1:= -2⨯ke1 -2⨯s⨯(s1) + ion;
d:=2⨯(c⨯⨯2⨯((s-1)⨯⨯2)+(s+1)⨯⨯2);
ke:=2⨯al⨯al⨯(c⨯⨯2⨯(ke0-ke1)+ke0+ke1);
pe:=al⨯(c⨯⨯2⨯(2⨯pe0-2⨯pe1+.625+cou)+2⨯c⨯(.625-cou)+2⨯pe0+2⨯pe1+.625+cou);
wein:=(pe+ke)/d
end wei;
lay:=false add(1 shift 5);
read(in,r1,dr,r2,delta,al,c);
gamma:=.577215665;
k:=abs entier((r2-r1)/dr);
systime(1,0,t1);
write(out,<:<10>Weinbaum beregning:>,lay,10,<<  dd dd dd>,
          systime(2,t1,t2),t2,<:<10>:>,lay,4,<:r:>,lay,5,<:alfa:>,
          lay,8,<:c:>,lay,10,<:Ep:>,lay,8,<:Ek:>,lay,8,
          <:E:>,lay,9,<:norm:>,lay,5,<:iterationer<10>:>);
begin
zone wei(10⨯256,1,stderror);
integer array tail(1:10);
open(wei,4,<:wein:>,0);
tail(1):=2;
if monitor(40,wei,0,tail)<>0 then write(out,<:<10>ingen reservation:>);
if monitor(50,wei,0,tail)<>0 then write(out,<:<10>ikke permanent:>);
outrec(wei,8);
for i:=1 step 1 until k do
 begin
 r:=r1+(i-1)⨯dr;
  min2(al,c,.25,delta,wein,itr);
  write(out,<:<10>:>,lay,4,<<dd>,r,<<  -d.dddddd>,al,c,pe/d,ke/d,(pe+ke)/d+(if r
>0 then 1/r else 0),
        1/sqrt(2⨯d),<<  ddd>,itr);
setposition(out,0,0);
wei((i-1)⨯8+1):=al; wei((i-1)⨯8+2):=c; wei((i-1)⨯8+3):=d;
wei((i-1)⨯8+4):=ke; wei((i-1)⨯8+5):=pe; wei((i-1)⨯8+6):=r;
wei((i-1)⨯8+7):=1.0⨯ itr;
end forstatement;
close(wei,true)
end
end;[end]