diff --git a/bedepe.f b/bedepe.f
new file mode 100644
index 0000000000000000000000000000000000000000..d3da260f84ab469c64c8da8d46f9019219106591
--- /dev/null
+++ b/bedepe.f
@@ -0,0 +1,437 @@
+ program BEDEPE
+
+* Oliver S. Kirsebom
+* 31 Nov 2010
+
+* Differs from the 0th version in that
+* FUNLXP and FUNLUX are used to generate
+* beta energies according to phase-space
+* distribution (the 0th version uses a flat
+* distribution from 0 to T0, and includes the
+* phase-space factor in the weight)
+
+
+* * * * * * * * * * *
+* DECLARATIONS *
+* * * * * * * * * * *
+
+ REAL pi
+ PARAMETER (pi = 3.1415927)
+ REAL*8 ran2
+ EXTERNAL triple,DummyIni,phasespace,ran2
+ character*80 histout,dataout
+ COMMON /PAWC/HMEMOR(8000000)
+ real malpha,me,QEC,mn,mp,mpar
+ parameter(mn = 939.565002441406)
+ parameter(mp = 938.271972656250)
+ parameter(malpha = 3727378.97610664)
+ parameter(amu = 931494.0)
+ PARAMETER(me = 510.99906)
+ REAL pacm(3),pbeta(3),pnu(3),prec(3)
+ + ,vacm(3),vfcm(3),vrec(3),valab(3),vflab(3)
+ + ,palab(3),pflab(3)
+ logical check
+ real T0,E0,va
+ common /COEFF/xi,Theta,a3
+ PARAMETER(Nmont=200,Nbig=1E7)
+ REAL Xbeta(Nbig),Ybeta(Nmont)
+ COMMON /PHSP/T0me
+
+
+
+* READ FROM INPUT FILE
+ open(91,file='input0',status='old')
+ read(91,*) Nevent
+ read(91,*) Qbeta
+ read(91,*) alpha_thres
+ read(91,*) Eexc
+ read(91,*) Ai
+ read(91,*) npar
+ read(91,*) a3
+ read(91,*) xi
+ read(91,*) Theta
+ close(91)
+
+
+* Check
+ if (Nevent.gt.Nbig) then
+ write(6,*) 'Number of events is too large'
+ write(6,*) 'Parameter NBIG must be increased (see code)'
+ STOP
+ endif
+
+
+
+* PARTICLE EMITTED
+ if (npar.eq.0) then !neutron
+ mpar = mn
+ Af = Ai-1.
+ elseif (npar.eq.1) then !proton
+ mpar = mp
+ Af = Ai-1.
+ elseif (npar.eq.2) then !alpha
+ mpar = malpha
+ Af = Ai-4.
+ endif
+
+* Comment: For historical reasons (the program was first
+* written only for alpha decay) many of the variables are
+* named alpha-something though they are generic, i.e. they
+* work for neutrons and protons as well.
+
+
+
+* INITIATE
+
+ Call DummyIni
+ call hlimit(8000000)
+ histout = 'bedepe.his'
+ dataout = 'bedepe.dat'
+ check = .false.
+
+
+
+* BOOK HISTOGRAMS
+
+ nhis_rec = 10
+ call hbook1(nhis_rec,'Recoil energy (keV)',250,0.,50.,0.)
+ nhis_eb = 11
+ call hbook1(nhis_eb,'beta kinetic energy (keV)',400,0.,20000.,0.)
+ nhis_ea = 12
+ call hbook1(nhis_ea,'alpha lab energy (keV)',400,0.,20000.,0.)
+ nhis_ef = 13
+ call hbook1(nhis_ef,'final nucleus lab energy (keV)'
+ + ,400,0.,20000.,0.)
+ nhis_ba = 14
+ call hbook1(nhis_ba,'lab angle between beta and '
+ + //'alpha (degree)',360,0.,180.,0.)
+ nhis_af = 15
+ call hbook1(nhis_af,'lab angle between alpha and '
+ + //'final nucleus (degree)',120,120.,180.,0.)
+
+
+
+*::::::::::::::::::::::::::::::::::::::::*
+* SIMULATION OF B8 BETA DECAY *
+*::::::::::::::::::::::::::::::::::::::::*
+*
+* Take direction of alpha in recoil frame (c.m.) as z-axis
+ Qbe8 = Eexc - alpha_thres
+ Ealpha = Af/Ai*Qbe8
+ rpacm = sqrt(2.*mpar*Ealpha)
+ pacm(1) = 0.
+ pacm(2) = 0.
+ pacm(3) = rpacm
+
+
+*=======================================================================*
+* Sample beta energy (E_beta)
+*=======================================================================*
+
+ T0 = Qbeta - Eexc ! end point of beta kinetic energy spectrum
+ T0me = T0/me
+ if (T0.lt.0.) then
+ write(6,*) 'ERROR : excitation energy exceeds maximum value'
+ STOP
+ endif
+ call funlxp(phasespace,Ybeta,0.,T0me)
+ call funlux(Ybeta,Xbeta,Nevent)
+
+
+ Wmax = 0.
+ sum = 0.
+ Wsum = 0.
+ Esum = 0.
+ ieve = 1
+ angle_min = 180.
+ open(41,file=dataout,status='unknown')
+
+ do while(ieve.le.Nevent)
+
+ IF (MOD(ieve,10000).EQ.0) THEN
+ WRITE(6,*) 'ieve ... ',ieve
+ + ,'(',INT(100.*real(ieve)/real(Nevent)),'%)'
+ ENDIF
+
+
+*=======================================================================*
+* Sample beta energy (E_beta) and direction (theta_beta,phi_beta)
+* and neutrino direction (theta_nu,phi_nu)
+*=======================================================================*
+
+* Beta energy
+ T_beta = Xbeta(ieve)*me
+ T_beta = max(0.,T_beta)
+ E_beta = me + T_beta
+
+* Neutrino energy
+ E_nu = T0 - T_beta
+
+* Beta direction
+ phi_beta = 2.*pi*ran2(10)
+ costh_beta = 2.*ran2(10)-1.
+ sinth_beta = sqrt(1.-costh_beta**2.)
+
+* Neutrino direction
+ phi_nu = 2.*pi*ran2(10)
+ costh_nu = 2.*ran2(10)-1.
+ sinth_nu = sqrt(1.-costh_nu**2.)
+
+* Beta momentum
+ rpbeta = sqrt( E_beta**2. - me**2. )
+ pbeta(1) = sinth_beta * cos(phi_beta) * rpbeta
+ pbeta(2) = sinth_beta * sin(phi_beta) * rpbeta
+ pbeta(3) = costh_beta * rpbeta
+
+* Neutrino momentum
+ rpnu = E_nu
+ pnu(1) = sinth_nu * cos(phi_nu) * rpnu
+ pnu(2) = sinth_nu * sin(phi_nu) * rpnu
+ pnu(3) = costh_nu * rpnu
+
+
+*=======================================================================*
+* Calculate recoil energy of daughter nucleus
+*=======================================================================*
+
+* Momentum of recoil
+ do i=1,3
+ prec(i) = -pbeta(i)-pnu(i)
+ enddo
+ precsq = prec(1)**2.+prec(2)**2.+prec(3)**2.
+
+* Energy of recoil
+ Erec = precsq/(2.*Ai*amu)
+
+
+*=======================================================================*
+* Calculate velocities of decay products in lab
+*=======================================================================*
+
+* Velocity of alpha and final nucleus in recoil frame
+ do i=1,3
+ vacm(i) = pacm(i)/mpar
+ vfcm(i) = -pacm(i)/(Af*amu)
+ enddo
+
+* Velocity of recoil
+ do i=1,3
+ vrec(i) = prec(i)/(Ai*amu)
+ enddo
+
+* Velocity and momenta of alphas in lab
+ do i=1,3
+ valab(i) = vacm(i) + vrec(i)
+ palab(i) = mpar*valab(i)
+ vflab(i) = vfcm(i) + vrec(i)
+ pflab(i) = Af*amu*vflab(i)
+ enddo
+ rpalab = sqrt(palab(1)**2.+palab(2)**2.+palab(3)**2.)
+ rpflab = sqrt(pflab(1)**2.+pflab(2)**2.+pflab(3)**2.)
+
+* Energies of alphas in lab
+ valabsq = valab(1)**2. + valab(2)**2. + valab(3)**2.
+ vflabsq = vflab(1)**2. + vflab(2)**2. + vflab(3)**2.
+ Ealab = 0.5 * mpar * valabsq
+ Eflab = 0.5 * Af*amu * vflabsq
+
+
+*=======================================================================*
+* Calculate weight of event
+*=======================================================================*
+
+* Angle between beta and neutrino
+ dot_bn = pbeta(1)*pnu(1) + pbeta(2)*pnu(2) + pbeta(3)*pnu(3)
+ cos_bn = dot_bn / (rpbeta*rpnu)
+
+* Angle between beta and alpha
+ dot_ba = pbeta(1)*pacm(1)+pbeta(2)*pacm(2)+pbeta(3)*pacm(3)
+ cos_ba = dot_ba / (rpbeta*rpacm)
+ angle_ba = acos(cos_ba)*180./pi
+
+* Angle between neutrino and alpha
+ dot_na = pnu(1)*pacm(1) + pnu(2)*pacm(2) + pnu(3)*pacm(3)
+ cos_na = dot_na / (rpnu*rpacm)
+
+* Triple correlation amplitude
+ va = rpacm/mpar
+ E0 = T0 + me
+ A = triple(Ai,rpbeta,cos_bn,cos_ba,cos_na,va,E0)
+
+* Weight of event
+ W = A
+ if (W.gt.Wmax) then
+ Wmax = W
+ endif
+
+
+*=======================================================================*
+* Calculate lab angle between alpha and final nucleus
+*=======================================================================*
+
+ dotprod = valab(1)*vflab(1) + valab(2)*vflab(2)
+ + + valab(3)*vflab(3)
+ cos_af = dotprod/sqrt(valabsq)/sqrt(vflabsq)
+ angle_af = acos(cos_af)*180./pi
+ if (angle_af.lt.angle_min) then
+ angle_min = angle_af
+ endif
+
+
+*=======================================================================*
+* Checks of self consistency
+*=======================================================================*
+
+* Check energy conservation
+ Etot = Ealab + Eflab + T_beta + E_nu
+ if (check) then
+ write(6,*) 'E_f-E_i and recoil energy (keV):'
+ + ,Etot-(Qbeta-alpha_thres),Erec
+ endif
+
+* Check momentum conservation
+ ptot1 = mpar*valab(1) + Af*amu*vflab(1) + pbeta(1) + pnu(1)
+ ptot2 = mpar*valab(2) + Af*amu*vflab(2) + pbeta(2) + pnu(2)
+ ptot3 = mpar*valab(3) + Af*amu*vflab(3) + pbeta(3) + pnu(3)
+ if (check) then
+ write(6,*) 'Ptot (MeV/c):'
+ + ,sqrt(ptot1**2.+ptot2**2.+ptot3**2.)
+ read(*,*)
+ endif
+
+
+*=======================================================================*
+* Plots and write data to files
+*=======================================================================*
+
+ call hfill(nhis_rec,real(Erec),0.,real(W))
+ call hfill(nhis_eb,real(T_beta),0.,real(W))
+ call hfill(nhis_ea,real(Ealab),0.,real(W))
+ call hfill(nhis_ef,real(Eflab),0.,real(W))
+ call hfill(nhis_ba,real(angle_ba),0.,real(W))
+ call hfill(nhis_af,real(angle_af),0.,real(W))
+ Wsum = Wsum + W
+ Esum = Esum + W*Erec
+ write(41,*) ieve, T_beta, Ealab, Eflab, angle_ba, angle_af, W
+
+
+ 99 ieve = ieve+1
+ enddo ! loop over ieve
+ write(41,*) 'Maximum weight:', Wmax
+ close(41)
+
+
+ write(6,*) 'Average recoil energy (keV):', Esum/Wsum
+ write(6,*) 'Minimum lab angle between alpha '
+ + //'and final nuclues (degree):', angle_min
+
+ call hrput(0,histout,'n')
+
+ END
+
+
+
+
+* Beta-decay phase space
+* T0me = total lepton kinetic energy in units of 511 keV
+* t = beta kinetic energy in units of 511 keV
+ real function phasespace(t)
+ COMMON /PHSP/T0me
+ q = T0me
+ phasespace = (t**2.+2.*t)**0.5*(q-t)**2.*(t+1.)
+ end
+
+
+* Triple-correlation function
+ real function triple(Ai,pb,cosbn,cosba,cosna,va,E0)
+ real pb,cosba,cosna
+ real me,Eb,M,E0,va
+ PARAMETER(me = 510.99906)
+ parameter(amu = 931494.0)
+ real r,v,bn,ab,an,xi,Theta,g2g1,g12g1,c,c0
+ common /COEFF/xi,Theta,a3
+
+ M = Ai*amu
+ Eb = sqrt( pb**2. + me**2.)
+
+ r = pb/Eb
+ v = va
+ bn = cosbn
+ ab = cosba
+ an = cosna
+
+ c = 2.*Eb/(M*v)
+ c0 = 2.*E0/(M*v)
+ c = max(0.,c)
+ c0 = max(0.,c0)
+
+ g2g1 = 1./3.*(2.*a3+1.)
+ g12g1 = 1./2.*Theta*(a3-1.)
+
+ d1 = -c
+ d2 = -c*g2g1
+ d3 = -(c0-c)
+ d4 = -(c0-c)*g2g1
+ d8 = 1./10.*xi*2.*c*(-g12g1)
+ d9 = 1./10.*xi*2.*(c0-c)*(-g12g1)
+
+ triple = 1. + g2g1*r*bn + g12g1*1./10.*xi*r*(ab*an-1./3.*bn)
+ + + (d1-1./5.*d9)*r*ab + (d2-2./5.*d8)*r**2.*ab*bn
+ + + (d3-1./5.*r**2.*d8)*an + (d4-2./5.*d9)*r*an*bn
+ + + d8*r**2.*ab**2.*an + d9*r*ab*an**2.
+
+ end
+
+
+
+* Dummy routines
+
+ subroutine dummyini
+ implicit none
+ real X,Y
+ real A(100)
+ integer i
+ real RNDM
+ external RNDM
+ do i=1,100
+ y = RNDM(X)
+ a(i) = y
+ enddo
+ call FLPSOR(A,20)
+ end
+
+ FUNCTION ran2(idummy)
+ IMPLICIT real*8 (A-H,O-Z)
+ parameter (IM1=2147483563,IM2=2147483399,AM=1d0/IM1,
+ *IMM1=IM1-1,
+ *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
+ *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2d-7,RNMX=1d0-EPS)
+ integer idum2,j,k,iv(NTAB),iy
+ SAVE iv,iy,idum2
+ DATA idum2/123456789/, iv/NTAB*0/, iy/0/
+ idum=idummy
+ if (idum.le.0) then
+ idum=max(-idum,1)
+ idum2=idum
+ do 11 j=NTAB+8,1,-1
+ k=idum/IQ1
+ idum=IA1*(idum-k*IQ1)-k*IR1
+ if (idum.lt.0) idum=idum+IM1
+ if (j.le.NTAB) iv(j)=idum
+11 continue
+ iy=iv(1)
+ endif
+ k=idum/IQ1
+ idum=IA1*(idum-k*IQ1)-k*IR1
+ if (idum.lt.0) idum=idum+IM1
+ k=idum2/IQ2
+ idum2=IA2*(idum2-k*IQ2)-k*IR2
+ if (idum2.lt.0) idum2=idum2+IM2
+ j=1+iy/NDIV
+ iy=iv(j)-idum2
+ iv(j)=idum
+ if(iy.lt.1)iy=iy+IMM1
+ ran2=min(AM*iy,RNMX)
+ return
+ end
+
diff --git a/description.pdf b/description.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..44def6aca6b5efe44fcfbd0ebd2d6969888e96ae
Binary files /dev/null and b/description.pdf differ
diff --git a/input0 b/input0
new file mode 100644
index 0000000000000000000000000000000000000000..9489b5be1858c655e8daa73256964f8466cc76b6
--- /dev/null
+++ b/input0
@@ -0,0 +1,9 @@
+ 100000
+ 13370.
+ 7366.6
+ 7654.
+ 12.
+ 2
+ 1.
+ 0.
+ 0.
diff --git a/spincoeff.f b/spincoeff.f
new file mode 100644
index 0000000000000000000000000000000000000000..05132f3b3bb0496840e5a74f751631deb80fe18b
--- /dev/null
+++ b/spincoeff.f
@@ -0,0 +1,56 @@
+ program SPINCOEFF
+
+* Oliver S. Kirsebom
+* 20 Jan 2011
+
+* This program calculates the spin-dependent
+* coefficients, Xi and Theta, which enter the
+* formula for the beta-nu-alpha triple
+* correlation.
+*
+* The subroutine tRACAW from CERNLIB is used
+* to calculate the Racah W-coefficient.
+*
+
+
+
+
+ REAL J,Ji,Jii,L
+
+* J = spin of parent nucleus
+* Ji = spin of daughter nucleus
+* Jii = spin of final nucleus
+* L = relative angular momentum in final breakup
+
+
+ write(6,*) 'Spin of parent nucleus? (J)'
+ read(*,*) J
+ write(6,*) 'Spin of daughter nucleus? (J`)'
+ read(*,*) Ji
+ write(6,*) 'Spin of final nucleus? (J``)'
+ read(*,*) Jii
+ write(6,*) 'Orbital angular momentum? (L)'
+ read(*,*) L
+
+ if (L.gt.0.) then
+ a = sqrt(L*(L+1.)*(2.*L+1.)/(2.*L-1.)/(2.*L+3.))
+ b = sqrt((2.*Ji-1.)*(2.*Ji+1.)*(2.*Ji+3.)/Ji/(Ji+1.))
+ c = rRACAW(2,Ji,L,Jii,Ji,L)
+ xi = 10.*a*b*c
+ else
+ xi = 0.
+ endif
+ write(6,*) 'Xi=',xi
+
+ if (Ji.gt.0.) then
+ a = (-1.)**(Ji-J)
+ b = sqrt(30.*Ji*(Ji+1.)*(2.*Ji+1.)/(2.*Ji-1.)/(2.*Ji+3.))
+ c = rRACAW(Ji,1.,Ji,1.,J,2.)
+ Theta = a*b*c
+ else
+ Theta = 0.
+ endif
+ write(6,*) 'Theta=',Theta
+
+
+ END
diff --git a/spincoeff.pdf b/spincoeff.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..822a1b97945a49ac9af8ba7b71b47dfdc3c04950
Binary files /dev/null and b/spincoeff.pdf differ