Adding ratios of mass constants to physical constants.

Also changed the version number.

src/groupr.f90

@@ -5581,12 +5581,6 @@ subroutine parts(mfd,mtd)
    use endf ! provides iverf
    ! externals
    integer::mfd,mtd
-   ! internals
-   real(kr),parameter::awr1=amassp/amassn
-   real(kr),parameter::awr2=amassd/amassn
-   real(kr),parameter::awr3=amasst/amassn
-   real(kr),parameter::awr4=amassh/amassn
-   real(kr),parameter::awr5=amassa/amassn
 
    law=0
    zap=1
@@ -5602,35 +5596,35 @@ subroutine parts(mfd,mtd)
          endif
       else if (mtd.ge.700.and.mtd.lt.720) then
          zap=1001
-         aprime=awr1
+         aprime=pnratio
          if (mfd.ne.31) then
             aprime=awr+1-aprime
             law=4
          endif
       else if (mtd.ge.720.and.mtd.lt.740) then
          zap=1002
-         aprime=awr2
+         aprime=dnratio
          if (mfd.ne.32) then
             aprime=awr+1-aprime
             law=4
          endif
       else if (mtd.ge.740.and.mtd.lt.760) then
          zap=1003
-         aprime=awr3
+         aprime=tnratio
          if (mfd.ne.33) then
             aprime=awr+1-aprime
             law=4
          endif
       else if (mtd.ge.760.and.mtd.lt.780) then
          zap=2003
-         aprime=awr4
+         aprime=hnratio
          if (mfd.ne.34) then
             aprime=awr+1-aprime
             law=4
          endif
       else if (mtd.ge.780.and.mtd.lt.800) then
          zap=2004
-         aprime=awr5
+         aprime=anratio
          if (mfd.ne.35) then
             aprime=awr+1-aprime
             law=4
@@ -5654,35 +5648,35 @@ subroutine parts(mfd,mtd)
          endif
       else if (mtd.ge.600.and.mtd.lt.650) then
          zap=1001
-         aprime=awr1
+         aprime=pnratio
          if (mfd.ne.31) then
             aprime=awr+1-aprime
             law=4
          endif
       else if (mtd.ge.650.and.mtd.lt.700) then
          zap=1002
-         aprime=awr2
+         aprime=dnratio
          if (mfd.ne.32) then
             aprime=awr+1-aprime
             law=4
          endif
       else if (mtd.ge.700.and.mtd.lt.750) then
          zap=1003
-         aprime=awr3
+         aprime=tnratio
          if (mfd.ne.33) then
             aprime=awr+1-aprime
             law=4
          endif
       else if (mtd.ge.750.and.mtd.lt.800) then
          zap=2003
-         aprime=awr4
+         aprime=hnratio
          if (mfd.ne.34) then
             aprime=awr+1-aprime
             law=4
          endif
       else if(mtd.ge.800.and.mtd.lt.850) then
          zap=2004
-         aprime=awr5
+         aprime=anratio
          if (mfd.ne.35) then
             aprime=awr+1-aprime
             law=4

src/heatr.f90

@@ -1878,11 +1878,6 @@ subroutine disbar(ee,ebar,dame,q,za,awr,ntape,matd,mtd,a,nwa)
       8.84675983E-03_kr, 6.50445797E-03_kr, 4.14703326E-03_kr,&
       1.78328072E-03_kr/)
    real(kr),parameter::edis=25.e0_kr
-   real(kr),parameter::am1=amassp/amassn
-   real(kr),parameter::am2=amassd/amassn
-   real(kr),parameter::am3=amasst/amassn
-   real(kr),parameter::am4=amassh/amassn
-   real(kr),parameter::am5=amassa/amassn
    real(kr),parameter::step=1.1e0_kr
    real(kr),parameter::small=1.e-10_kr
    real(kr),parameter::zero=0
@@ -1908,16 +1903,16 @@ subroutine disbar(ee,ebar,dame,q,za,awr,ntape,matd,mtd,a,nwa)
       endif
       thresh=((awr+1)/awr)*(-q)
       awp=1
-      if (iverf.le.5.and.mtd.ge.700.and.mtd.lt.718) awp=am1
-      if (iverf.le.5.and.mtd.ge.720.and.mtd.lt.738) awp=am2
-      if (iverf.le.5.and.mtd.ge.740.and.mtd.lt.758) awp=am3
-      if (iverf.le.5.and.mtd.ge.760.and.mtd.lt.778) awp=am4
-      if (iverf.le.5.and.mtd.ge.780.and.mtd.lt.798) awp=am5
-      if (iverf.ge.6.and.mtd.ge.600.and.mtd.lt.649) awp=am1
-      if (iverf.ge.6.and.mtd.ge.650.and.mtd.lt.699) awp=am2
-      if (iverf.ge.6.and.mtd.ge.700.and.mtd.lt.749) awp=am3
-      if (iverf.ge.6.and.mtd.ge.750.and.mtd.lt.799) awp=am4
-      if (iverf.ge.6.and.mtd.ge.800.and.mtd.lt.849) awp=am5
+      if (iverf.le.5.and.mtd.ge.700.and.mtd.lt.718) awp=pnratio
+      if (iverf.le.5.and.mtd.ge.720.and.mtd.lt.738) awp=dnratio
+      if (iverf.le.5.and.mtd.ge.740.and.mtd.lt.758) awp=tnratio
+      if (iverf.le.5.and.mtd.ge.760.and.mtd.lt.778) awp=hnratio
+      if (iverf.le.5.and.mtd.ge.780.and.mtd.lt.798) awp=anratio
+      if (iverf.ge.6.and.mtd.ge.600.and.mtd.lt.649) awp=pnratio
+      if (iverf.ge.6.and.mtd.ge.650.and.mtd.lt.699) awp=dnratio
+      if (iverf.ge.6.and.mtd.ge.700.and.mtd.lt.749) awp=tnratio
+      if (iverf.ge.6.and.mtd.ge.750.and.mtd.lt.799) awp=hnratio
+      if (iverf.ge.6.and.mtd.ge.800.and.mtd.lt.849) awp=anratio
       afact=awp/((awr+1)**2)
       arat=awp/(awr+1-awp)
       en=0

