function nucrx_tdep_speed(tk,sc,irx) % These are Bosch-Hale cross sections format long e % tk is T in keV; n is the order of the quadrature, sc = scale fac avog=6.02214129e23; boltz1=0.086173324; boltz2=1.3806488e-10; conv=1.6021e-33; % T(d,n)He3 if irx==1 m1=3.0168492; m2=2.0141078; bb=34.3827; a=[6.927e4 7.454e8 2.050e6 5.2002e4 0]; b=[6.38e1 -9.95e-1 6.981e-5 1.728e-4]; else end %He3(d,p)He4 if irx==2 m1=3.0160293; m2=2.0141078; bb=68.7508; a=[5.7501e6 2.5226e3 4.5566e1 0 0]; b=[-3.1995e-3 -8.553e-6 5.9014e-8 0]; else end %D(d,p)T if irx==3 m1=2.0141078; m2=2.0141078; bb=31.3970; a=[5.5576e4 2.1054e2 -3.2638e-2 1.4987e-6 1.8181e-10]; b=[0 0 0 0]; else end %D(d,n)He3 if irx==4 m1=2.0141078; m2=2.0141078; bb=31.3970; a=[5.3701e4 3.3027e2 -1.2706e-1 2.9327e-5 -2.5151e-9]; b=[0 0 0 0]; else end mu=(m1*m2/(m1+m2))/avog; t6=tk/boltz1; s=0.; for n=4:2:30 ptwt=abspeed2(n); x=ptwt(:,1); w=ptwt(:,2); bhat=(bb/sqrt(tk))*ones(1,n); bigw=w./(x.*exp(-x.*x)); ws=sc*bigw; xs=sc*x; e=xs.^2*tk; sigma1=a(1)+e.*(a(2)+e.*(a(3)+e.*(a(4)+e.*a(5)))); sigma2=1+e.*(b(1)+e.*(b(2)+e.*(b(3)+e*b(4)))); sigma=sigma1./sigma2; sigma=1e-3*sigma; fx=xs.*exp(-xs.*xs-bhat'./xs).*sigma; xn=2*sum(ws.*fx); rate=2*conv*sqrt(2/(pi*mu*boltz2*t6))*xn; fprintf('%i %16.4e\n',n,rate) end end function ptwt = abspeed2(m) load abp1.dat; n=80; a=abp1(1:m,1); b=abp1(2:m,2); rtb=sqrt(b); t=diag(rtb,-1)+diag(a)+diag(rtb,1); [f,lambda]=eig(t); pt=diag(lambda); for i=1:m wt(i)=f(1,i)^2/2; end ptwt=[pt,wt']; end