|
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: 11520 (0x2d00)
Types: TextFile
Names: »tppf1«
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
└─⟦4334b4c0b⟧
└─⟦e3cb726bb⟧ »test«
└─⟦this⟧
program ppf1
c implicit real*8(a-h,o-z)
external in,out
zone in,out
c external *rbplot
dimension n1(9),n3(9)
c double precision a(6)
real a(6)
real sumu , cpu,time,time1,base
real xins(1100),xz(1100)
real xzz(1100),xans(1100)
complex b(6)
common /tvline/length,isy1,isy2
c read(in,31)n2,n0
call uforeadi(in,n2)
call uforeadi(in,n0)
read(in,32)n1
do 711 ib=1,8
711 n3(ib)=n1(ib)
n3(9)=n1(9)+1
32 format(9i3)
write(out,31)n2
n22=n2
if(n2-1)52,52,53
52 n22=2
31 format(3i4,i10)
c 53 read(in,55)xacc,g
53 call uforead(in,xacc)
call uforead(in,g)
write(out,55)xacc,g
55 format(f10.8,f10.6,f10.4)
c read(in,9)xl,d,rf,e,brod
call uforead(in,xl)
call uforead(in,d)
call uforead(in,rf)
call uforead(in,e)
call uforead(in,brod)
write(out,9)xl,d,rf,e,brod
i4=-4
bred=brod*0.466882*10.**i4
xmax=0.700323
del=xmax/1100.
do 49 jy=1,1100
49 xins(jy)=0.
c if the time measure is short you schuold
c compensate for the time used to call systime
c The cpu time will depend somewhat on activities
c on other process ( but yo will hardly see it )
c The real time used is highly dependant on
c other process.
c in order to obtain a accuracy of 0.1 millisecond
c the time is messured relative to base.
cpu = systime(1,0.0,base)
cpu = systime(1,base,time)
c if(n0)111,112,112
c 111 read(11,30)xins
c 111 call uforead(in,xins)
c go to 95
c 112 continue
do 211 ll=1,9
theto=float(ll-1)*10.
do 20 jj=1,n3(ll)
theta=float(jj-1)*10./float(n1(ll))+theto
theta=theta*3.14159/180.
si=sin(theta)
co=cos(theta)
do 20 kk=1,n2
phi=float(kk-1)*90./float(n22-1)
phi=phi*3.14159/180.
co2=cos(2.*phi)
si2=sin(2.*phi)
co1=cos(phi)
si1=sin(phi)
9 format(6f10.4)
do 11 i=1,6
alfa=0.
l=0
lim=0
x=0.
19 xh=0.01*rf*2.*d/(rf+2.*d)
call froq(a,b,x,d,e,si,co2)
a0=a(i)
x=x+xh
call froq(a,b,x,d,e,si,co2)
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 if(a0-a1)70,71,70
71 write(out,72)
72 format(14h first warning)
70 x=xh*(a0-rf)/(a0-a1)+x-xh
call froq(a,b,x,d,e,si,co2)
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-1100)23,23,11
23 call intens(xmax,b,xacc,x,d,e,si,co,si1,co1,co2,si2,xi,i)
alfa=x
xins(n)=xi/float(n3(ll))+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 froq(a,b,x,d,e,si,co2)
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-1100)24,24,11
24 call intens(xmax,b,xacc,x,d,e,si,co,si1,co1,co2,si2,xi,i)
alfa=x
xins(n)=xi/float(n3(ll))+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 froq(a,b,x,d,e,si,co2)
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-1100)25,25,11
25 call intens(xmax,b,xacc,x,d,e,si,co,si1,co1,co2,si2,xi,i)
alfa=x
xins(n)=xi/float(n3(ll))+xins(n)
902 l=l+1
if(l-6)39,11,11
39 x=x+d/2.
go to 19
8 a1=a(i)
if(a0-a1)73,74,73
74 write(out,75)
75 format(15h second warning)
x=x-xh
go to 101
73 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)101,61,61
101 call froq(a,b,x,d,e,si,co2)
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=lim+10
l=5
go to 6
28 if(x-2.*xmax)6,11,11
18 n=x/del
xalfa=x-alfa
if(xalfa-2.*xacc)62,62,905
905 if(n-1100)26,26,11
26 call intens(xmax,b,xacc,x,d,e,si,co,si1,co1,co2,si2,xi,i)
alfa=x
xins(n)=xi/float(n3(ll))+xins(n)
62 l=l+1
if(l-6)45,11,11
45 x=x+d/2.
go to 19
11 continue
20 continue
211 continue
c of time messurement procedure
cpu = systime(1,base,time1) - cpu
time = time1 - time
write(out,300) cpu,time
300 format(10x,'cpu= ',f9.4,' real time = ',f9.4)
if(n0)95,97,95
95 continue
do 42 jc=1,1100
42 xz(jc)=float(jc)*del
do 40 jh=1,1100
sumu=0.
do 41 jq=1,1100
xqh=(xz(jq)-xz(jh))/bred
axqh=1.+xqh*xqh
41 sumu=sumu+xqh*xins(jq)/(axqh*axqh)
40 xans(jh)=sumu
do 22 ixi=1,1100
22 xzz(ixi)=xz(ixi)*100./(4.66882*g)
if(n0)98,97,98
97 write(out,30)xins
go to 99
30 format(10e10.3)
98 jyy=1100
c call time1
c call rbplot(xans,jyy,xzz)
c call tvbgn('gt-40','lpenh')
c call tvrng('user',0.,ymin,7.5,ymax)
c call tvmap('y',xans,jyy,5)
c call tvtext(xz,xans,tekst,30)
c call tvgrph(xz,xans,jyy,'magnetic field in kg',
c >'intensity')
c call tvend
99 stop
end
subroutine froq(a,b,x,d,e,si,co2)
c implicit double precision(c-h,o-z)
c double precision a(6)
real a(6)
complex b(6)
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*co2)
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*co2
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
c
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)
szcc=sqrt(zcc)
szcb=sqrt(zcb)
z1=-za/4.+szcc
z2=-za/4.-szcc
z3=za/4.+szcb
z4=za/4.-szcb
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,si1,co1,co2,
>si2,xi,i)
c implicit double precision(a-h,o-z)
complex u(4),um(4),b(6),sz,sp,sm,sumu,qumu,snor1,snor2
complex to,s3,solk
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*co1*(1.+3.*(e1-d)*e
>/((e1-d)*(e1-d)-9.*co*co*x*x/4.))/xna
alfir=x*si*si1*(-1.)*(1.-3.*e*(e1-d)/((e1-d)
>*(e1-d)-9.*x*x*co*co/4.))/xna
u(2)=cmplx(alfrr,alfir)
x1r=-sqrt(3.)/2.*x*si*(co1*alfrr+si1*(-1.)*alfir)/(e01-e1)-
>sqrt(3.)*e/(e01-e1)
x1i=-sqrt(3.)/2.*x*si*(co1*alfir-si1*(-1.)*alfrr)/
>(e01-e1)
x4r=-(sqrt(3.)/2.*x*si*co1+sqrt(3.)*e*alfrr)/(e04-e1)
x4i=-(sqrt(3.)/2.*x*si*si1*(-1.)+sqrt(3.)*e*alfir)/(e04-e1)
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*co1*(1.+3.*(e2-d)*e
>/((e2-d)*(e2-d)-9.*co*co*x*x/4.))/xnb
alfii=x*si*si1*(-1.)*(1.-3.*e*(e2-d)/((e2-d)
>*(e2-d)-9.*x*x*co*co/4.))/xnb
um(2)=cmplx(alfri,alfii)
x1mr=-sqrt(3.)/2.*x*si*(co1*alfri+si1*(-1.)*alfii)/(e01-e2)-
>sqrt(3.)*e/(e01-e2)
x1mi=-sqrt(3.)/2.*x*si*(co1*alfii-si1*(-1.)*alfri)/
>(e01-e2)
um(1)=cmplx(x1mr,x1mi)
x4mr=-(sqrt(3.)/2.*x*si*co1+sqrt(3.)*e*alfri)/(e04-e2)
x4mi=-(sqrt(3.)/2.*x*si*si1*(-1.)+sqrt(3.)*e*alfii)/(e04-e2)
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))
apm=aimag(sp*conjg(sm))
amp=aimag(sm*conjg(sp))
rlpm=real(sm+sp)
alpm=aimag(sm-sp)
rz=real(sz)
xi=(0.5*si*si*rzz+(1.+co*co)/8.*(rmm+rpp)-
>si*si*co2/8.*(rmp+rpm)-si*si*si2/8.*(amp-apm)-
>0.5*si*co*rz*co1*rlpm-0.5*si*co*rz*si1*alpm)*si
return
end
▶EOF◀