c **********************************************************************
c Argonne v18pq potential package 
c
c prepared 4 Aug 2003 by R. B. Wiringa, Physics Division, 
c Argonne National Laboratory, Argonne, IL 60439
c e-mail: wiringa@anl.gov 
c 
c references: 
c  "Quadratic momentum dependence in the nucleon-nucleon interaction"
c   R. B. Wiringa, A. Arriaga, and V. R. Pandharipande
c   submitted to Physical Review C on 4 Jun 2003
c   http://arxiv.org/abs/nucl-th/0306018
c   
c  "Accurate nucleon-nucleon potential with charge-independence breaking"
c   R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla, 
c   Physical Review C51, 38 (1995) - WSS95.
c
c this file contains 4 subroutines:
c   subroutine av18pqpw(lpot,l,s,j,t,t1z,t2z,r,v,w,wp,wdp)
c   subroutine av18pqop(lpot,r,vnn,vpp)
c   subroutine empot(lpot,r,vem)
c   subroutine consts(hc,mpi0,mpic,mp,mn,alpha,mup,mun)
c
c   av18pqpw gives the full potential in a particular partial wave and
c            charge state including separate local and p-dependent parts
c            and 1st and 2nd derivatives of the p-dependent parts
c   av18pqop gives the strong interaction part in operator format
c            including separate local and p-dependent parts
c            and 1st and 2nd derivatives of the p-dependent parts
c   empot  gives the electromagnetic part in operator format
c   consts gives values of fundamental constants and masses used
c
c notes: 
c 1) Argonne v18pq is a variant of the Argonne v18 potential in which
c    L**2 operators have been generalized to p**2 operators, and the
c    quadratic spin-orbit operator (L.S)**2 has been changed to the
c    alternate form Q12 = 3L.S - (3+s1.s2)L**2 + 6(L.S)**2 ;
c    analytic 1st and 2nd derivatives of the radial functions 
c    associated with the p**2 operators are provided
c
c 2) empot does not include the energy-dependence of the Coulomb
c    interaction used in eq.(4) of WSS95, i.e., it uses alpha, 
c    not alpha'.
c
c 3) the vacuum polarization in empot is a short-range approximation 
c    to eq.(7) suitable for bound states, but not for scattering.
c    it is taken from eq.(3.13) of Rev. Mod. Phys. 44, 48 (1972)
c
c    8/28/97 error in this formula detected and corrected:
c    should be -(gamma+5/6) instead of printed (-gamma+5/6)
c
c 4) these subroutines should be compiled with a compiler option
c    that forces all floating point constants to be evaluated at
c    real*8 significance, e.g., on an IBM RS6000 the xlf compiler
c    option qdpc=e should be used; on SGI machines, the -r8 option
c    should be used; on a Cray no action is needed.
c    if such an option is not available and the default precision is 
c    real*4 (32 bits), then all constants should be explicitly 
c    converted to double precision by appending a D0.
c
c *id* av18pqpw ********************************************************
c subroutine for partial-wave projection of argonne v18pq potential
c calls subroutines av18pqop, empot
c ----------------------------------------------------------------------
c arguments for av18pqpw
c lpot: switch for potential choice
c       = 1 : v18pq
c l:    orbital angular momentum of pair (0,1,2,...)
c s:    total spin of pair (0 or 1)
c j:    total angular momentum of pair (0,1,2,...)
c t:    total isospin of pair (0 or 1)
c t1z:  isospin of particle 1 (1 for p, -1 for n)
c t2z:     "    "     "     2 (1 for p, -1 for n)
c r:    separation in fm
c v:    local potential in MeV (2x2 array)
c       (all local strong and em terms, including L**2 part of p**2)
c w:    p**2 potential in MeV-fm**2 (2x2 array)
c wp:   1st derivative of w in MeV-fm
c wdp:  2nd derivative of w in MeV
c ----------------------------------------------------------------------
c order of terms in v(l,m), w(l,m), etc.:
c      single channel                 coupled channel (l=j-1,s=1)
c      v(1,1) = v(l,s,j,t,t1z,t2z)    v(1,1) = v(l,s,j,t,t1z,t2z)
c      v(2,1) = 0                     v(2,1) = v(l<->l+2)
c      v(1,2) = 0                     v(1,2) = v(l<->l+2)
c      v(2,2) = 0                     v(2,2) = v(l+2,s,j,t,t1z,t2z)
c ----------------------------------------------------------------------
      subroutine av18pqpw(lpot,l,s,j,t,t1z,t2z,r,v,w,wp,wdp)
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      integer l,s,j,t,t1z,t2z,s1ds2,t1dt2,t12
      dimension vnn(18),vpp(12),vem(14),v(2,2),w(2,2),wp(2,2),wdp(2,2)
