9 implicit real*8(a-
h,o-
z)
29 +
i2/1,1,1,2,3,3,1,2/,
40 SUBROUTINE radgen(e1,q2l,nul,ys,xs,phil,PhRAD,q2tr,anutr,WEIGHT)
57 include
"radgenkeys.inc"
59 REAL phrad(*),q2tr,anutr
60 REAL e1,q2l,nul,phil,xs,ys,
weight
74 write(*,*)
"radgen: ",ys, q2l, xs, nul, phil, e1
77 if(kill_elas_res.ne.2)
then
78 call
itafun(ys,xs,plrun,pnrun,ixytest,itagen)
87 if(itagen.ne.0.and.abs(ixytest).ne.2)
then
89 call
mpolrad(e1,ys,xs,1.,plrun,pnrun,itagen)
90 if(ixytest.eq.1)ixytest=-2
92 call
radgam_pol(e1,ys,xs,phil,ixytest,itagen,q2tr,anutr)
95 phrad(1)=
sngl(dplabg(1))
96 phrad(2)=
sngl(dplabg(2))
97 phrad(3)=
sngl(dplabg(3))
100 phrad(4)=sqrt(phrad(1)*phrad(1)+phrad(2)*phrad(2)+
104 if(kill_elas_res.eq.2)
then
131 dimension db(3),dx(31),
dy(31),
dz(31)
147 if(
i.eq.1.or.
j.ne.1)
k=
k+1
165 + ixytest,itagen,q2tr,anutr)
172 include
"tailcom.inc"
173 include
"amf2com.inc"
174 include
"density.inc"
182 rr1=
r1*(tbor+tine+tpro+tnuc)
184 if (rr1.gt.(tine+tpro+tnuc))
then
186 elseif(rr1.gt.(tpro+tnuc))
then
189 call
conk2(dble(erad),dble(xrad),dble(yrad),1)
190 elseif(rr1.gt.tnuc)
then
193 call
conk2(dble(erad),dble(xrad),dble(yrad),1)
197 call
conk2(dble(erad),dble(xrad),dble(yrad),2)
203 if(ita.eq.1)sicurr=
r1*tine
204 if(ita.eq.2)sicurr=
r1*tnuc
205 if(ita.eq.3)sicurr=
r1*tpro
207 call
conk2(dble(erad),dble(xrad),dble(yrad),2)
209 call
conk2(dble(erad),dble(xrad),dble(yrad),1)
228 r1=
r1*dsitkm(ndxtkm(ita),ita)
231 if (
r1.le.dsitkm(ktk,ita))
then
233 dtk1=
r1/dsitkm(ktk,ita)
235 dtk1=(
r1-dsitkm(ktk-1,ita))/
236 + (dsitkm(ktk,ita)-dsitkm(ktk-1,ita))
238 dtk=dcvtkm(ktk,ita)+(dtk1-0.5d0)*ddetkm(ktk,ita)
242 dtk=dcvtkm(ndxtkm(ita),ita)
247 taout=(sx-sqly*cos(dtk))/ap2
255 rmax=(w2-amc2)/(1.+taout)
257 ddeler(krr,ktk)=(rmax-rcal)/nrr
258 drcurr(krr,ktk)=rcal+ddeler(krr,ktk)*(krr-0.5)
261 dsigmr(krr,ktk)=
podinl(drcurr(krr,ktk))
263 dsigmr(krr,ktk)=dsigmr(krr-1,ktk)+
podinl(drcurr(krr,ktk))
268 dsircu=
r2*dsigmr(nrr,ktk)
270 if (dsircu.le.dsigmr(krr,ktk))
then
272 drr1=dsircu/dsigmr(krr,ktk)
274 drr1=(dsircu-dsigmr(krr-1,ktk))
275 + /(dsigmr(krr,ktk)-dsigmr(krr-1,ktk))
277 drr=drcurr(krr,ktk)+(drr1-0.5d0)*ddeler(krr,ktk)
285 dom=(sx-
y)/ap/(1.+taout)
300 dphpoi=
r3*dsumph(kmp)
302 if (dphpoi.le.dsumph(kph))
then
306 temp = dsumph(kph) - dsumph(kph-1)
307 if (temp .eq. 0) temp = 1
e-20
308 dphk1 = (dphpoi - dsumph(kph-1)) / temp
309 dphk=dphi(kph)+(dphk1-1.0d0)*ddephi(kph)
316 if (r4.gt.0.5) dphk=-dphk
324 sigma_total=
sngl(tbor+tine+tpro+tnuc)
334 dgpx=dgpxy*dcos(dphig)
335 dgpy=dgpxy*dsin(dphig)
342 dgplx=-dgpz*dsts+dgpx*dcts
344 dgplz=dgpz*dcts+dgpx*dsts
345 dcphi=dcos(dble(genphi))
346 dsphi=dsin(dble(genphi))
347 dplabg(1)=dgplx*dcphi-dgply*dsphi
348 dplabg(2)=dgplx*dsphi+dgply*dcphi
353 q2tr=2d0*amh*erad*xrad*yrad
378 implicit real*8(a-
h,o-
z)
382 include
"tailcom.inc"
383 include
"amf2com.inc"
405 dimension dx(31),
dy(31),
dz(31)
410 dh=(dup-dlow)/float(ndim-1)
412 daux=dlow+dh*float(
i-1)
418 dsum2=dsum2+.5d0*(dx(
i)-dx(
i-1))*(
dy(
i)+
dy(
i-1))
422 elseif (ndim.eq.1)
then
436 + ,uncurr,plcurr,pncurr,itagen)
437 implicit real*8(a-
h,o-
z)
442 include
"tailcom.inc"
444 include
"radgenkeys.inc"
445 include
"deltacom.inc"
446 include
"density.inc"
449 real e1curr,yscurr,xscurr,uncurr,plcurr,pncurr
451 dimension tai(5),si(2,3),si0(2,3),tls(2,3,4)
460 if (abs(pncurr).eq.2.)
then
468 call
conk2(e1,xs,ys,1)
473 call
deltas(delta,delinf,tr)
511 if (ntx.eq.ntpho .and. (ita.eq.2.or.ita.eq.3))
then
525 if(ita.eq.2.and.kill_elas_res.eq.1)goto30
526 if(ita.eq.3.and.ire.le.1)
then
540 if(ita.eq.2) call
conk2(e1,xs,ys,ita)
542 if(kill_elas_res.ne.2) call
qqt(sib,tai(ita))
550 if(ita.eq.2) call
conk2(e1,xs,ys,1)
568 + +tai(1)+(tai(2)*extai2+tai(3)*extai3)/rtara
570 si(il,
in)=si(il,
in)+sig
571 si0(il,
in)=si0(il,
in)+sib
572 tls(il,
in,1)=tls(il,
in,1)+tai(1)
573 tls(il,
in,2)=tls(il,
in,2)+tai(2)*extai2/rtara
574 tls(il,
in,3)=tls(il,
in,3)+tai(3)*extai3/rtara
583 if(itagen.eq.-1.or.itagen.eq.1)
then
588 if(itagen.eq.-1.or.itagen.eq.2)
then
593 if(itagen.eq.-1.or.itagen.eq.3)
then
607 if(kill_elas_res.eq.2)sigcor=1.
655 implicit real*8(a-
h,o-
z)
665 call
conkin(e1,xs*amh/amp,ys)
679 implicit real*8(a-
h,o-
z)
709 allm=dlog((sqlm+
y)/(sqlm-
y))/sqlm
714 an=2.*alfa**2/sqls*axy*barn*amh/amp
718 dcts=(
s*sx + ap2*
y)/sqly/sqls
721 dsts=sin( acos(min(
max(dcts,-1.0),+1.0)) )
731 if ((mcset_ptarget(1:1).eq.
'L').or.
732 + (mcset_ptarget(1:2).eq.
'DT'))
then
736 elseif(mcset_ptarget(1:1).eq.
'T')
then
737 sqn=dsqrt(
s*
x*
y-aly*aml2-amp2*
y*
y)
738 ae=(-
s*
x+ap2*ym)/sqls/sqn/2.
740 ce=-(
s*
y+al2*sx)/sqls/sqn/2.
744 apn=(
y+4.*aml2)*(ae+be)+ce*sxp
745 dk2ks=as*ym+al2*bs+
cs*
x
746 dksp1=as*
s+bs*
x+
cs*ap2
747 dapks=2.*(al2*(as*ae+bs*be)+ap2*
cs*ce+ym*(as*be+bs*ae)
748 + +
s*(as*ce+
cs*ae)+
x*(bs*ce+
cs*be))
768 implicit real*8(a-
h,o-
z)
774 include
"tailcom.inc"
776 dimension sfm0(8),tm(8)
778 call
strf(0d0,0d0,sfm0)
780 tm(2)=(-(amp2*
y-
s*
x))/(2.*amp2)
781 tm(3)=(2.*(apq*dk2ks-dapks*
y)*aml)/amp
782 tm(4)=apq/amp*(-(dk2ks*sx-2.*dksp1*
y)*aml)/amp2
783 tm(7)=(-(4.*aml2+3.*apn**2-3.*apq**2+
y))/2.
784 tm(8)=apq/amp*(-3.*(apn*sxp-apq*sx))/(2.*ap)
785 ek=(3.*apq**2-
y)/amp2
789 do 1 isf=isf1,isf2,isf3
791 if(isf.eq.3.or.isf.eq.4)ppol=
pl*pn
793 if(isf.ge.5)ppol=qn/6
794 ssum=ssum+tm(isf)*sfm0(isf)*ppol
817 implicit real*8(a-
h,o-
z)
823 include
"deltacom.inc"
829 data am2/.26110d-6,.111637d-1,3.18301d0/
831 del1=-ym*(alm*allm**2/2.+2.*
fspen(2d0*sqlm/(
y+sqlm))-pi2/2.)/sqlm
832 del2=(3.*
y/2.+4.*aml2)*allm-2.
837 sqlmi=dsqrt(
y*
y+2.*a2*
y)
838 allmi=dlog((sqlmi+
y)/(sqlmi-
y))/sqlmi
839 suml=suml+2.*(
y+a2)*allmi/3.-10./9.+4.*a2*(1.-a2*allmi)/3./
y
845 elseif(
y.lt.64d0)
then
854 sumh = -(aaa+bbb*
log(1.+ccc*
y)) *2*
pi/alfa
866 deltai=aj0*dlog((w2-amc2)/aml/dsqrt(w2))
875 allss=dlog((sqlss+ss)/(-sqlss+ss))/sqlss
876 allxx=dlog((sqlxx+xx)/(-sqlxx+xx))/sqlxx
878 sfpr=dlm**2/2.-dlm*dlog(ss*xx/aml2/w2)
879 + -(dlog(ss/xx))**2/2.+
fspen((
s*
x-
y*amp2)/ss/xx)-pi2/3.
880 delta0=(ss*allss+xx*allxx)/2.+sfpr
881 delta=deltai+delta0+del1+del2+
sum
882 delinf=(dlm-1.)*dlog((w2-amc2)**2/ss/xx)
887 small_old=alfa/
pi*(pi2/6.-
fspen(1.-amp2*
y/
s/
x)
896 redfac=exp(-alfa/
pi*(dlm-1.)*
log(
s*
x/(4.*amp2*demin**2) ))
897 delta5=-(dlm-1.)*
log((w2-amc2)**2*
s*
x/(4.*amp2*demin**2*ss*xx))
917 implicit real*8(a-
h,o-
z)
924 dimension am2(3),a(5),b(5),
c(5)
925 data am2/.26110d-6,.111637d-1,3.18301d0/
926 data a/0d0,0d0,0d0,1.2227d-3,1.64178d-3/
927 data b/2.2877d-3,2.51507d-3,2.79328d-3,2.96694d-3,2.92051d-3/
928 data c/4.08041425d0,3.09624477d0,2.07463133d0,1d0,1d0/
933 sqlmi=dsqrt(
y*
y+2.*a2*
y)
934 allmi=dlog((sqlmi+
y)/(sqlmi-
y))/sqlmi
935 suml=suml+2.*(
y+a2)*allmi/3.-10./9.+4.*a2*(1.-a2*allmi)/3./
y
940 elseif(
y .lt. 16d0)
then
942 elseif(
y .lt. 100d0)
then
944 elseif(
y .lt. 8317.44d0)
then
946 elseif(
y .ge. 8317.44d0)
then
949 stop ' Y<0 in VACPOL'
967 implicit real*8(a-
h,o-
z)
991 implicit real*8(a-
h,o-
z)
1023 implicit real*8(a-
h,o-
z)
1025 include
"cmpcom.inc"
1026 include
"sxycom.inc"
1027 include
"ppicom.inc"
1028 include
"tailcom.inc"
1029 include
"radgen.inc"
1038 if(ita.eq.1 .or. ita.eq.5)
then
1039 tade=(w2-amc2)/(ap*demin) -1.d0
1040 costkm=(sx-ap2*tade)/sqly
1041 if(costkm.ge.1d0)
then
1045 dtkmax=acos(
max(-1d0,costkm ))
1049 call
intbtk2(dbtk,nbtk,dtkmax)
1053 call
inttk2(isumtk,dbtk(itk),dbtk(itk+1),dsumtk,derrtk)
1057 if(ita.le.3)ndxtkm(ita)=isumtk
1069 implicit real*8(a-
h,o-
z)
1071 include
"cmpcom.inc"
1072 include
"polcom.inc"
1073 include
"sxycom.inc"
1074 include
"ppicom.inc"
1075 include
"radgen.inc"
1076 include
"bseocom.inc"
1078 dimension tm(8,6),ajm2(2),ajm3(3),ii(8)
1079 data ii/1,2,3,4,1,2,5,6/
1082 b2=(-aly*ta+sxp*sx*ta+2.*sxp*
y)/2.
1083 b1=(-aly*ta-sxp*sx*ta-2.*sxp*
y)/2.
1084 c1=-(4.*(amp2*ta**2-sx*ta-
y)*aml2-(
s*ta+
y)**2)
1085 c2=-(4.*(amp2*ta**2-sx*ta-
y)*aml2-(ta*
x-
y)**2)
1089 bi12=(sxp*(sx*ta+2.*
y))/(sc1*sc2*(sc1+sc2))
1090 bi1pi2=1./sc2+1./sc1
1091 bis=-
b1/sc1/c1+b2/sc2/
c2
1092 bir=b2/sc2/
c2+
b1/sc1/c1
1094 b11i=(3.*
b1**2-aly*c1)/2./aly**2/sqly
1099 sqrtmb=((ta-tamin)*(tamax-ta)*(
s*
x*
y-(
y**2*amp2)-(aml2*aly)))
1100 if (sqrtmb.ge.0)
then
1103 write(*,*)
'radgen: ta,tamin,tamax,s,x,y,aml2,aly',
1104 & ta,tamin,tamax,
s,
x,
y,aml2,aly
1105 write(*,*)
"radgen: sqrtmb=",sqrtmb
1109 z1=(
y*sxp+ta*(
s*sx+ap2*
y)-ap*cos(phipoi)*sqrtmb)/aly
1110 z2=(
y*sxp+ta*(
x*sx-ap2*
y)-ap*cos(phipoi)*sqrtmb)/aly
1114 bis=bb/z2**2+bb/z1**2
1115 bir=bb/z2**2-bb/z1**2
1121 ccpe=(ae-be)*ta+2.*ce
1122 ccps=(as-bs)*ta+2.*
cs
1123 sis=(2.*bi1pi2*sps+bir*sps*ta+bis*ccps)/2.
1124 sir=( (2.*bi12*sps*ta+bir*ccps+bis*sps*ta))/2.
1125 si12=(bi12*ccps+bi1pi2*sps)/2.
1126 eis=(2.*bi1pi2*spe+bir*spe*ta+bis*ccpe)/2.
1127 eir=( (2.*bi12*spe*ta+bir*ccpe+bis*spe*ta))/2.
1128 ei12=(bi12*ccpe+bi1pi2*spe)/2.
1129 ois=((2.*bi1pi2+bir*ta)*(ccpe*sps+ccps*spe)+(ccpe*ccps+
1130 + spe*sps*ta**2)*bis+8.*bb*spe*sps+4.*bi12*spe*sps*ta**2)/4.
1131 oir=( ((2.*bi12+bis)*(ccpe*sps+ccps*spe)*ta+(ccpe*ccps+
1132 + spe*sps*ta**2)*bir+4.*bi1pi2*spe*sps*ta))/4.
1133 oi12=((ccpe*ccps+spe*sps*ta**2)*bi12+(ccpe*sps+ccps*spe)*
1134 + bi1pi2+4.*bb*spe*sps)/4.
1135 eeis=((ccpe**2+spe**2*ta**2)*bis+8.*bb*spe**2+4.*bi12*spe
1136 + **2*ta**2+4.*bi1pi2*ccpe*spe+2.*bir*ccpe*spe*ta)/4.
1137 eeir=( ((ccpe**2+spe**2*ta**2)*bir+4.*bi12*ccpe*spe*ta+4.
1138 + *bi1pi2*spe**2*ta+2.*bis*ccpe*spe*ta))/4.
1139 eei12=((ccpe**2+spe**2*ta**2)*bi12+4.*bb*spe**2+2.*bi1pi2
1141 ei1pi2=(4.*bb*spe+bi12*spe*ta**2+bi1pi2*ccpe)/2.
1142 eei1i2=((ccpe**2+spe**2*ta**2)*bi1pi2+4.*(2.*ccpe-spe*ta)
1143 + *bb*spe+8.*b1i*spe**2+2.*bi12*ccpe*spe*ta**2)/4.
1144 eb=((ccpe-spe*ta)*bb+2.*b1i*spe)/2.
1145 eeb=((ccpe-spe*ta)**2*bb+4.*(ccpe-spe*ta)*b1i*spe+4.*b11i
1147 call
ffu(1,bb,bis,bir,bi12,bi1pi2,sir,sis,si12
1148 + ,eis,eir,ei12,ei1pi2,ta)
1149 call
ffu(2,eb,eis,eir,ei12,ei1pi2,oir,ois,oi12
1150 + ,eeis,eeir,eei12,eei1i2,ta)
1151 call
ffu(3,eeb,eeis,eeir,eei12,eei1i2,0d0,0d0,0d0
1152 + ,0d0,0d0,0d0,0d0,ta)
1155 ajm3(1)=(
y-3.*apq**2)/amp2
1164 if(
i.eq.4.or.
i.eq.8)ajk=ajm2(
k)
1165 if(
i.eq.5.or.
i.eq.6)ajk=ajm3(
k)
1167 tm(
i,
j)=tm(
i,
j)+tm3(ii(
i),
j-
k+1,
k)*ajk
1168 if((
i.eq.5.or.
i.eq.6).and.
k.eq.2)
1169 + tm(
i,
j)=tm(
i,
j)+tm3(ii(
i),
j-
k+1,1)*ta/amp2
1183 subroutine ffu(n,bb,bis,bir,bi12,bi1pi2,sir,sis,si12
1184 + ,eis,eir,ei12,ei1pi2,ta)
1185 implicit real*8(a-
h,o-
z)
1187 include
"cmpcom.inc"
1188 include
"polcom.inc"
1189 include
"sxycom.inc"
1190 include
"ppicom.inc"
1191 include
"bseocom.inc"
1193 hi2=aml2*bis-ym*bi12
1194 shi2=aml2*sis-ym*si12
1195 ehi2=aml2*eis-ym*ei12
1196 ohi2=aml2*ois-ym*oi12
1199 tm3(3,1,
n)=(8.*(apq*dk2ks-dapks*
y)*aml*hi2)/amp
1200 tm3(3,2,
n)=(-2.*((2.*(bi12*dk2ks*ta-2.*shi2)*apq+(2.*shi2-
1201 + sir*
y+sis*ym)*apn+4.*dapks*hi2*ta)-4.*((2.*ei12-eis)*
1202 + dk2ks-(si12-sis)*apn)*aml2)*aml)/amp
1203 tm3(3,3,
n)=(2.*(((2.*si12+sir-sis)*apn*ta-2.*dk2ks*ei12*ta
1204 + -6.*ohi2-oir*
y+ois*ym)-4.*aml2*oi12)*aml)/amp
1205 tm3(3,4,
n)=(2.*(2.*oi12-oir+ois)*aml*ta)/amp
1206 tm3(5,1,
n)=-2.*(4.*aml2+3.*apn**2-3.*apq**2+
y)*hi2
1207 tm3(5,2,
n)=-2.*(6.*aml2*apn*eir-3.*apn**2*bi12*ta+3.*apn*
1208 + apq*bi1pi2+6.*apq*ehi2+hi2*ta)
1209 tm3(5,3,
n)=-(24.*aml2*eei12-6.*apn*ei1pi2-6.*apq*ei12*ta-
1212 tm3(4,1,
n)=(-4.*(dk2ks*sx-2.*dksp1*
y)*aml*hi2)/amp2
1213 tm3(4,2,
n)=(((2.*(sxp-2.*sx)*shi2+2.*bi12*dk2ks*sx*ta+8.*
1214 + dksp1*hi2*ta-sir*sxp*
y+sis*sxp*ym)-4.*(2.*bi12*dk2ks-bis*
1215 + dk2ks-si12*sxp+sis*sxp)*aml2)*aml)/amp2
1216 tm3(4,3,
n)=((((sxp*ta-ym)*sis-(sxp*ta-
y)*sir+2.*bi12*dk2ks
1217 + *ta+6.*shi2-2.*si12*sxp*ta)+4.*aml2*si12)*aml)/amp2
1218 tm3(4,4,
n)=(-(2.*si12-sir+sis)*aml*ta)/amp2
1219 tm3(6,1,
n)=(-3.*(apn*sxp-apq*sx)*hi2)/amp
1220 tm3(6,2,
n)=(-3.*(2.*(apn*bir+eir*sxp)*aml2-(2.*bi12*sxp*ta
1221 + -bi1pi2*sx)*apn+(bi1pi2*sxp+2.*hi2)*apq+2.*ehi2*sx))/(2.*
1223 tm3(6,3,
n)=(-3.*(8.*aml2*ei12-apn*bi1pi2-apq*bi12*ta-ei12*
1224 + sx*ta-ei1pi2*sxp))/(2.*amp)
1226 tm3(1,1,
n)=-4.*(2.*aml2-
y)*hi2
1227 tm3(1,2,
n)=4.*hi2*ta
1228 tm3(1,3,
n)=-2.*(2.*bb+bi12*ta**2)
1229 tm3(2,1,
n)=(((sxp**2-sx**2)-4.*amp2*
y)*hi2)/(2.*amp2)
1230 tm3(2,2,
n)=(2.*aml2*bir*sxp-4.*amp2*hi2*ta-bi12*sxp**2*ta+
1231 + bi1pi2*sxp*sx+2.*hi2*sx)/(2.*amp2)
1232 tm3(2,3,
n)=(2.*(2.*bb+bi12*ta**2)*amp2+4.*aml2*bi12-bi12*
1233 + sx*ta-bi1pi2*sxp)/(2.*amp2)
1248 implicit real*8(a-
h,o-
z)
1250 include
"cmpcom.inc"
1251 include
"sxycom.inc"
1252 include
"tailcom.inc"
1253 include
"ppicom.inc"
1254 include
"amf2com.inc"
1256 dimension sfm(8),iel(8)
1257 data iel/0,2,1,2,2,4,2,3/
1259 if(ita.ne.5) call
strf(taa ,
r,sfm)
1261 do 11 isf=isf1,isf2,isf3
1263 if(isf.eq.3.or.isf.eq.4)ppol=
pl*pn
1265 if(isf.ge.5)ppol=qn/6
1266 if(ita.eq.2)ppol=ppol*(amt/amp)**iel(isf)
1267 irrend=
i1(isf)+
i2(isf)-1
1268 if(ita.eq.5)irrend=1
1271 pres=-sfm0(isf)/2./
y**2/
r
1274 if(irr.eq.1.and.ita.eq.4)
pp=
pp-sfm0(isf)*(1.+
r*taa/
y)**2
1275 pres=
pp*
r**(irr-2)/(
y+
r*taa)**2/2.
1293 double precision function rv2(ta)
1297 implicit real*8(a-
h,o-
z)
1300 include
"cmpcom.inc"
1301 include
"sxycom.inc"
1302 include
"tailcom.inc"
1303 include
"amf2com.inc"
1304 include
"ppicom.inc"
1305 include
"radgen.inc"
1308 call
strf(0d0,0d0,sfm0)
1312 rmax=(w2-amc2)/(1.+ta)
1327 drcurr(irr,itkcur) = rcur
1328 ddeler(irr,itkcur) = dr
1330 dsigmr(irr,itkcur) = dr*dsum1
1335 elseif(ita.eq.2 .or. ita.eq.3)
then
1338 res=
podinl((sx-
y/aa)/(1d0+ta/aa))/(1.+ta/aa) /aa**2
1339 elseif(ita.eq.4)
then
1343 elseif(ita.eq.5)
then
1367 implicit real*8(a-
h,o-
z)
1369 include
"cmpcom.inc"
1370 include
"sxycom.inc"
1371 include
"tailcom.inc"
1372 include
"radgenkeys.inc"
1373 include
"pypars.inc"
1380 if(ita.eq.1 .or. ita.eq.4 .or. ita.eq.5)
then
1390 itara=int(rtara+0.1)
1391 itarz=int(rtarz+0.1)
1399 call
mkf2(
t,aks,itara,itarz,
f2,ff1)
1400 if (
mstp(199).eq.0)
then
1401 f1=amp*(1.+anu**2/
t)*
f2/anu/(1.+dr)
1411 if ((plrun.ne.0).and.(pnrun.ne.0))
then
1412 call
mkasym(
t,aks,itara,itarz,asym1,asym2)
1417 g1=(asym1+ga*asym2)/(1+ga**2)*f1*
fdilut(
t,aks,itara)
1418 g2=(asym2-ga*asym1)/ga/(1+ga**2)*f1*
fdilut(
t,aks,itara)
1423 elseif((plrun.eq.0).and.(abs(pnrun).eq.2))
then
1424 call
mkasym(
t,aks,itara,itarz,asym1,asym2)
1460 f1s= 0.71/(0.783**2+
t)-0.64/(1.02**2+
t)-0.13/(1.80**2+
t)
1461 f2s=-0.11/(0.783**2+
t)+0.13/(1.02**2+
t)-0.02/(1.80**2+
t)
1463 f1v= 0.05/(1.210**2+
t)-0.52/(2.45**2+
t)+0.28/(2.95**2+
t)
1464 f2v=-1.99/(1.210**2+
t)+0.20/(2.45**2+
t)+0.19/(2.95**2+
t)
1466 f1rho=(0.955+0.090/(1+
t/0.355)**2)/(1+
t/0.536)/2.
1467 f2rho=(5.335+0.962/(1+
t/0.268) )/(1+
t/0.603)/2.
1481 gen=f1n-
t*f2n/4./amp2
1488 f1=2.*amp2*
tau*gmn**2
1491 g2=2.*amp2*
tau**2*gmn*(
gen-gmn)/tau1
1496 elseif(ire.eq.1)
then
1497 f1=2.*amp2*
tau*gmp**2
1498 f2=4.*amp2*
tau*(gep**2+
tau*gmp**2)/tau1
1499 g1=amp2*
tau*2.*gmp*(gep+
tau*gmp)/tau1
1500 g2=2.*amp2*
tau**2*gmp*(gep-gmp)/tau1
1505 elseif (ire.eq.2)
then
1506 call
ffdeu(
t,fcdeu,fmdeu,fqdeu)
1510 f1=4./3.*
tau*tau1*fm**2
1513 g2=-1./3.*
tau**2*fm*(2.*
tau*fq+6*
fc-3.*fm)
1516 + 16./3.*
tau**3/tau1*(
tau*fq+3.*
fc-3.*fm)*fq
1517 b3=-4./3.*(3.*
tau+2.)*fm**2*
tau**2+
1518 + 16./9.*
tau**3/tau1*(
tau*fq+3.*
fc-3.*fm)*fq
1519 b4=4./3.*(6.*
tau+1.)*fm**2*
tau**2+
1520 + 16./9.*
tau**3/tau1*(
tau*fq+3.*
fc+3.*(3.*
tau+2.)*fm)*fq
1521 elseif (ire.eq.3)
then
1523 f1=2.*amt**2*
tau*gm**2
1524 f2=4.*amt**2*
tau*(ge**2+
tau*gm**2)/tau1
1525 g1=amt**2*
tau*2.*gm*(ge+
tau*gm)/tau1
1526 g2=2.*amt**2*
tau**2*gm*(ge-gm)/tau1
1531 elseif(ire.eq.4)
then
1535 f2=4.*amp2*
tau*(rtarz*ge)**2
1542 elseif((ire.eq.14).or.(ire.eq.84).or.(ire.eq.20))
then
1551 f2=4.*amp2*
tau*(rtarz*ff)**2
1559 elseif (ita.eq.3)
then
1562 call
ffquas(
t,geun,gmun,gepo,gmpo)
1563 f1=2.*amp2*
tau*gmun**2
1564 f2=4.*amp2*
tau*(geun**2+
tau*gmun**2)/tau1
1565 g1=amp2*
tau*2.*gmpo*(gepo+
tau*gmpo)/tau1
1566 g2=2.*amp2*
tau**2*gmpo*(gepo-gmpo)/tau1
1575 sfm(2)=epsi*(
f2+qn/6.*b2)
1579 sfm(6)=epsi**3*(b2/3.+b3+b4)
1580 sfm(7)=epsi*(b2/3.-b3)
1581 sfm(8)=epsi**2*(b2/3.-b4)
1593 implicit real*8(a-
h,o-
z)
1595 include
"cmpcom.inc"
1597 gep=1.2742/(1.+
t/0.6394**2)-.2742/(1.+
t/1.582**2)
1598 gmp=(1.3262/(1.+
t/0.6397**2)-.3262/(1.+
t/1.3137**2))*amm
1611 implicit real*8(a-
h,o-
z)
1613 parameter(c2i3 = 2/3)
1614 dimension a(4),b(4),
c(4)
1615 dimension al2ar(4),be2ar(4),ga2ar(4)
1616 data amd/1.8756280d0/chbar/0.19732858d0/
1617 data amp/0.9382796d0/
1618 data dmu/0.857406d0/dqu/25.84d0/
1624 gd2=1d0/(1.d0+
t/4./0.89852**2)**4
1626 gd2e=gd2/(2d0*
eta+1d0)
1630 al2ar(1)=1.8591*chbar**2
1631 al2ar(4)=2d0*amd*0.58327d0
1632 be2ar(1)=19.586*chbar**2
1633 be2ar(4)=2d0*amd*0.1d0
1634 ga2ar(1)=1.0203*chbar**2
1635 ga2ar(4)=2d0*amd*0.17338d0
1637 al2ar(1)=1.5250*chbar**2
1638 al2ar(4)=2d0*amd*0.24086d0
1639 be2ar(1)=43.677*chbar**2
1640 be2ar(4)=2d0*amd*0.029138d0
1641 ga2ar(1)=1.8706*chbar**2
1642 ga2ar(4)=2d0*amd*0.42693d0
1646 al2ar(
i)=al2ar(1)+(al2ar(4)-al2ar(1))/3d0*(
i-1)
1647 be2ar(
i)=be2ar(1)+(be2ar(4)-be2ar(1))/3d0*(
i-1)
1648 ga2ar(
i)=ga2ar(1)+(ga2ar(4)-ga2ar(1))/3d0*(
i-1)
1652 a(1)=2.4818*chbar**2
1653 a(2)=-10.850*chbar**2
1654 a(3)=6.4416*chbar**2
1656 a(1)=1.5706*chbar**2
1657 a(2)=12.238*chbar**2
1658 a(3)=-42.046*chbar**2
1660 a(4)=al2ar(4)*(1d0-a(2)/al2ar(2)-a(3)/al2ar(3)-a(1)/al2ar(1))
1670 bzn=1d0/be2ar(4)-1d0/be2ar(3)
1671 bbb=(2d0-dmu*amd/amp)/2d0/sqrt(2d0)/amd
1672 b(3)=(b(1)/be2ar(1)+b(2)/be2ar(2)-b(1)/be2ar(4)-b(2)/be2ar(4)
1674 b(4)=-(b(1)/be2ar(1)+b(2)/be2ar(2)-b(1)/be2ar(3)-b(2)/be2ar(3)
1676 ccc=(1d0-dmu*amd/amp-dqu)/4./amd**2
1684 znc2=ga2ar(1)*(ga2ar(3)*ga2ar(4)-ga2ar(2)*ga2ar(3)
1685 + +ga2ar(2)**2-ga2ar(2)*ga2ar(4))
1686 c(2)=-ga2ar(2)/znc2*(
c(1)*(
1687 + ga2ar(3)*ga2ar(4)-ga2ar(1)*ga2ar(3)+ga2ar(1)**2
1688 + -ga2ar(1)*ga2ar(4))-ccc*ga2ar(3)*ga2ar(4)*ga2ar(1) )
1689 znc3=ga2ar(1)*(ga2ar(3)-ga2ar(2))*(ga2ar(4)-ga2ar(3))
1690 c(3)=ga2ar(3)/znc3*(
c(1)*(
1691 + ga2ar(2)*ga2ar(4)-ga2ar(1)*ga2ar(4)+ga2ar(1)**2
1692 + -ga2ar(1)*ga2ar(2))-ccc*ga2ar(2)*ga2ar(4)*ga2ar(1) )
1693 znc4=ga2ar(1)*(ga2ar(4)-ga2ar(2))*(ga2ar(4)-ga2ar(3))
1694 c(4)=-ga2ar(4)/znc4*(
c(1)*(
1695 +ga2ar(2)*ga2ar(3)-ga2ar(1)*ga2ar(3)+ga2ar(1)**2-ga2ar(1)*ga2ar(2)
1696 + )-ccc*ga2ar(2)*ga2ar(3)*ga2ar(1) )
1703 g00=g00+a(
i)/(al2ar(
i)+
t)
1704 gp0=gp0+sqt*b(
i)/(be2ar(
i)+
t)
1705 gpm=gpm+
t*
c(
i)/(ga2ar(
i)+
t)
1708 gc=gd2e*( (1d0-c2i3*
eta)*g00+4d0*c2i3*sq2e*gp0
1709 + +c2i3*(2d0*
eta-1d0)*gpm)
1710 gm=gd2e*(2d0*g00+2d0*(2d0*
eta-1d0)/sq2e*gp0-2d0*gpm)
1711 gq=gd2e*(-g00+2d0/sq2e*gp0-(1d0+1d0/
eta)*gpm)
1725 implicit real*8(a-
h,o-
z)
1727 include
"cmpcom.inc"
1740 f0=
ddexp(-a**2*qf**2) - b**2*qf**2*
ddexp(-
c**2*qf**2)
1741 fm=
ddexp(-am**2*qf**2) - bm**2*qf**2*
ddexp(-cm**2*qf**2)
1742 df=d*
ddexp(-((qf-q0)/
p)**2)
1744 gm=fm*rtara/rtarz * (-2.13)
1757 implicit real*8(a-
h,o-
z)
1759 include
"cmpcom.inc"
1761 dimension dr(20),dq(20)
1762 DATA dr/.2d0,.6d0,.9d0,1.4d0,1.9d0,2.3d0,2.6d0
1763 : ,3.1d0,3.5d0,4.2d0,4.9d0,9*0.d0/
1764 DATA dq/.034724d0,.430761d0,.203166d0,.192986d0,.083866d0
1765 : ,.033007d0,.014201d0,0.d0,.00686d0,0.d0,.000438d0,9*0.d0/
1766 DATA dgam2/.666666666667d0/
1775 d2rg=2.d0*dr(
i)**2/dgam2
1778 IF(dqri.GE.1.d-6) dj0=dsin(dqri)/dqri
1779 dai=dq(
i)/(1.d0+d2rg)*(dcos(dqri)+d2rg*dj0)
1782 delast=dexp(-dq2f*dgam2/4.d0) * dabs(delast)
1793 include
"tailcom.inc"
1799 iq=max0(int(q2/q2bin)+1,1)
1800 IF(iq+1.LT.nqbin)
forgetp=ffnuc(iq)+
1801 1 (ffnuc(iq+1)-ffnuc(iq))*(q2/q2bin-float(iq-1))
1812 include
"tailcom.inc"
1813 real q2max,rmin,rmax,
eps,hbc
1819 DATA q2max,rmin,rmax,
eps,nqbin1/.3,1.
e-6,20.,1.
e-6,600/
1825 q2bin=q2max/float(nqbin-1)
1829 q2=q2bin*float(iq-1)
1831 qfor=sqrt(q2)*1000./hbc
1835 IF(iq.GT.1)ffnuc(iq)=ffnuc(iq)/ffnuc(1)
1851 include
"tailcom.inc"
1854 IF (qfor.GT.1.
e-5)
THEN
1873 include
"tailcom.inc"
1874 include
"cmpcom.inc"
1875 real r,zpara,cpara,wpara
1884 1 (1.0+exp((
r-cpara)/zpara))
1885 elseif (ire.eq.20)
then
1892 1 (1.0+exp((
r-cpara)/zpara))
1893 elseif (ire.eq.84)
then
1899 cpara=1.12*rtara**(1./3.)-0.86/rtara**(1./3.)
1900 forrhop=1./(1.+exp((
r-cpara)/zpara))
1907 cpara=1.07*rtara**(1./3.)
1908 forrhop=1./(1.+exp((
r-cpara)/zpara))
1921 implicit real*8(a-
h,o-
z)
1923 include
"cmpcom.inc"
1924 include
"tailcom.inc"
1936 if(
t.le.(2.d0*fermom)**2)
then
1937 sqrat=dsqrt(
t)/fermom
1938 supele=0.75*sqrat-sqrat**3/16.
1947 geun=gep*dsqrt(supele*rtarz)
1949 gmun=gep*dsqrt(supmag*(rtarz*amm**2+rtarn*amn**2))
1956 gmpo=gep*dsqrt(supmag*(rtarn*amn**2))
1963 implicit real*8(a-
h,o-
z)
1964 data chbar/.197328d0/
1973 delff=(datan(sqtf/.93d0)-2.*datan(sqtf/3.19d0)+
1974 + datan(sqtf/5.45d0))*1.580d0/sqtf
1975 delff=dmax1(0.d0,delff)
1977 supst=dmax1(0.d0,supele)
1988 implicit real*8(a-
h,o-
z)
2006 double precision xl,xu,
y,a,b,
c,fct
2010 c=.49863193092474078d0*b
2011 y=.35093050047350483d-2*(fct(a+
c)+fct(a-
c))
2012 c=.49280575577263417d0*b
2013 y=
y+.8137197365452835d-2*(fct(a+
c)+fct(a-
c))
2014 c=.48238112779375322d0*b
2015 y=
y+.12696032654631030d-1*(fct(a+
c)+fct(a-
c))
2016 c=.46745303796886984d0*b
2017 y=
y+.17136931456510717d-1*(fct(a+
c)+fct(a-
c))
2018 c=.44816057788302606d0*b
2019 y=
y+.21417949011113340d-1*(fct(a+
c)+fct(a-
c))
2020 c=.42468380686628499d0*b
2021 y=
y+.25499029631188088d-1*(fct(a+
c)+fct(a-
c))
2022 c=.39724189798397120d0*b
2023 y=
y+.29342046739267774d-1*(fct(a+
c)+fct(a-
c))
2024 c=.36609105937014484d0*b
2025 y=
y+.32911111388180923d-1*(fct(a+
c)+fct(a-
c))
2026 c=.33152213346510760d0*b
2027 y=
y+.36172897054424253d-1*(fct(a+
c)+fct(a-
c))
2028 c=.29385787862038116d0*b
2029 y=
y+.39096947893535153d-1*(fct(a+
c)+fct(a-
c))
2030 c=.25344995446611470d0*b
2031 y=
y+.41655962113473378d-1*(fct(a+
c)+fct(a-
c))
2032 c=.21067563806531767d0*b
2033 y=
y+.43826046502201906d-1*(fct(a+
c)+fct(a-
c))
2034 c=.16593430114106382d0*b
2035 y=
y+.45586939347881942d-1*(fct(a+
c)+fct(a-
c))
2036 c=.11964368112606854d0*b
2037 y=
y+.46922199540402283d-1*(fct(a+
c)+fct(a-
c))
2038 c=.7223598079139825d-1*b
2039 y=
y+.47819360039637430d-1*(fct(a+
c)+fct(a-
c))
2040 c=.24153832843869158d-1*b
2041 y=b*(
y+.48270044257363900d-1*(fct(a+
c)+fct(a-
c)))
2053 implicit real*8 (a-
h,o-
z)
2058 include
"sxycom.inc"
2059 include
"cmpcom.inc"
2060 include
"ppicom.inc"
2061 include
"radgen.inc"
2068 dwids=dsqrt(ap*aml/
s)/dc
2069 dwidp=dsqrt(ap*aml/
x)/dc
2073 dts=acos( min(1.0d0,(
s*sx+ap2*
y)/sqls/sqly) )
2074 dtp=acos( min(1.0d0,(
x*sx-ap2*
y)/sqlx/sqly) )
2081 dbtk(2)=dmin1(4.d0*dc*dwids,dts/3.0d0)
2082 dbtk(3)=dmax1(dts-dwids,dts/1.5d0)
2083 dbtk(4)=dmin1(dts+dwids,dts+(dtp-dts)/3.0d0)
2084 dbtk(5)=dmax1(dtp-dwidp,dts+(dtp-dts)/1.5d0)
2085 dbtk(6)=dmin1(dtp+dwidp,dtp+(dpi-dtp)/3.0d0)
2086 dbtk(7)=.98d0*dbtk(6)+.02d0*dpi
2094 if(dbtk(nbtk+1).ge.dtkmax) goto 20
2097 20 dbtk(nbtk+1)=dmin1(dbtk(nbtk+1),dtkmax)
2111 subroutine inttk2(isumtk,dbtk1,dbtk2,dsumtk,derrtk)
2112 implicit real*8 (a-
h,o-
z)
2115 include
"radgen.inc"
2116 include
"sxycom.inc"
2117 include
"cmpcom.inc"
2118 include
"tailcom.inc"
2119 include
"ppicom.inc"
2121 ddtk=(dbtk2-dbtk1)/ntk
2133 ta=(sx-sqly*cos(dtk))/ap2
2134 dsigma=an*alfa/
pi*
rv2(ta) * sin(dtk)*sqly/ap2
2138 dsitkm(isumtk,ita) = dsumtk + ddtk * dsum
2139 dcvtkm(isumtk,ita) = dtk
2140 ddetkm(isumtk,ita) = ddtk
2144 if (itk.gt.2.and.itk.lt.ntk)
then
2145 derr1=derr1+dabs(dsigma-dold)
2146 elseif(itk.eq.2.or.itk.eq.ntk)
then
2147 derr1=derr1+1.5d0*dabs(dsigma-dold)
2156 dsumtk=dsumtk+ddtk*dsum
2157 derrtk=derrtk+ddtk*derr1
2169 subroutine itafun(ys,xs,pl,pn,ixytest,itagen)
2171 include
"cmpcom.inc"
2172 include
"radgen.inc"
2173 include
"density.inc"
2174 include
"xytabcom.inc"
2176 dimension xym(2),na(2),a(2*
nt)
2180 if (ys.gt.
y(nny)) goto 10
2182 5
print *,
' ys=',ys,
' out of the region'
2184 10
if (nny.eq.
nt) goto 5
2186 if(ixytest.eq.1)
then
2187 stop .eq.
' ixytest1 is not supported in the version'
2202 if (ntx.eq.ntdis)
then
2203 tbor=fint(2,xym,na,a,tbor_udis)
2204 sigrad=fint(2,xym,na,a,sigrad_udis)
2205 sig1g=fint(2,xym,na,a,sig1g_udis)
2206 vac=fint(2,xym,na,a,vac_udis)
2207 vertex=fint(2,xym,na,a,vertex_udis)
2208 small=fint(2,xym,na,a,small_udis)
2209 redfac=fint(2,xym,na,a,redfac_udis)
2211 tbor=fint(2,xym,na,a,tbor_upho)
2212 sigrad=fint(2,xym,na,a,sigrad_upho)
2213 sig1g=fint(2,xym,na,a,sig1g_upho)
2214 vac=fint(2,xym,na,a,vac_upho)
2215 vertex=fint(2,xym,na,a,vertex_upho)
2216 small=fint(2,xym,na,a,small_upho)
2217 redfac=fint(2,xym,na,a,redfac_upho)
2233 if (ntx.eq.ntdis)
then
2234 tine=fint(2,xym,na,a,tine_udis)
2235 tnuc=fint(2,xym,na,a,tnuc_udis)
2236 tpro=fint(2,xym,na,a,tpro_udis)
2238 tine=fint(2,xym,na,a,tine_upho)
2239 tnuc=fint(2,xym,na,a,tnuc_upho)
2240 tpro=fint(2,xym,na,a,tpro_upho)
2243 if(rr1.gt.sigrad-tbor)
then
2251 if(rr1.gt.(tpro+tnuc))
then
2255 elseif(rr1.gt.tnuc)
then
2279 if (ntx.le.ntdis)
then
2280 dsitkm(ii,itagen)=fint(2,xym,na,a,densdis(1,1,ii,itagen))
2282 dsitkm(ii,itagen)=fint(2,xym,na,a,denspho(1,1,ii,itagen))
2300 if (ntx.le.ntdis)
then
2301 ddetkm(ii,itagen)=fint(2,xym,na,a,widdis(1,1,ii,itagen))
2303 ddetkm(ii,itagen)=fint(2,xym,na,a,widpho(1,1,ii,itagen))
2310 do itoo=1,iittkk/ntk
2311 ssuumm=ssuumm+ddetkm(itoo,itagen)*ntk
2313 dcvtkm(iittkk,itagen)=ssuumm+ddetkm((iittkk-1)/ntk+1,itagen)
2314 + *(mod(iittkk,ntk)-0.5)
2327 subroutine xytabl(tname,e1,plrun,pnrun,ixytest,ire)
2329 include
"cmpcom.inc"
2330 include
"radgen.inc"
2331 include
"density.inc"
2332 include
"xytabcom.inc"
2333 include
"mc_set.inc"
2338 open (unit=lun,
file=tname)
2340 if(ixytest.eq.0)
then
2343 y(iy)=radgen_ymin+(radgen_ymax-radgen_ymin)/(nty-1)*(iy-1)
2353 elseif (ix.le.27)
then
2356 elseif (ix.le.36)
then
2359 elseif (ix.le.48)
then
2360 if (mod(ix,2).gt.0)
then
2374 ylim=(bmc2-bmp**2)/(2.*bmp*e1*(1d0-
x(ix)))
2377 call
mpolrad(e1,ys,
x(ix),1.,plrun,pnrun,-1)
2379 sig1g_u(ix,iy)=sig1g
2380 sigrad_u(ix,iy)=sigrad
2386 small_u(ix,iy)=
small
2387 redfac_u(ix,iy)=redfac
2392 write(lun,
'(1x,14g14.4)')
x(ix),
y(iy)
2393 + ,sig1g_u(ix,iy),sigrad_u(ix,iy),tbor_u(ix,iy)
2394 + ,tine_u(ix,iy),tnuc_u(ix,iy),tpro_u(ix,iy)
2395 + ,vac_u(ix,iy),vertex_u(ix,iy),small_u(ix,iy)
2399 width(ix,iy,itkm+1,1)=ddetkm(itkm*ntk+1,1)
2400 width(ix,iy,itkm+1,2)=ddetkm(itkm*ntk+1,2)
2402 +
width(ix,iy,itkm+1,3)=ddetkm(itkm*ntk+1,3)
2406 if(ire.eq.1)ittmax=2
2407 write(lun,*)((
width(ix,iy,itkm+1,itt)
2408 + ,itkm=0,6),itt=1,ittmax),ndxtkm
2416 if (dsitkm(itkm,1).lt.1.
e-24 .or.
2417 & dsitkm(ntk*7,1).lt.1
e-24 )
then
2418 denstk(ix,iy,itkm,1)=1.0
2420 denstk(ix,iy,itkm,1)=dsitkm(itkm,1)/dsitkm(ntk*7,1)
2423 if (dsitkm(itkm,2).lt.1.
e-24 .or.
2424 & dsitkm(ntk*7,2).lt.1
e-24 )
then
2425 denstk(ix,iy,itkm,2)=1.0
2427 denstk(ix,iy,itkm,2)=dsitkm(itkm,2)/dsitkm(ntk*7,2)
2438 if (dsitkm(itkm,3).lt.1.
e-24 .or.
2439 & dsitkm(ntk*7,3).lt.1
e-24 )
then
2440 denstk(ix,iy,itkm,3)=1.0
2442 denstk(ix,iy,itkm,3)=dsitkm(itkm,3)/dsitkm(ntk*7,3)
2446 write(lun,
'(35g14.4)')(denstk(ix,iy,itkm,1)
2448 write(lun,
'(35g14.4)')(denstk(ix,iy,itkm,2)
2451 +
write(lun,
'(35g14.4)')(denstk(ix,iy,itkm,3)
2457 elseif(ixytest.gt.0)
then
2463 if (ntx.eq.ntdis)
then
2464 read(lun,*)
x(ix),
y(iy)
2465 + ,sig1g_udis(ix,iy),sigrad_udis(ix,iy),tbor_udis(ix,iy)
2466 + ,tine_udis(ix,iy),tnuc_udis(ix,iy),tpro_udis(ix,iy)
2467 + ,vac_udis(ix,iy),vertex_udis(ix,iy)
2468 + ,small_udis(ix,iy),redfac_udis(ix,iy)
2474 if(ire.eq.1)ittmax=2
2476 + ((widdis(ix,iy,ite,itt),ite=1,7),itt=1,ittmax),ndxtkm
2477 read(lun,*)(densdis(ix,iy,itkm,1)
2479 read(lun,*)(densdis(ix,iy,itkm,2)
2482 +
read(lun,*)(densdis(ix,iy,itkm,3)
2485 read(lun,*)
x(ix),
y(iy)
2486 + ,sig1g_upho(ix,iy),sigrad_upho(ix,iy),tbor_upho(ix,iy)
2487 + ,tine_upho(ix,iy),tnuc_upho(ix,iy),tpro_upho(ix,iy)
2488 + ,vac_upho(ix,iy),vertex_upho(ix,iy)
2489 + ,small_upho(ix,iy),redfac_upho(ix,iy)
2492 if(ire.eq.1)ittmax=2
2494 + ((widpho(ix,iy,ite,itt),ite=1,7),itt=1,ittmax),ndxtkm
2495 read(lun,*)(denspho(ix,iy,itkm,1)
2497 read(lun,*)(denspho(ix,iy,itkm,2)
2500 +
read(lun,*)(denspho(ix,iy,itkm,3)
2507 stop 'xytabl -- ixytest < 0'