Numerical Calculation of the rate of the radiative pion decay
p
+->e+ n e
g
plot_pienug_x.kumac
Parameter:
1 g
2 Threshold Energy
3 Threshold Angle
macro plot_pienug_x 1=0. 2=5. 3=55
* Macro to plot the terms, contributing to the radiative decay
* Formula from A.Stetz
* factorization after Bryman, units 10e-6
hi/del *
ve/del *
Application comis quit
subroutine calc_terms(gamma,x_thresh,Y_THRESH)
IMPLICIT NONE
int i,j,k,n_max,s,t
real x,y,l,lx,sum_y,y_thresh,gamma,y_energ,x_thresh,r_all
real hx,hy,xstart,ystart
REAL SD_p(30,30),SD_n(30,30),IB(30,30),INT_p(30,30),INT_n(30,30)
vector rate_IB(1),rate_SD_P(1),rate_SD_N(1),rate_INT(1)
VECTOR rINT_SD_p(30),rINT_SD_n(30),rINT_IB(30),
& rINT_INT(30)
data N_MAX,s,t /30,0,0/
xstart=x_thresh/139.57*2.
ystart=y_thresh/139.57*2.
hx=(1-xstart)/(n_max*1.-2.)
hy=(1-ystart)/(n_max*1.-2.)
print *,gamma,hx,hy
do i=1,N_MAX-1
x=xstart+(i-1)*hx
do j=1,N_MAX-1
y=ystart+(j-1)*hy
l=x+y-1.
lx=(y-1.)/x
if (l.gt.0.) then
IB(i,j)=(1.-y)/x/x*((x-1.)*(x-1.)+1.)/l
SD_p(i,j)=(1.-x)*l*l*(1.+gamma)**2
SD_n(i,j)=(1.-x)*(1.-y)*(1.-y)*(1.-gamma)**2
INT_p(i,j)=lx*(1.-x)*(1.+gamma)
INT_n(i,j)=lx*(x-1.-x*x/l)/x*(1.-gamma)
endif
enddo
enddo
rate_IB(1)=0.
rate_sd_p(1)=0.
rate_sd_n(1)=0.
rate_int(1)=0.
DO i=1,N_MAX-1 ! threshhold
rINT_IB(i)=0.
rINT_SD_p(i)=0.
rINT_SD_n(i)=0.
RINT_INT(I)=0.
do j=1,n_max-1
if (s.ne.4) then
s=4
else
s=2
endif
if (j.eq.1.or.j.eq.n_max-1) then
s=1
endif
rINT_IB(i)=rint_IB(i)+s*ib(i,j)
rINT_SD_p(i)=rint_SD_p(i)+s*SD_p(i,j)
rINT_SD_n(i)=rint_SD_n(i)+s*SD_n(i,j)
rint_int(I)=rint_INT(i)+s*INT_n(i,j)+s*INT_P(I,J)
enddo
if (t.ne.4) then
t=4
else
t=2
endif
if (i.eq.1.or.i.eq.n_max-1) then
t=1
endif
rate_IB(1)=rate_IB(1)+t*rint_IB(i)*hy/3.
rate_SD_P(1)=rate_SD_P(1)+t*rint_SD_P(i)*hy/3.
rate_SD_N(1)=rate_SD_N(1)+t*rint_SD_N(i)*hy/3.
rate_INT(1)=rate_INT(1)+t*rint_INT(i)*hy/3.
enddo
rate_IB(1)=.1436*rate_IB(1)*hx/3.*1.e-6
rate_SD_P(1)=.2106*rate_SD_P(1)*hx/3.*1.e-5
rate_SD_N(1)=.2106*rate_SD_N(1)*hx/3.*1.e-5
rate_INT(1)=.004*rate_INT(1)*hx/3.*1.e-6
r_all=rate_IB(1)+rate_SD_P(1)+rate_SD_N(1)+rate_int(1)
print *,r_all,rate_IB,rate_SD_P,rate_SD_N,rate_int
do i=1,30
rINT_IB(i)=.1436*rint_IB(i)*1.e-6*hy/3.
rINT_SD_p(i)=2.106*rint_SD_P(i)*1.e-6*hy/3.
rINT_SD_n(i)=2.106*rint_SD_n(i)*1.e-6*hy/3.
* rint_INT_p(i)=0.004*rint_INT_p(i)*hy/3.
* rint_int_n(i)=0.004*rint_INT_n(i)
rint_INT(I)=0.004*rint_INT(I)*1.e-6*hy/3.
enddo
end
quit
call calc_terms($RSIGMA([1]),$RSIGMA([2]),$RSIGMA([3]))
ve/cr gam_en(29) R
do i=1,29
rtemp=[2]+[i]/29.*(139.57/2.-[2])
ve/in gam_en([i]) [rtemp]
enddo
t= Limits [1], [2]MeV, [3]MeV
null [2] 70. 1.e-10 2.e-6
atitle E[g] 'Decay Rate'
itx 10 3.e-6 $quote([t])