c ------------------------
c strong interaction terms
c ------------------------
      call av18pqop(lpot,r,vnn,vpp)
      s1ds2=4*s-3
      t1dt2=4*t-3
      t12=3*t1z*t2z-t1dt2
      vc=vnn(1)+t1dt2*vnn(2)+s1ds2*vnn(3)+s1ds2*t1dt2*vnn(4)
     &  +t12*vnn(15)+s1ds2*t12*vnn(16)+(t1z+t2z)*vnn(18)
      vt=vnn(5)+t1dt2*vnn(6)+t12*vnn(17)
      vls=vnn(7)+t1dt2*vnn(8)
      vl2=vnn(9)+t1dt2*vnn(10)+s1ds2*vnn(11)+s1ds2*t1dt2*vnn(12)
      vq=vnn(13)+t1dt2*vnn(14)
      vp2=vpp(1)+t1dt2*vpp(2)+s1ds2*vpp(3)+s1ds2*t1dt2*vpp(4)
      vp2p=vpp(5)+t1dt2*vpp(6)+s1ds2*vpp(7)+s1ds2*t1dt2*vpp(8)
      vp2dp=vpp(9)+t1dt2*vpp(10)+s1ds2*vpp(11)+s1ds2*t1dt2*vpp(12)
c ---------------------
c electromagnetic terms
c ---------------------
      call empot(lpot,r,vem)
      if (t1z+t2z) 10,20,30
   10 vc=vc+s1ds2*vem(7)
      vt=vt+vem(10)
      go to 40
   20 vc=vc+vem(5)+s1ds2*vem(8)
      vt=vt+vem(11)
      vls=vls+vem(14)
      go to 40
   30 vc=vc+vem(1)+vem(2)+vem(3)+vem(4)+s1ds2*vem(6)
      vt=vt+vem(9)
      vls=vls+vem(12)
   40 continue
c ---------------------
c partial wave assembly
c ---------------------
   50 ncc=1
      if (s.eq.1.and.j.gt.l) ncc=2
      if (ncc.eq.1) then
        s12=0.
        if (s.eq.1.and.l.eq.j) s12=2.
        if (l.eq.(j+1)) s12=-2.*(j+2.)/(2.*j+1.)
        ls=(j*(j+1)-l*(l+1)-s*(s+1))/2
        lq=3*ls-4*s*l*(l+1)+6*ls**2
        v(1,1)=vc+s12*vt+ls*vls+l*(l+1)*vl2+lq*vq
        v(2,1)=0
        v(1,2)=0
        v(2,2)=0
      else if (ncc.eq.2) then
        s12m=-2.*(j-1.)/(2.*j+1.)
        s12=sqrt(36.*j*(j+1))/(2.*j+1.)
        s12p=-2.*(j+2.)/(2.*j+1.)
        lsm=j-1
        lsp=-(j+2)
        lqm=3*lsm-4*s*l*(l+1)+6*lsm**2
        lqp=3*lsp-4*s*(l+2)*(l+3)+6*lsp**2
        v(1,1)=vc+s12m*vt+lsm*vls+l*(l+1)*vl2+lqm*vq
        v(2,1)=s12*vt
        v(1,2)=s12*vt
        v(2,2)=vc+s12p*vt+lsp*vls+(l+2)*(l+3)*vl2+lqp*vq
      end if
c ---------------------
c p**2 & derivatives
c ---------------------
      w(1,1)=vp2
      w(2,1)=0
      w(1,2)=0
      w(2,2)=(ncc-1)*vp2
      wp(1,1)=vp2p
      wp(2,1)=0
      wp(1,2)=0
      wp(2,2)=(ncc-1)*vp2p
      wdp(1,1)=vp2dp
      wdp(2,1)=0
      wdp(1,2)=0
      wdp(2,2)=(ncc-1)*vp2dp
      return
      end