src/phys.f90

@@ -40,6 +40,12 @@ module physics
    real(kr),parameter,public::amassh=3.014932234673e0_kr !hellion (3)
    real(kr),parameter,public::amassa=4.001506179127e0_kr  !alpha
    real(kr),parameter,public::amasse=5.485799090e-4_kr  !electron
+   real(kr),parameter,public::pnratio=amassp/amassn ! proton/neutron mass
+   real(kr),parameter,public::dnratio=amassd/amassn ! deuteron/neutron mass
+   real(kr),parameter,public::tnratio=amasst/amassn ! triton/neutron mass
+   real(kr),parameter,public::hnratio=amassh/amassn ! helion/neutron mass
+   real(kr),parameter,public::anratio=amassa/amassn ! alpha/neutron mass
+
    real(kr),parameter,public::epair=amasse*amu*clight*clight/ev
 end module physics
 

src/samm.f90

@@ -1295,27 +1295,27 @@ subroutine ppdefs(ier)
          kza(i,ier)=0
          lpent(i,ier)=1
       else if (mt(i,ier).eq.103) then  ! proton
-         ema(i,ier)=amassp/amassn
+         ema(i,ier)=pnratio
          spina(i,ier)=0.5e0_kr
          kza(i,ier)=1
          lpent(i,ier)=1
       else if (mt(i,ier).eq.104) then  ! deuteron
-         ema(i,ier)=amassd/amassn
+         ema(i,ier)=dnratio
          spina(i,ier)=1
          kza(i,ier)=1
          lpent(i,ier)=1
       else if (mt(i,ier).eq.105) then  ! triton
-         ema(i,ier)=amasst/amassn
+         ema(i,ier)=tnratio
          spina(i,ier)=0.5e0_kr
          kza(i,ier)=1
          lpent(i,ier)=1
       else if (mt(i,ier).eq.106) then  ! he3
-         ema(i,ier)=amassh/amassn
+         ema(i,ier)=hnratio
          spina(i,ier)=0.5e0_kr
          kza(i,ier)=2
          lpent(i,ier)=1
       else if (mt(i,ier).eq.107) then  ! alpha
-         ema(i,ier)=amassa/amassn
+         ema(i,ier)=anratio
          spina(i,ier)=0
          pa(i,ier)=0
          kza(i,ier)=2
@@ -1719,16 +1719,15 @@ subroutine fxradi(ier)
    integer::ippx,kgroup,ichan
    real(kr)::ff,twomhb,etac
    real(kr)::factor,alabcm,aa,redmas,z
-   real(kr),parameter::emneut=amassn               !neutron mass in amu
    real(kr),parameter::hbarrr=hbar/ev              !hbar in eV-s
    real(kr),parameter::amuevv=amu*clight*clight/ev !amu in eV
    real(kr),parameter::cspeed=clight/100           !c in m/s (not cm/s!)
    real(kr),parameter::fininv=1.e0_kr/finstri      !fine structure constant
    real(kr),parameter::zero=0
 
    ff=1.e+15_kr
-   twomhb=sqrt(2*emneut*amuevv)/(hbarrr*ff*cspeed)
-   etac=fininv*amuevv/(hbarrr*ff*cspeed)*emneut
+   twomhb=sqrt(2*amassn*amuevv)/(hbarrr*ff*cspeed)
+   etac=fininv*amuevv/(hbarrr*ff*cspeed)*amassn
 
    do ippx=1,nppm(ier)
       if (ema(ippx,ier).eq.zero) then

src/thermr.f90

@@ -2300,7 +2300,7 @@ subroutine calcem(temp,itemp,iold,inew,ne,nex)
       write(nsyso,'(/''            mu            theta      dsigma/dmu'')')
       do i=1,nmu
          write(nsyso,'(i5,1x,f15.8,1x,f12.4,1x,1p,e14.7)') i,uj(i),&
-           acos(uj(i))*180/pi,sj(i)/2
+           acos(uj(i))*180.0/pi,sj(i)/2.0
       enddo
    endif
 

src/vers.f90

@@ -3,7 +3,7 @@ module version
    ! These values are updated during the NJOY revision-control process.
    implicit none
    private
-   character(8),public::vers='2016.40'
-   character(8),public::vday='20Aug18'
+   character(8),public::vers='2016.41'
+   character(8),public::vday='12Sep18'
 end module version