|
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: 9984 (0x2700)
Types: TextFile
Names: »nnppf2«
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
└─⟦4334b4c0b⟧
└─⟦e3cb726bb⟧ »test«
└─⟦this⟧
c implicit real(a-h,o-z)
program ppf2
real a(6),xins(1500)
complex b(6)
zone in,out
external in,out
real a0,a1,alfa,ca,co,d,del,e,f,
> rf,si,theta,x,x1,x2,xacc,xalfa,xh,xi,xl,
> xm,xmax,xnum,xtra,cpu,time,time1,base
do 49 jy=1,1500
49 xins(jy)=0.
c unformated read
call uforead(in,xacc)
call uforead(in,theta)
c read(in,55)xacc,theta
a write(out,56)xacc,theta
56 format(10x,11h accuracy =,f10.8,8h theta =,f6.2)
55 format(f10.8,2f10.4)
c unformated read
call uforead(in,xl)
call uforead(in,d)
call uforead(in,rf)
call uforead(in,e)
c read(in,9)xl,d,rf,e
a write(out,57)d,e,rf
a write(out,59)
59 format(3h )
57 format(10x,4h d =,f10.4,4h e =,f10.4,12h frequency =,
>f10.4)
xmax=0.700323
xtra=5000./0.466882
del=xmax/1500.
c call time0
c if the time messure is short you schuld compensate
c for the time used to call systime. the cpu time
c will depend somewhat on the activities of other
c process. the real time used is highly dependant
c on other processes.
c in order to octain a accuracy of 0.1 millisecond
c the time is messured relative to a base.
cpu=systime(1,0.0,base)
cpu=systime(1,base,time)
theta=theta*3.14159/180.
si=sin(theta)
co=cos(theta)
9 format(6f10.4)
do 11 i=1,6
alfa=0.
a write(out,92)i
l=0
lim=0
x=0.
19 xh=0.01*rf*2.*d/(rf+2.*d)
call freq(a,b,x,d,e,si)
a0=a(i)
x=xh+x
call freq(a,b,x,d,e,si)
a1=a(i)
if(a0-rf)15,21,13
21 x=x-xh
go to 4
13 if(a1-a0)12,16,16
15 if(a1-a0)16,16,12
16 l=l+1
if(l-6)122,12,11
122 x=x+d/2.
go to 19
12 x=xh*(a0-rf)/(a0-a1)+x-xh
call freq(a,b,x,d,e,si)
ca=abs(rf-a(i))
if(ca-xacc*xmax)4,4,6
4 n=x/del
xalfa=x-alfa
if(xalfa-2.*xacc)910,910,911
911 if(n-1500)23,23,11
23 call intens(xmax,b,xacc,x,d,e,si,co,xi,i)
alfa=x
xins(n)=xi+xins(n)
xm=xtra*x
a write(out,30)xm,xins(n)
910 l=l+1
if(l-6)36,11,11
36 x=x+d/2.
6 xh=xh/xl
lim=lim+1
x=x-xh
call freq(a,b,x,d,e,si)
a write(out,89)a,x
ca=abs(rf-a(i))
if(ca-xacc*xmax)10,10,7
10 n=x/del
xalfa=x-alfa
if(xalfa-2.*xacc)900,900,901
901 if(n-1500)24,24,11
24 call intens(xmax,b,xacc,x,d,e,si,co,xi,i)
alfa=x
xins(n)=xi+xins(n)
xm=xtra*x
a write(out,30)xm,xins(n)
900 l=l+1
if(l-6)37,11,11
37 x=x+d/2.
go to 19
7 a0=a(i)
x=x+2.*xh
call freq(a,b,x,d,e,si)
a write(out,89)a,x
ca=abs(rf-a(i))
if(ca-xacc*xmax)14,14,8
14 n=x/del
xalfa=x-alfa
if(xalfa-2.*xacc)902,902,903
903 if(n-1500)25,25,11
25 call intens(xmax,b,xacc,x,d,e,si,co,xi,i)
alfa=x
xins(n)=xi+xins(n)
xm=xtra*x
a write(out,30)xm,xins(n)
902 l=l+1
if(l-6)39,11,11
39 x=x+d/2.
go to 19
8 a1=a(i)
x1=x
x=2.*xh*(a0-rf)/(a0-a1)+x-2.*xh
x2=x
xnum=abs((x1-x2)/(d+x1+x2))
if(xnum-0.2)124,61,61
124 call freq(a,b,x,d,e,si)
a write(out,89)a,x
a write(out,92)lim,l
89 format(7e10.4)
if(lim-20)60,61,61
61 lim=0
x=2.6*d
l=l+1
if(l-6)19,11,11
60 ca=abs(rf-a(i))
if(ca-xacc*xmax)18,18,27
27 if(x)123,28,28
123 x=2.6*d
lim=10+lim
l=5
go to 6
28 if(x-2.*xmax)6,11,11
18 n=x/del
xalfa=x-alfa
if(xalfa-2.*xacc)904,904,905
905 if(n-1500)26,26,11
26 call intens(xmax,b,xacc,x,d,e,si,co,xi,i)
alfa=x
xins(n)=xi+xins(n)
xm=xtra*x
a write(out,30)xm,xins(n)
a write(out,92)i,l
92 format(2i3)
904 l=l+1
if(l-6)45,11,11
45 x=x+d/2.
go to 19
11 continue
c call of time procedure
c call time1
cpu = systime(1,base,time1) - cpu
time = time1 - time
write(out,300) cpu,time
300 format(10x,'cpu= ',f9.4,' real time = ',f9.4)
30 format(10x,13h h in gauss =,e10.4,12h intensity =,e10.4)
stop
end
subroutine freq(a,b,x,d,e,si)
c implicit real(a-h,o-z)
complex b(6)
real a(6)
real c2,c3,c4,cxx,d,dd,dsi,dt,e,e10,esi,
> o,q,r,si,x,x1,xh,xx,z,z1,z2,z3,z4,
> za,zaa,zb,zbb,zc,zcb,zcc,zzcb,zzcc
integer h
xh=x/2.
c2=-(2.*d*d+6.*e*e+10.*xh*xh)
c3=-(16.*xh*d*xh-24.*xh*xh*d*si*si+24.*xh*xh*e
>*si*si)
c4=d*d*d*d+6.*d*d*e*e+9.*e*e*e*e
>-10.*xh*xh*d*d+12.*xh*xh*d*d
>*si*si+9.*xh*xh*xh*xh+18.*e*e*xh*xh-36.
>*e*e*xh*xh*si*si+24.*xh*xh*d*e*si*si
q=-(12.*c4+c2*c2)/9.
r=(-72.*c2*c4+27.*c3*c3+2.*c2*c2*c2)/54.
dd=sqrt(-q)
xx=r/(dd*q)
cxx=abs(xx)
if(1.-cxx)12,12,11
12 dt=acos(-1.)
go to 13
11 dt=acos(r/(dd*q))
13 dt=dt/3.
x1=-2.*dd*cos(dt)+c2/3.
zaa=abs(4.*x1-4.*c2)
zbb=abs(x1*x1-4.*c4)
za=sqrt(zaa)
zb=sqrt(zbb)/2.
zc=-(x1+c2)/4.
if(c3)1,1,2
2 zb=-zb
1 zcc=abs(zc-zb)
zcb=abs(zc+zb)
zzcc=sqrt(zcc)
zzcb=sqrt(zcb)
z1=-za/4.+zzcc
z2=-za/4.-zzcc
z3=za/4.+zzcb
z4=za/4.-zzcb
a(1)=abs(z1-z2)
a(2)=abs(z1-z3)
a(3)=abs(z1-z4)
a(4)=abs(z2-z3)
a(5)=abs(z2-z4)
a(6)=abs(z3-z4)
b(1)=cmplx(z1,z2)
b(2)=cmplx(z1,z3)
b(3)=cmplx(z1,z4)
b(4)=cmplx(z2,z3)
b(5)=cmplx(z2,z4)
b(6)=cmplx(z3,z4)
return
end
subroutine intens(xmax,b,xacc,x,d,e,si,co,xi,i)
c implicit real(a-h,o-z)
complex u(4),um(4),b(6),sz,sp,sm,sumu,qumu,snor1,snor2
complex to,s3,solk
real a,a1,a2,a3,a4,aimag,alfii,alfir,alfri,
> alfrr,axna,axnb,b1,b2,b3,b4,
> co,d,e,e01,e02,e03,e04,e1,
> e2,flk,h,o,qrt,rlpm,rmm,rmp,rpm,
> rpp,rz,rzz,sh3,si,sis,
> x,x1i,x1mi,x1mr,x1r,x4i,x4mi,x4mr,
> x4r,xacc,xi,xm1r,xm4r,xmax,xna,xnb,xucc,z
xucc=xacc/10.
do 2 lj=1,4
u(lj)=cmplx(0.,0.)
2 um(lj)=cmplx(0.,0.)
e1=real(b(i))
e2=aimag(b(i))
e01=1.5*x*co+d
e02=0.5*x*co-d
e03=-0.5*x*co-d
e04=-1.5*x*co+d
a1=abs(e1-e01)
a2=abs(e1-e02)
a3=abs(e1-e03)
a4=abs(e1-e04)
b1=abs(e2-e01)
b2=abs(e2-e02)
b3=abs(e2-e03)
b4=abs(e2-e04)
if(a1-xucc*xmax)4,4,3
4 u(1)=cmplx(1.,0.)
go to 11
3 if(a2-xucc*xmax)6,6,5
6 u(2)=cmplx(1.,0.)
go to 11
5 if(a3-xucc*xmax)8,8,7
8 u(3)=cmplx(1.,0.)
go to 11
7 if(a4-xucc*xmax)10,10,9
10 u(4)=cmplx(1.,0.)
go to 11
9 u(3)=cmplx(1.,0.)
xna=3.*x*x*si*si/(4.*(e1-e01))+3.*e*e/(e1-e04)+e02-e1
sis=abs(si)
i4=-4
if(sis-10.**i4)31,31,30
31 axna=abs(xna)
axna=axna*e*e
i8=-8
if(axna-10.**i8)41,41,40
41 u(3)=cmplx(0.,0.)
u(1)=cmplx(0.,0.)
u(2)=cmplx(1.,0.)
x4r=-sqrt(3.)*e/(e04-e1)
u(4)=cmplx(x4r,0.)
go to 42
40 u(2)=cmplx(0.,0.)
u(4)=cmplx(0.,0.)
x1r=-sqrt(3.)*e/(e01-e1)
u(1)=cmplx(x1r,0.)
go to 42
30 alfrr=-x*si*(1.+3.*(e1-d)*e
>/((e1-d)*(e1-d)-9.*co*co*x*x/4.))/xna
alfir=0.0
u(2)=cmplx(alfrr,alfir)
x1r=-sqrt(3.)/2.*x*si*alfrr/(e01-e1)-
>sqrt(3.)*e/(e01-e1)
x1i=0.0
x4r=-(sqrt(3.)/2.*x*si+sqrt(3.)*e*alfrr)/(e04-e1)
x4i=0.0
u(1)=cmplx(x1r,x1i)
u(4)=cmplx(x4r,x4i)
42 sumu=cmplx(0.,0.)
do 21 jak=1,4
21 sumu=sumu+u(jak)*conjg(u(jak))
snor1=cmplx(1.,0.)/csqrt(sumu)
do 22 juk=1,4
22 u(juk)=snor1*u(juk)
11 if(b1-xucc*xmax)14,14,13
14 um(1)=cmplx(1.,0.)
go to 12
13 if(b2-xucc*xmax)16,16,15
16 um(2)=cmplx(1.,0.)
go to 12
15 if(b3-xucc*xmax)18,18,17
18 um(3)=cmplx(1.,0.)
go to 12
17 if(b4-xucc*xmax)20,20,19
20 um(4)=cmplx(1.,0.)
go to 12
19 um(3)=cmplx(1.,0.)
xnb=3.*x*x*si*si/(4.*(e2-e01))+3.*e*e/(e2-e04)+e02-e2
sis=abs(si)
i4=-4
if(sis-10.**i4)34,34,33
34 axnb=abs(xnb)
axnb=axnb*e*e
i8=-8
if(axnb-10.**i8)43,43,44
43 um(3)=cmplx(0.,0.)
um(1)=cmplx(0.,0.)
um(2)=cmplx(1.,0.)
xm4r=-sqrt(3.)*e/(e04-e2)
um(4)=cmplx(xm4r,0.)
go to 45
44 um(2)=cmplx(0.,0.)
um(4)=cmplx(0.,0.)
xm1r=-sqrt(3.)*e/(e01-e2)
um(1)=cmplx(xm1r,0.)
go to 45
33 alfri=-x*si*(1.+3.*(e2-d)*e
>/((e2-d)*(e2-d)-9.*co*co*x*x/4.))/xnb
alfii=0.0
um(2)=cmplx(alfri,alfii)
x1mr=-sqrt(3.)/2.*x*si*alfri/(e01-e2)-
>sqrt(3.)*e/(e01-e2)
x1mi=0.0
um(1)=cmplx(x1mr,x1mi)
x4mr=-(sqrt(3.)/2.*x*si+sqrt(3.)*e*alfri)/(e04-e2)
x4mi=0.0
um(4)=cmplx(x4mr,x4mi)
45 qumu=cmplx(0.,0.)
do 23 jik=1,4
23 qumu=qumu+um(jik)*conjg(um(jik))
snor2=cmplx(1.,0.)/csqrt(qumu)
do 24 jok=1,4
24 um(jok)=snor2*um(jok)
12 sz=cmplx(0.,0.)
do 1 lk=1,4
flk=float(lk)
solk=cmplx(2.5-flk,0.)
1 sz=sz+solk*u(lk)*conjg(um(lk))
sh3=sqrt(3.)
s3=cmplx(sh3,0.)
to=cmplx(2.,0.)
sm=s3*conjg(um(2))*u(1)+to*conjg(um(3))*u(2)+
>s3*conjg(um(4))*u(3)
sp=s3*conjg(um(1))*u(2)+to*conjg(um(2))*u(3)+
>s3*conjg(um(3))*u(4)
rzz=real(sz*conjg(sz))
rmm=real(sm*conjg(sm))
rpp=real(sp*conjg(sp))
rpm=real(sp*conjg(sm))
rmp=real(sm*conjg(sp))
rlpm=real(sm+sp)
rz=real(sz)
c xi=(si*si*rzz+co*co/4.*(rmm+rpp)+co*co/4.*(rmp+rpm)
c >-si*co*rz*rlpm)*4.0
xi=(0.5*si*si*rzz+(1.+co*co)/8.*(rmm+rpp)-si*si/8.*(rmp+rpm)
>-0.5*si*co*rz*rlpm)*4.0
return
end
▶EOF◀