c *id* av18pqop ********************************************************
c subroutine for strong interaction part of argonne v18pq potential 
c calls subroutine consts
c ----------------------------------------------------------------------
c this subroutine returns two sets of radial functions: vnn and vpp
c
c vnn are the radial functions corresponding to the static and 
c L-dependent operators, including vnn(9:12) which accompany the 
c L**2 operator embedded in p**2; these are essentially the same
c operators as used in av18, except (L.S)**2 is replaced by Q12
c
c vpp are the radial functions that accompany the p-dependent
c operators, including vpp(1:4) = 0.5*r**2*vnn(9:12);
c the vpp(5:8) are the first derivatives of these radial functions
c and vpp(9:12) are the second derivatives
c ----------------------------------------------------------------------
c arguments for av18pot
c lpot: switch for potential choice
c       = 1 : v18pq
c r:    separation in fm
c vnn:  output potential in MeV (18 component array)
c vpp:  p**2 output potential in MeV-fm**2 and its 1st derivative in
c       MeV-fm and 2nd derivative in MeV (12 component array)
c ----------------------------------------------------------------------
c order of operators l in vnn(l):
c l:    1=1                              2=t1.t2
c       3=s1.s2                          4=(s1.s2)(t1.t2)
c       5=S12 [=3(s1.r)(s2.r)-s1.s2]     6=S12(t1.t2)
c       7=L.S                            8=L.S(t1.t2)
c       9=L**2                          10=L**2(t1.t2)
c      11=L**2(s1.s2)                   12=L**2(s1.s2)(t1.t2)
c      13=Q12                           14=Q12(t1.t2)
c      15=T12 [=3*t1z*t2z-t1.t2]        16=(s1.s2)T12
c      17=S12*T12                       18=t1z+t2z
c where s1=sigma_1, t1=tau_1, t1z=tau_1(z), etc.
c ----------------------------------------------------------------------
c order of operators l in vpp(l):
c l:    1=p**2                           2=p**2(t1.t2)
c       3=p**2(s1.s2)                    4=p**2(s1.s2)(t1.t2)
c       5=[p**2]'                        6=[p**2(t1.t2)]'
c       7=[p**2(s1.s2)]'                 8=[p**2(s1.s2)(t1.t2)]'
c       9=[p**2]"                       10=[p**2(t1.t2)]"
c      11=[p**2(s1.s2)]"                12=[p**2(s1.s2)(t1.t2)]"
c where [ ]' denotes a 1st derivative and [ ]" a 2nd derivative
c ----------------------------------------------------------------------
      subroutine av18pqop(lpot,r,vnn,vpp)
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension vnn(18),vpp(12)
      real*8 mpi0,mpic,mp,mn,mup,mun
      real*8 mpi,mu0,muc,mu
      data small/1e-4/
      call consts(hc,mpi0,mpic,mp,mn,alpha,mup,mun)
      mpi=(mpi0+2.*mpic)/3.
      mu0=mpi0/hc
      muc=mpic/hc
      mu=mpi/hc
      fsq=.075
      cpi=2.1
      rws=.5
      aiws=5.
      x=mu*r
      x0=mu0*r
      xc=muc*r
      if (r.le.small) then
        tpi=3*cpi**2*r/mu**3
        ypi0=(mpi0/mpic)**2*(mpi0/3)*cpi*r/mu0
        ypic=(mpic/3)*cpi*r/muc
        tpi0=3*cpi*ypi0/mu0**2
        tpic=3*cpi*ypic/muc**2
      else
        expcut=exp(-cpi*r*r)
        rcut=1-expcut
        rcutp=2*cpi*r*expcut
        rcutdp=2*cpi*expcut*(1-2*cpi*r**2)
        ypi0=(mpi0/mpic)**2*(mpi0/3)*exp(-x0)*rcut/x0
        ypic=(mpic/3)*exp(-xc)*rcut/xc
        tpi0=(1+(3+3/x0)/x0)*ypi0*rcut
        tpic=(1+(3+3/xc)/xc)*ypic*rcut
        y=exp(-x)/x
        yp=-mu*(1+1/x)*y
        ydp=mu**2*(1+2*(1+1/x)/x)*y
        t=1+3*(1+1/x)/x
        tp=-mu*3*(1+2/x)/x**2
        tdp=mu**2*6*(1+3/x)/x**3
        tpi=t*y*rcut**2
        tpip=(tp*y+t*yp)*rcut**2+2*t*y*rcut*rcutp
        tpidp=(tdp*y+2*tp*yp+t*ydp)*rcut**2+4*(tp*y+t*yp)*rcut*rcutp
     &       +2*t*y*(rcut*rcutdp+rcutp**2)
      end if
      ypi0=fsq*ypi0
      ypic=fsq*ypic
      tpi0=fsq*tpi0
      tpic=fsq*tpic
      tpi2=tpi*tpi
      tpi2p=2*tpi*tpip
      tpi2dp=2*(tpi*tpidp+tpip**2)
      expws=exp((r-rws)*aiws)
      ws=1./(1.+expws)
      wsp=-aiws*expws*ws**2
      wsdp=aiws*wsp+2*wsp**2/ws
      ws0=1./(1+exp(-rws*aiws))
      wsp0=-aiws*exp(-rws*aiws)*ws0**2
      wsm=ws*(1.-wsp0*x/(mu*ws0))
      wsmp=wsp-wsp0*(x*wsp+mu*ws)/(mu*ws0)
      wsmdp=wsdp-wsp0*(x*wsdp+2*mu*wsp)/(mu*ws0)
      wsx=ws*x
      wsx2=wsx*x
      wsx2p=x**2*wsp+2*mu*x*ws
      wsx2dp=x**2*wsdp+4*mu*x*wsp+2*mu**2*ws
      dypi00=(mpi0/mpic)**2*(mpi0/3)*cpi/mu0
      dypic0=(mpic/3)*cpi/muc
      ypi0p=ypi0-fsq*dypi00*ws*r/ws0
      ypicp=ypic-fsq*dypic0*ws*r/ws0
c ----------------------------------------
      p11pp=  -9.882847*tpi2+2589.7742*wsm+2952.3910*wsx2+ypi0p
      p11np=  -9.882847*tpi2+2587.9836*wsm+2952.3910*wsx2-ypi0p+2*ypicp
      p11nn=  -9.882847*tpi2+2586.1930*wsm+2952.3910*wsx2+ypi0p
      pt1pp=   1.420069*tpi2 -453.5357*wsx -837.3820*wsx2+tpi0
      pt1np=   1.420069*tpi2 -453.5357*wsx -837.3820*wsx2-tpi0+2*tpic
      pt1nn=   1.420069*tpi2 -453.5357*wsx -837.3820*wsx2+tpi0
      pls1=   -1.749197*tpi2 -493.8470*wsm+1533.0637*wsx2
      pl211=  -0.008159*tpi2 +132.7694*wsm -169.8510*wsx2
      pls21=   0.135181*tpi2  -17.7975*wsm  -46.2542*wsx2
      p10=    -8.351808*tpi2+2325.5929*wsm +957.8091*wsx2-ypi0p-2*ypicp
      pt0=     1.327862*tpi2-1170.8528*wsx +580.5596*wsx2-tpi0-2*tpic
      pls0=    0.060223*tpi2  +58.3208*wsm -126.0235*wsx2
      pl210=  -0.023577*tpi2   +1.8164*wsm +127.2921*wsx2
      pls20=   0.000589*tpi2  -25.1123*wsm   -4.6897*wsx2
      p01pp= -10.518030*tpi2+2836.0715*wsm +651.1945*wsx2-3*ypi0p
      p01np= -10.812190*tpi2+2816.4190*wsm+1002.5300*wsx2
     &                                               -3*(-ypi0p+2*ypicp)
      p01nn= -10.518030*tpi2+2832.4903*wsm +651.1945*wsx2-3*ypi0p
      pl201=   0.134747*tpi2   -9.4691*wsm
      p00=    -4.739629*tpi2+1121.2225*wsm+2764.3395*wsx2
     &                                               -3*(-ypi0p-2*ypicp)
      pl200=  -0.227084*tpi2 +166.5629*wsm
c ----------------------------------------
      p11=(p11pp+p11nn+p11np)/3
      p11cd=(.5*(p11pp+p11nn)-p11np)/6
      p11cs=(p11pp-p11nn)/4
      pt1=(pt1pp+pt1nn+pt1np)/3
      pt1cd=(.5*(pt1pp+pt1nn)-pt1np)/6
      pt1cs=(pt1pp-pt1nn)/4
      p01=(p01pp+p01nn+p01np)/3
      p01cd=(.5*(p01pp+p01nn)-p01np)/6
      p01cs=(p01pp-p01nn)/4
c ----------------------------------------
      vnn(1)=.0625*(9*p11+3*p10+3*p01+p00)
      vnn(2)=.0625*(3*p11-3*p10  +p01-p00)
      vnn(3)=.0625*(3*p11  +p10-3*p01-p00)
      vnn(4)=.0625*(  p11  -p10  -p01+p00)
      vnn(5)=.25*(3*pt1+pt0)
      vnn(6)=.25*(  pt1-pt0)
      vnn(7)=.25*(3*pls1+pls0)
      vnn(8)=.25*(  pls1-pls0)
      vnn(9)= .0625*(9*pl211+3*pl210+3*pl201+pl200)
      vnn(10)=.0625*(3*pl211-3*pl210+  pl201-pl200)
      vnn(11)=.0625*(3*pl211+  pl210-3*pl201-pl200)
      vnn(12)=.0625*(  pl211-  pl210-  pl201+pl200)
      vnn(13)=.25*(3*pls21+pls20)
      vnn(14)=.25*(  pls21-pls20)
      vnn(15)=.25*(3*p11cd+p01cd)
      vnn(16)=.25*(  p11cd-p01cd)
      vnn(17)=pt1cd
      vnn(18)=p01cs
c -----------------------
c vpsq & derivative terms
c -----------------------
      pl2p11=   -0.008159*tpi2p  +132.7694*wsmp  -169.8510*wsx2p
      pl2p10=   -0.023577*tpi2p    +1.8164*wsmp  +127.2921*wsx2p
      pl2p01=    0.134747*tpi2p    -9.4691*wsmp
      pl2p00=   -0.227084*tpi2p  +166.5629*wsmp
      pl2dp11=  -0.008159*tpi2dp +132.7694*wsmdp -169.8510*wsx2dp
      pl2dp10=  -0.023577*tpi2dp   +1.8164*wsmdp +127.2921*wsx2dp
      pl2dp01=   0.134747*tpi2dp   -9.4691*wsmdp
      pl2dp00=  -0.227084*tpi2dp +166.5629*wsmdp
c ----------------------------------------
      psq11=.5*r**2*pl211
      psq10=.5*r**2*pl210
      psq01=.5*r**2*pl201
      psq00=.5*r**2*pl200
      psqp11=r*pl211+.5*r**2*pl2p11
      psqp10=r*pl210+.5*r**2*pl2p10
      psqp01=r*pl201+.5*r**2*pl2p01
      psqp00=r*pl200+.5*r**2*pl2p00
      psqdp11=pl211+2.*r*pl2p11+.5*r**2*pl2dp11 
      psqdp10=pl210+2.*r*pl2p10+.5*r**2*pl2dp10
      psqdp01=pl201+2.*r*pl2p01+.5*r**2*pl2dp01
      psqdp00=pl200+2.*r*pl2p00+.5*r**2*pl2dp00
      vpp(1)=.0625*(9*psq11+3*psq10+3*psq01+psq00)
      vpp(2)=.0625*(3*psq11-3*psq10  +psq01-psq00)
      vpp(3)=.0625*(3*psq11  +psq10-3*psq01-psq00)
      vpp(4)=.0625*(  psq11  -psq10  -psq01+psq00)
      vpp(5)=.0625*(9*psqp11+3*psqp10+3*psqp01+psqp00)
      vpp(6)=.0625*(3*psqp11-3*psqp10+  psqp01-psqp00)
      vpp(7)=.0625*(3*psqp11+  psqp10-3*psqp01-psqp00)
      vpp(8)=.0625*(  psqp11-  psqp10-  psqp01+psqp00)
      vpp(9)= .0625*(9*psqdp11+3*psqdp10+3*psqdp01+psqdp00)
      vpp(10)=.0625*(3*psqdp11-3*psqdp10+  psqdp01-psqdp00)
      vpp(11)=.0625*(3*psqdp11+  psqdp10-3*psqdp01-psqdp00)
      vpp(12)=.0625*(  psqdp11-  psqdp10-  psqdp01+psqdp00)
      return
      end
c *id* empot ***********************************************************
c subroutine for electromagnetic part of Argonne v18 potential
c for avn' models returns pp Coulomb only
c calls subroutine consts
c ----------------------------------------------------------------------
c arguments for empot
c lpot: switch for potential choice
c       = 1 : full EM potential
c       > 1 : C1(pp) only
c r:    input separation in fm
c vem:  output potential in MeV (14 component array)
c ----------------------------------------------------------------------
c order of operators in vem(l)
c l:    1=C1    (pp)          2=DF    (pp)          3=C2      (pp)
c       4=VP    (pp)                                5=C1      (np)
c       6=s1.s2 (pp)          7=s1.s2 (nn)          8=s1.s2   (np)
c       9=S12   (pp)         10=S12   (nn)         11=S12     (np)
c      12=L.S   (pp)         13=L.S   (nn)         14=L.S     (np)
c C1 = one-photon-exchange Coulomb with form factor
c C2 = two-photon-exchange Coulomb
c DF = Darwin-Foldy
c VP = vacuum polarization (short-range approximation)
c all other terms from magnetic moment (MM) interactions
c ----------------------------------------------------------------------
      subroutine empot(lpot,r,vem)
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension vem(14)
      real*8 mpi0,mpic,mp,mn,mup,mun
      real*8 kr,me,mr
      data small/1e-5/
      do 5 l=1,14
        vem(l)=0
    5 continue
      call consts(hc,mpi0,mpic,mp,mn,alpha,mup,mun)
      b=4.27
      br=b*r
      pi=acos(-1.)
      me=0.510999
      mr=mp*mn/(mp+mn)
      gamma=0.577216
      beta=.0189
      if (r.lt.small) then
        fcoulr=5*b/16
        ftr3=b**3*br**2/720
        flsr3=b**3/48
        kr=me*small/hc
      else
        fcoulr=(1-(1+11*br/16+3*br**2/16+br**3/48)*exp(-br))/r
        ftr3=(1-(1+br+br**2/2+br**3/6+br**4/24+br**5/144)*exp(-br))/r**3
        flsr3=(1-(1+br+br**2/2+7*br**3/48+br**4/48)*exp(-br))/r**3
        kr=me*r/hc
      end if
      fivp=-gamma-5./6.+abs(log(kr))+6*pi*kr/8
      fdelta=b**3*(1+br+br**2/3)*exp(-br)/16
      fnpr=b**3*(15+15*br+6*br**2+br**3)*exp(-br)/384
      vem(1)=alpha*hc*fcoulr
c ------------------------
      if (lpot.ge.2) return
c ------------------------
      vem(2)=-alpha*hc**3*fdelta/(4*mp**2)
      vem(3)=-vem(1)**2/mp
      vem(4)=2*alpha*vem(1)*fivp/(3*pi)
      vem(5)=alpha*hc*beta*fnpr
      vem(6)=-alpha*hc**3*mup**2*fdelta/(6*mp**2)
      vem(7)=-alpha*hc**3*mun**2*fdelta/(6*mn**2)
      vem(8)=-alpha*hc**3*mup*mun*fdelta/(6*mn*mp)
      vem(9)=-alpha*hc**3*mup**2*ftr3/(4*mp**2)
      vem(10)=-alpha*hc**3*mun**2*ftr3/(4*mn**2)
      vem(11)=-alpha*hc**3*mup*mun*ftr3/(4*mp*mn)
      vem(12)=-alpha*hc**3*(4*mup-1)*flsr3/(2*mp**2)
      vem(13)=0
      vem(14)=-alpha*hc**3*mun*flsr3/(2*mn*mr)
      return
      end
c *id* consts **********************************************************
c subroutine for constants in av18 potential
c ----------------------------------------------------------------------
c arguments for consts
c hc:    output value for hbar*c (MeV-fm)
c mpi0:    "      "    "  neutral pion mass (MeV)
c mpic:    "      "    "  charged pion mass (MeV)
c mp:      "      "    "  proton mass (MeV)
c mn:      "      "    "  neutron mass (MeV)
c alpha:   "      "    "  electromagnetic constant alpha
c mup:     "      "    "  proton magnetic moment (nm)
c mun:     "      "    "  neutron magnetic moment (nm)
c ----------------------------------------------------------------------
      subroutine consts(hc,mpi0,mpic,mp,mn,alpha,mup,mun)
      real*8 hc,mpi0,mpic,mp,mn,alpha,mup,mun
      hc=197.327053
      mpi0=134.9739
      mpic=139.5675
      mp=938.27231
      mn=939.56563
      alpha=1./137.035989
      mup= 2.7928474
      mun=-1.9130427
      return
      end
