c May 2006 c PRECOM2.for contains the most interconnected subprograms c ALPHA NUTRA c ASHELL OMEGA c COLLEC SINGLE c EMMY TCOR c FACT WELL c c******************************************************************** c function alpha(ppi,hpi,pnu,hnu,ee,ze,gpi,gnu) c A simplified version of OMEGA for use in TCOR for calculating c the correction function to the state densities with good isospin c c written in 1990 c c Calculates 2-component particle-hole state densities in the c simple equi-spacing model c (no pairing, shell or finite well depth corrections). c c called from: TCOR c c c * I am alpha and omega, the beginning and the end. c * I will give unto him that is athirst of the fountain of the c * water of life freely. (Rev. 21.6) c common /state/ fac(52), zee pour = 1.e+35 alpha = 0. ze = ee nppi = ppi nhpi = hpi npnu = pnu nhnu = hnu if (nppi.lt.0) go to 20 if (nhpi.lt.0) go to 20 if (npnu.lt.0) go to 20 if (nhnu.lt.0) go to 20 nppi1 = nppi + 1 nhpi1 = nhpi + 1 npnu1 = npnu + 1 nhnu1 = nhnu + 1 nph = nppi + nhpi + npnu + nhnu if (nph.lt.1) go to 20 rho = gpi**ppi / fac(nppi1) rho = rho * gpi**hpi / fac(nhpi1) rho = rho * gnu**pnu / fac(npnu1) rho = rho * gnu**hnu / fac(nhnu1) ph = nph c -------------------------------------- c 1. Calculate pauli correction function c -------------------------------------- qpi = amax1(ppi,hpi) epaul = ppi*(ppi+1.) + hpi*(hpi+1.) epaul = qpi*qpi - epaul*0.25 epaul = epaul/gpi thresh = qpi*qpi/gpi qnu = amax1(pnu,hnu) epauli = pnu*(pnu+1.) + hnu*(hnu+1.) epauli = qnu*qnu - epauli*0.25 epauli = epauli/gnu epauli = epauli + epaul thresh = qnu*qnu/gnu + thresh ze = ee - epauli edif = ee - thresh if (edif.gt.0.) then term3 = 0. act = 0. if (nppi.gt.0) act = 1. if (nhpi.gt.0) act = act+1. if (npnu.gt.0) act = act+1. if (nhnu.gt.0) act = act+1. act1 = 1./act act2 = act1 + act1 act3 = act2 + act1 if (qpi.gt.0.) then den = 5.7 + 0.6*act den = den * gpi * edif * act / ph den = 9.35 + den t3 = 0. if (nppi.gt.0) t3 = (ppi-act2)*(ppi-act3) + act1 if (nhpi.gt.0) t3 = (hpi-act2)*(hpi-act3) + act1 + t3 term3 = t3/(gpi*den) end if if (qnu.gt.0.) then den = 5.7 + 0.6*act den = den * gnu * edif * act / ph den = 9.35 + den t3 = 0 if (npnu.gt.0) t3 = (pnu-act2)*(pnu-act3) + act1 if (nhnu.gt.0) t3 = (hnu-act2)*(hnu-act3) + act1 + t3 term3 = t3/(gnu*den) + term3 end if ze = ze - term3 end if c -------------------------- c 2. Calculate state density c -------------------------- if (edif.le.0.) go to 20 if (ze.le.0.) go to 20 alpha = (ph-1.)*alog10(ze) - alog10(fac(nph)) alpha = rho * 10.**alpha if (alpha.lt.0.1) alpha=0. if (alpha.gt.pour) alpha=0. 20 return end c c********************************************************************* c function ashell(kp,api0,anu0,temp) c Shell corrected equilibrium level density parameters using simple c staircase weighting functions. c Always evaluated with the maximum weight in the center of the c shell gap (as for closed shell nuclei). c For other nuclei, take a weighted average with a0. c c Written in 1996 c c Called from EMMY c dimension asp(2), asp0(2) common /shell/ db(2), ds(7,2), epauli(7,14,2), dodd(7,2), fspan, 1 ghole(2), gpart(2), idif(7,2), ishell, mdnear(7,2), mdfar(7,2), 2 ndwn, nup2, range(2), rose(7), thresh(7,14,2) asp0(1) = api0 asp0(2) = anu0 span = fspan * temp do 12 ia=1,2 xmovf = range(ia) c ** Range is the number of mass units over which shell effects c ** wash out. idifa = iabs(idif(kp,ia)) difa = idifa if (difa.ge.xmovf) then asp(ia) = asp0(ia) go to 12 end if db2 = db(ia) * 0.5 wtmax = span / ds(kp,ia) ninter = wtmax + 1. c ** With this ninter, the last interval will have a fractional weight sum = 0. wtsum = 0. weight = wtmax + 1. do 10 j=1,ninter xj = j ehigh = xj * ds(kp,ia) elow = ehigh - ds(kp,ia) weight = weight - 1. if (weight.le.0.) go to 10 wtsum = wtsum + weight if (elow.gt.db2) then sum = sum + weight c ** Bins in this region have one s.p. state in them else if (elow.eq.db2) then sum = sum + (xj+1)*0.5*weight c ** Degen level is at low end of bin. Bin has 1 other s.p. st. else if (ehigh.eq.db2) then sum = sum + xj*0.5*weight c ** Degen level is at top of bin. Add in half of it. else if (ehigh.gt.db2) then sum = sum + xj*weight c ** This bin has the degenerate level inside it; maybe 1 other end if c ** The cases where ehigh 0.446, it is the energy of the phase transition. c ** In either case delta and c are set to zero for e(ia) < ephase. c ** This leads to a zero value for coll. if (xnrat.le.0.446) then ephase = ckp * (3.23*xnrat - 1.57*xnrat*xnrat) else ephase = ckp * (0.716 + 2.44*xnrat**2.17) end if if (e(ia).le.ephase) go to 20 erat = e(ia) / ckp delrat = 0.996 - 1.76*(xnrat**1.6)/(erat**0.68) c = ckp * delrat * delrat x = q / gres(kc,kp,ia) y = x*x + 4.*c/gres(kc,kp,ia) coll = coll - c + q*sqrt(y) - q*x 20 continue collec = coll return end c c******************************************************************** c subroutine emmy(kc,kp,multi) c c Written in 1982; revised 1990 c Constant temperature eqb state densities at low E added 1996 c c Particle emission rates from particle-hole states c and general two-component equilibrium evaporation rates. c Includes sum over allowed exit channel isopins in c an isopin conserved calculation. c Constant temperature part assumed spin cutoff factor is c constant below the match point. c c called from: PRECO c calls to: ASHELL, OMEGA, WELL c dimension prod(51), ref(3), ur(51) common /emit/ clos(13,13), fiss(3), fweis, fweist, fteq, 1 mxdime, mxdimp, mxdim2, ndwn1, nup1, preeq1(6,51), preeq2(2,51), 2 ptot, rho(13,13), ttot(6,51), weiss(6,51), weist(6,51) common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3), 1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6), 2 nsigin, nphd(3,7,2), rzz, sigin(7,51) common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin, 1 xncrit(7,2) common /shell/ db(2), ds(7,2), epauli(7,14,2), dodd(7,2), fspan, 1 ghole(2), gpart(2), idif(7,2), ishell, mdnear(7,2), mdfar(7,2), 2 ndwn, nup2, range(2), rose(7), thresh(7,14,2) common /state/ fac(52), ze common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3), 1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff common /welld/ ehole(2),ehbar,epbar,etest(2), fwd(2), vfull common /wmain/ w(3,13,14,51), strong(9,10,51), seq(6,51), 1 seq2(2,51), steq(6,51), steq2(2,51), xiso(9,10,51), xisoeq(2,51), 2 spre(2,51) flow = 1.e-18 spill = 1.e+18 umax = e - ben(kc,kp) if (umax.le.0.) go to 38 mkp = kp if (kp.gt.3) mkp = kp-3 jout = jpout(kp) + jnout(kp) iout = ipout(kc) + inout(kc) ap = jout apr = ap * (acom(kc)-ap) / acom(kc) sd = nsd(kp) rel = sd * apr * 3.72e-04 c -------------------------------- c 1. equilibrium evaporation rates c -------------------------------- eeh = e c ** EEH should be E-ROSE(7) but with 2ndary emiss, the 2nd CN can c ** have low E and need a constant temperature state density, c ** which extends down to zero excitation energy. c ** Since composite state densities cancel out in equilibrium c ** branching ratios, this makes no difference and the constant c ** temperature form is not used there. They are, however, used in c ** the residual nucleus state densities. c ** But be sure that same CN state density is used in calculating c ** the fission rates! api0 = 1.645 * gres(kc,kp,1) anu0 = 1.645 * gres(kc,kp,2) aspr0 = api0 + anu0 asp = 1.645 * (gres(kc,7,1)+gres(kc,7,2)) umatch = 2.5 + 150./(acom(kc)-ap) temp = sqrt(aspr0/umatch) - 1.25/umatch temp = 1. / temp c ** Matching energy determines transition to constant temp state c ** densities for residual nuclei in evaporation rates. if (ishell.gt.0) then aspr = ashell(kp,api0,anu0,temp) temp = sqrt(aspr/umatch) - 1.25/umatch temp = 1. / temp else aspr = aspr0 end if const = 2.* sqrt(aspr*umatch) x = (aspr*umatch)**0.25 * umatch const = exp(const) / x x = exp((umatch+rose(kp))/temp) const = const / x ref(1) = 2. * sqrt(asp*eeh) imax=0 if (iso.gt.0) then imax=2 if (ct(kc,kp,2).le.0.) imax=1 ee = eeh - esym(kc,7,1) ref(2) = 0. if (ee.gt.0.) ref(2) = 2. * sqrt(asp*ee) ref(3) = ref(2) end if c ** i=0 is for the T-mixed part of the evaporation rates. c ** i=1,2 are for the two resid. isospins in the T-conserved part do 10 ne=2,neps1(kc,kp) ee = eeh wrs = 0. wrt = 0. ur(ne) = umax - eps(ne) prod(ne) = eps(ne) * sigin(kp,ne) * rel do 6 i=0,imax u = ur(ne) if (i.eq.1) then ee = ee - esym(kc,7,1) if (ee.le.0.) go to 8 u = u - esym(kc,kp,1) else if (i.eq.2) then u = u + esym(kc,kp,1) - esym(kc,kp,2) end if if (u.le.0.) go to 8 u = u - rose(kp) if (u.ge.umatch) then quo = (u/ee)**1.25 * (aspr/asp)**0.25 prex = prod(ne) / quo ex = 2. * sqrt(aspr*u) else u = u + rose(kp) prex = prod(ne) * ee**1.25 * asp**0.25 prex = const * prex ex = u / temp end if ex = ex-ref(i+1) wr = prex * exp(ex) if (wr.lt.flow) wr=0. if (wr.gt.spill) wr=0. if (i.eq.0) then wrs = wr else x = wr*ct(kc,kp,i) wrt = wrt + x c ** If multiple preeq emission and T consv are both active, c ** then set fraction of primary emission going to lower T. if (multi.eq.0) go to 6 if (kc.gt.1) go to 6 if (kp.gt.2) go to 6 if (i.eq.1) then xisoeq(kp,ne) = 1. else if (wrt.gt.0.) xisoeq(kp,ne) = 1. - x/wrt end if end if 6 continue 8 if (iso.gt.0) weist(kp,ne) = wrt 10 weiss(kp,ne) = wrs c -------------------------------- c 2. preequilibrium emission rates c -------------------------------- if (kc.gt.1) then if (e.le.esym(kc,7,1)) go to 38 end if c ** This bypasses 2ndary preeq. emission with T conservation c ** when energetically forbidden. npilo1 = max0(jpin-ipout(kc),jpout(kp)) + 1 do 36 ne=2,neps1(kc,kp) u = ur(ne) i = 1 16 ee = ur(ne) - esym(kc,kp,i) if (ee.le.0.) go to 36 do 32 np1=ndwn+1,nup1 c ** np1 is the particle number index in the emitting nucleus. c ** np (below) is the particle number in the residual nucleus. c ** nh is the same in both nuclei, so nph is for the residual. mp1 = np1 if (kp.gt.3) mp1 = mxdimp-np1+1 if (iso.gt.0) then c ** Set indices for XISO, used in secondary particle emission. c ** Indices reflect configuration in primary residual states. jp1 = np1-1 if (kp.eq.2) jp1 = mxdim2-np1+2 end if np = np1-1-jout p = np if (np.lt.0) go to 32 nh = np1-1 - jin+iout nph = np+nh if (nph.lt.1) go to 32 npihi1 = np1 - max0(jnin-inout(kc),jnout(kp)) c ** Now set indices for the packed array of emission rates. do 30 nppi1=npilo1,npihi1 wr = 0. sig = 0. if (rho(np1,nppi1).le.0.) go to 30 mppi1 = nppi1 if (kp.gt.3) mppi1=mxdimp-nppi1+2 if (iso.gt.0) then jppi1 = nppi1 if (kp.eq.2) jppi1=mxdim2-nppi1+3 end if c ** Calculate residual nucleus state density and emission rates. nppi = nppi1-jpout(kp)-1 npnu = np - nppi nhpi = nppi1-1 - jpin + ipout(kc) nhnu = np1 - nppi1 - jnin + inout(kc) if (nh.eq.1) then v = ehole(nhnu+1) else if (nh.eq.2) then v = vfull end if c ** For the initial composite nucleus, allow auxilliary configurations c ** in the residual nucleus. These correpond to CN states where one c ** hole (or more than one), by chance, occupies a single particle c ** state next to the Fermi level and that, when the Fermi level c ** moves down during particle emission, is no longer a hole state c ** but a particle state. Thus the residual nucleus has one less c ** degree of freedom than the "main" residual configurations. j2lim = nhpi j4lim = nhnu if (kc.eq.1) then j2lim = min0(nhpi,nppi) j4lim = min0(nhnu,npnu) end if do 29 j2 = j2lim,nhpi do 29 j4 = j4lim,nhnu nhj = j2 + j4 h = nhj nphj = np + nhj rhol = 0. if (nphj.lt.1) go to 29 ph = nphj c ** Calculate finite well depth corrections to the residual c ** nucleus state density. fwd(1) = 1. fwd(2) = 1. if (nhj.le.0) go to 28 if (nh.gt.2) then go to 28 else if (nh.eq.2) then if (ee.le.vfull) go to 28 else if (ee.le.etest(nhnu+1)) go to 28 end if fwd(1) = well(nphj,nhj,ee,v,vfull) fwd(2) = well(nphj+1,nhj,ee,v,vfull) 28 epbar = ee/ph ehbar = epbar if (fwd(1).gt.0.) epbar = epbar * fwd(2) / fwd(1) if (nhj.gt.0) ehbar = (ee - p*epbar) / h rhol = omega(kc,nppi,j2,npnu,j4,ee,kp,i,0) * ct(kc,kp,i) * fwd(1) 29 sig = sig + rhol wr = prod(ne) * sig / rho(np1,nppi1) if (wr.lt.flow) wr=0. if (wr.gt.spill) wr=0. c ** If multiple preeq emission and T consvervation are both active, c ** then set the fraction of primary emission going to the lower T. if (iso.eq.0) go to 30 if (multi.eq.0) go to 30 if (kc.gt.1) go to 30 if (kp.gt.2) go to 30 if (np1.gt.mxdim2+1) go to 30 if (i.eq.1) then xiso(jp1,jppi1,ne) = 1. else x = wr + w(mkp,mp1,mppi1,ne) if (x.gt.0) xiso(jp1,jppi1,ne) = 1. - wr/x end if 30 w(mkp,mp1,mppi1,ne) = w(mkp,mp1,mppi1,ne) + wr 32 continue if (iso.le.0) go to 36 if (i.gt.1) go to 36 if (ct(kc,kp,2).le.0.) go to 36 i=2 go to 16 36 continue c ** Next stmts zero primary emission rates for states with nh=0 c ** with A_a>1. They are already zero for incident nucleons. c ** [Emission rates are calculated to allow generation of xiso c ** for secondary emission following direct stripping.] c ** To allow particle emission to begin at (p,h)=(Aa,0), c ** comment all stmts below up to (but NOT including) stmt 38. if (kc.gt.1) go to 38 if (jin.eq.1) go to 38 do 34 ne=2,neps1(1,kp) mp1 = ndwn + 1 mppi1 = npilo1 if (kp.gt.3) then mp1 = mxdimp - mp1 + 1 mppi1 = mxdimp - mppi1 + 2 end if 34 w(mkp,mp1,mppi1,ne) = 0. 38 return end c c******************************************************************* c function fact(x) c c Calculates factorials using Sterling's Approximation with a c correction factor for better accuracy at low x c c Corr'n factor changed (Feb 2007) from (12*x-1) to 12*x in denom. c This does slightly better at low x. c c Called from: SINGLE c pi = 3.14159265 e = 2.71828183 fact = 1. if (x.le.1.) go to 10 rat = (x/e)**x root = sqrt(2.*pi*x) corr = 1. + 1./(12*x) fact = rat * root * corr 10 return end c c********************************************************************* c subroutine nutra(kp,multi,mxdim2) c c Written 1979; revised for shell effects 1982; c Revised again 1990, 2000 and 2005. c c Phenomenological energy spectra for nucleon transfer. c Includes isospin conservation, collective pairing and finite c well depth effects. c c called from: PRECO c calls to: OMEGA, WELL c common /direct/ betalam(10), coul(7), elam(10), elast(51), fracel, 1 lam(10), lamax, scoll(51), sigbar(7), sigcni, siggeom(51), 2 snock(6,51), snutra(6,51), wcoll(10) common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3), 1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6), 2 nsigin, nphd(3,7,2), rzz, sigin(7,51) common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin, 1 xncrit(7,2) common /shell/ db(2), ds(7,2), epauli(7,14,2), dodd(7,2), fspan, 1 ghole(2), gpart(2), idif(7,2), ishell, mdnear(7,2), mdfar(7,2), 2 ndwn, nup2, range(2), rose(7), thresh(7,14,2) common /state/ fac(52), ze common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3), 1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff common /welld/ ehole(2),ehbar,epbar,etest(2), fwd(2), vfull common /wmain/ w(3,13,14,51), strong(9,10,51), seq(6,51), 1 seq2(2,51), steq(6,51), steq2(2,51), xiso(9,10,51), xisoeq(2,51), 2 spre(2,51) flow9 = 1.001e-09 spill = 1.e+18 xin = jin kin = jin + jpin xpin = jpin jout = jpout(kp) + jnout(kp) xout = jout atarg = acom(1) - xin ecom = e - ben(1,7) elabi = ecom * acom(1) / atarg ares = acom(1) - xout umax = e - ben(1,kp) vab = 0.25 * 50. * xin xkab = 12. if (jin.eq.4) then if(ecom.gt.20.) xkab = 12. - 11.*(ecom-20.)/ecom end if if (jin.eq.2) xkab=1. if (jin.eq.3) xkab=1. if (jout.eq.2) xkab=1. if (jout.eq.3) xkab=1. if (jin.eq.jout) xkab=1. c ** xkab is an enhancement factor for (nucln,He4) and (He4,nucln) velo = xin / (elabi+vab) x = jout * nsd(kp) * xkab y = 6. if (jin.eq.1) y=2. if (jin.eq.4) y=4. c ** y is Aproj*(2s+1) xhold = 0.0125 * x * ct(1,7,1) / (y*ecom*xin) if (jin.gt.jout) xhold = xhold * sqrt(ecom) / 7.25 c if (jin.eq.jout) xhold = xhold * sqrt(ecom) / 15.5 [WRONG may 06] if (jin.eq.jout) xhold = xhold * sqrt(ecom) / 14.5 c ** trann (tranp) is number of transferred neutrons (protons) tranp = iabs(jpout(kp)-jpin) za = (rzz+rzz)/atarg trann = iabs(jnout(kp)-jnin) ph = trann + tranp nph = ph nex = 0 if (nph.gt.0) go to 4 c ** nex = 1(2) allows for exchange of a proton(neutron) in the c ** inelastic channel for projectiles with A = 2 to 4. nex = 1 ph = 2. nph = 2 tranp = 2. 4 jp = jpout(kp) - jpin jn = jnout(kp) - jnin jpickp = max0(jp,0) jstrip = jpickp - jp jpickn = max0(jn,0) jstrin = jpickn - jn nppi = jstrip if (nex.eq.1) nppi=1 nhpi = jpickp if (nex.eq.1) nhpi=1 npnu = jstrin if (nex.eq.2) npnu=1 nhnu = jpickn if (nex.eq.2) nhnu=1 x = xhold * (3800./ares)**nph if (kin.eq.1) x = xhold * (5500./ares)**nph x = x * velo**(2.*ph) pow = 2.*(2.+xpin)*nhpi + npnu+npnu x = x * za**pow v = ehole(1) c ** Currently the only projectile for which there is a different c ** ehole for p and n excitation is neutrons. Pickup reactions c ** for incident neutrons must involve a picked up proton so c ** ehole(1) is the appropriate Veff to use. c ** c ** The factor 'reduce' (below) reduces the intensity of c ** configurations with an extra p-h pair in the residual nucleus. c ** Redux is set to this value when these configs are calculated. reduce = nppi*nppi + 1.5*nhpi*nhpi + npnu*npnu + nhnu*nhnu reduce = reduce * sqrt(elabi/xin) reduce = reduce * 7 / (v * atarg*atarg) if (nppi+nhpi.eq.0) v = v * za c ** Above stmt reduces the effective well depth for pure n transfer c ** when calculating finite well depth corrections in the c ** state densities. do 24 ne=2,neps1(1,kp) i = 1 u = umax-eps(ne) sig = 0. 8 ee = u-esym(1,kp,i) if (ee.le.0.) go to 22 c ** Next statements allow for extra configurations c ** (base + 1p1h or 2p2h). c ** nexmax is the maximum number of extra p-h pairs allowed. nexmax = 3 do 18 nexp = 0,nexmax kppi = nppi + nexp khpi = nhpi + nexp do 18 nexn = 0,nexmax-nexp kpnu = npnu + nexn khnu = nhnu + nexn kpart1 = kppi + kpnu + 1 next = nexp + nexn redux = 1. if (next.gt.0) redux = reduce**next c ** Below are new loops added to allow for additional residual c ** configurations corresponding to transfer at the Fermi level. if (next+nex.eq.0) then nppilo = 0 nhpilo = 0 npnulo = 0 nhnulo = 0 else nppilo = kppi nhpilo = khpi npnulo = kpnu nhnulo = khnu end if do 18 j1=nppilo,kppi do 18 j2=nhpilo,khpi do 18 j3=npnulo,kpnu do 18 j4=nhnulo,khnu np = j1 + j3 nh = j2 + j4 nphtemp = np + nh if (nphtemp.eq.0) go to 18 p = np h = nh ph = nphtemp c ** End of creation of new loops. do 12 ia=1,2 gpart(ia) = gres(1,kp,ia) 12 ghole(ia) = gres(1,kp,ia) fwd(1) = 1. fwd(2) = 1. if (nh.le.0) go to 14 c ** Finite well depth correction results are evaluated with veff c ** for the first interaction even though nh may be higher than 1. fwd(1) = well(nphtemp,nh,ee,v,vfull) fwd(2) = well(nphtemp+1,nh,ee,v,vfull) 14 epbar = ee/ph ehbar = epbar if (fwd(1).gt.0.) epbar=epbar*fwd(2)/fwd(1) if (nh.gt.0) ehbar = (ee-p*epbar)/h rhol = omega(1,j1,j2,j3,j4,ee,kp,i,1) * redux*fwd(1)*ct(1,kp,i) sig = rhol + sig c ** Next stmts store (complex,N) strength into strong(np1,nppi1,ne) c ** to make it available for secondary preeq emission. Only c ** strength to T< states is stored here. c ** Strength to T> states (if any) is included in repent. if (multi.eq.0) go to 18 if (kp.gt.2) go to 18 if (i.gt.1) go to 18 rhol = rhol * eps(ne) * sigin(kp,ne) * x rhol = rhol / (sigcni*ct(1,7,1)) eint = (eps(ne+1)-eps(ne-1)) * 0.5 if (kp.eq.1) then strong(kpart1,kppi+1,ne)= strong(kpart1,kppi+1,ne)+rhol*eint else np1 = mxdim2 - kpart1 + 1 nppi1 = mxdim2 - kppi + 1 strong(np1,nppi1,ne) = strong(np1,nppi1,ne) + rhol*eint end if 18 continue c ** Recyle to new residual isospin if appropriate. if (iso.le.0) go to 22 if (i.gt.1) go to 22 if (ct(1,kp,2).le.0.) go to 22 i = 2 go to 8 22 sig = sig * eps(ne) * sigin(kp,ne) * x if (sig.lt.flow9) sig=0. if (sig.gt.spill) sig=0. 24 snutra(kp,ne) = snutra(kp,ne) + sig if (kin.ne.kp) go to 26 c ** Next stmts change exchanged pair for inelastic channel. if (nex.eq.2) go to 26 nex = 2 tranp = 0. trann = 2. ph = 2. go to 4 26 return end c c******************************************************************** c function omega(kc,nppi,nhpi,npnu,nhnu,ee,kp,it,in) c c written in 1990 c c Calculate 2-component particle-hole state densities in the c shell-shifted equi-spacing model c normal pairing is included in epauli c collective pairing corrections added here c isospin conservation corrections also added here c finite well depth corrections added in calling routine c c called from: COMPOUND, NUTRA, EMMY c calls to: COLLEC, TCOR, SINGLE c c c * I am alpha and omega, the beginning and the end. c * I will give unto him that is athirst of the fountain of the c * water of life freely. (Rev. 21.6) c common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin, 1 xncrit(7,2) common /shell/ db(2), ds(7,2), epauli(7,14,2), dodd(7,2), fspan, 1 ghole(2), gpart(2), idif(7,2), ishell, mdnear(7,2), mdfar(7,2), 2 ndwn, nup2, range(2), rose(7), thresh(7,14,2) common /state/ fac(52), ze common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3), 1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff common /welld/ ehole(2),ehbar,epbar,etest(2), fwd(2), vfull pour = 1.e+35 c ** Set up particle and hole number related quantities nppi1 = nppi+1 nhpi1 = nhpi+1 npnu1 = npnu+1 nhnu1 = nhnu+1 ppi = nppi hpi = nhpi pnu = npnu hnu = nhnu omega = 0. nph = nppi + nhpi + npnu + nhnu if (nph.lt.1) go to 20 nq1 = max0(nppi1,nhpi1) nq2 = max0(npnu1,nhnu1) ph=nph c ** Main terms in the Pauli correction term coll = 0. if (icol.gt.0) coll = collec(kc,nppi,nhpi,npnu,nhnu,kp,in) adjust = 0. if (in.eq.0) go to 8 c ** IN=1 means omega is called from nutra. ADJUST (calculated c ** below) cancels out simple pairing corrections for direct c ** transfer (in nutra) of a pair of neutrons or protons, since c ** pairs are not broken and no pairing correction is needed in c ** that case. if (nppi+nhpi.eq.2) then if (nppi.ne.1) adjust = cpre(kc,kp,1) end if if (npnu+nhnu.eq.2) then if (npnu.ne.1) adjust = adjust + cpre(kc,kp,2) end if 8 hold = thresh(kp,nq1,1) + thresh(kp,nq2,2) + coll - adjust if (ee.lt.hold) go to 20 c ** Set single particle state densities if (ishell.gt.0) then call single(kc,nppi,nhpi,npnu,nhnu,ee,kp,coll,adjust) else do 10 ia=1,2 gpart(ia) = gres(kc,kp,ia) 10 ghole(ia) = gres(kc,kp,ia) end if c ** Calculate non-energy dependent factors c ** ACT is the number of active classes of excitons. rho = 1. act=0. if (nppi1.gt.1) then rho = rho * (gpart(1)**nppi) / fac(nppi1) act = 1. end if if (nhpi1.gt.1) then rho = rho * (ghole(1)**nhpi) / fac(nhpi1) act = act+1. end if if (npnu1.gt.1) then rho = rho * (gpart(2)**npnu) / fac(npnu1) act = act+1. end if if (nhnu1.gt.1) then rho = rho * (ghole(2)**nhnu) / fac(nhnu1) act = act+1. end if c ** The above tests are needed in the S2-ESM at low excitation c ** energies because the g's can be zero giving 0**0. c ** Now calculate the energy dependent part of the Pauli corr'n. gpi = gres(kc,kp,1) gnu = gres(kc,kp,2) qpi = nq1-1 qnu = nq2-1 edif = ee-hold term3 = 0. act1 = 1./act act2 = act1+act1 act3 = act2+act1 if(qpi.gt.0.) then den = 5.7 - 0.6*act den = den * gpi * edif * act / ph den = 9.35 + den t3 = 0. if (nppi.gt.0) t3 = (ppi-act2)*(ppi-act3) + act1 if (nhpi.gt.0) t3 = (hpi-act2)*(hpi-act3) + act1 + t3 term3 = t3/(gpi*den) end if if(qnu.gt.0.) then den = 5.7 - 0.6*act den = den * gnu * edif * act / ph den = 9.35 + den t3 = 0. if (npnu.gt.0) t3 = (pnu-act2)*(pnu-act3) + act1 if (nhnu.gt.0) t3 = (hnu-act2)*(hnu-act3) + act1 + t3 term3 = t3/(gnu*den) + term3 end if c ** Now put it all together ze = ee - epauli(kp,nq1,1)-epauli(kp,nq2,2) - coll+adjust - term3 if (ze.le.0) go to 20 omega = (ph-1.)*alog10(ze) - alog10(fac(nph)) omega = rho * 10.**omega if (omega.lt.0.1) omega=0. if (omega.gt.pour) omega=0. if (iso.gt.0) omega = omega * tcor(kc,ppi,hpi,pnu,hnu,ee,kp,it) 20 return end c c********************************************************************** c subroutine single(kc,nppi,nhpi,npnu,nhnu,e,kp,coll,adjust) c c Effective single particle state densities in the S2-ESM c c Written 1994 from stand-alone program shell.for c Replaces 1982 subroutine (revised in 1990,1992) c c called from: COMPOUND and OMEGA c dimension del(2), fe(2), fs(2), ish(2), nq(2), np(2), nh(2), 1 nhole(2), nhlo(2) common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin, 1 xncrit(7,2) common /shell/ db(2), ds(7,2), epauli(7,14,2), dodd(7,2), fspan, 1 ghole(2), gpart(2), idif(7,2), ishell, mdnear(7,2), mdfar(7,2), 2 ndwn, nup2, range(2), rose(7), thresh(7,14,2) c ------------------------------------------------------------ c 1. Initialize variables and calculate preliminary quantities c ------------------------------------------------------------ do 170 ia = 1,2 fe(ia) = 1. fs(ia) = 0. gpart(ia) = gres(kc,kp,ia) 170 ghole(ia) = gres(kc,kp,ia) npi = nppi + nhpi nnu = npnu + nhnu nph = npi+nnu ph = nph nqpi1 = max0(nppi,nhpi) + 1 nqnu1 = max0(npnu,nhnu) + 1 nq(1) = nqpi1 - 1 nq(2) = nqnu1 - 1 np(1) = nppi np(2) = npnu nh(1) = nhpi nh(2) = nhnu ihlotot=0 ihtot=0 ex = e - thresh(kp,nqpi1,1)-thresh(kp,nqnu1,2) - coll + adjust if (ex.lt.0.) go to 190 exfree = ex phfree = ph c ** Exfree is the permutable excitation energy, i.e. that not c ** tied up in exiting excitons across a shell gap. c ** Phfree is the number of excitons not trapped between the Fermi c ** level and a shell gap. c ** Both are refined in the next loop. do 172 ia = 1,2 nhole(ia) = nh(ia) nhlo(ia) = 0. del(ia) = 0. dbig = db(ia) dsmall = ds(kp,ia) c ** Now calculate the fraction of shell effects remaining. c ** If idifa.ge.range, shell effects are assumed to be washed out i = 0 if (dbig.gt.dsmall) then xmovf = range(ia) idifa = abs(idif(kp,ia)) difa = idifa if (difa.lt.xmovf) i=1 end if ish(ia) = i c ** Ish(ia)=0 means no shell effects; 1 means shell effects active if (i.eq.0) go to 172 if (nq(ia).le.0) go to 172 fe(ia) = difa / xmovf fs(ia) = 1. - fe(ia) c ** Fs is the fraction of shell effects still present c ** Fe is the fraction washed out with distance from shell closure if (idif(kp,ia).gt.0) then ih = nh(ia) else ih = np(ia) nhole(ia) = ih end if iq = nq(ia) ihlo = min0(idifa,iq) - iq + ih if (ihlo.le.0) go to 172 nhlo(ia) = ihlo c ** Ihlo and nhlo(ia) are the number of excitons of type ia c ** "trapped" between the Fermi level and a shell gap at threshold. ihlotot = ihlotot + ihlo ihtot = ihtot + ih c ** Ihlotot (& hlotot) are the total number of trapped excitons c ** (neutron+proton) at threshold. c ** Ihtot is the total number of excitons on the shell gap side of c ** the Fermi level when there are trapped excitons. hlo = ihlo dlim = (ex+dsmall*0.5) / dbig limdel = dlim limdel = amin0(limdel,ihlo) dlim = limdel c ** Dlim is the maximum number of holes which can be excited across c ** the shell gap. exn = ex * hlo / ph deltg = amin1(dlim,exn/dbig) if (ih.lt.nph) then if (ih.gt.ihlo) then x = amin1(difa*dsmall/dbig,0.5) deltg = deltg * (1.-x) end if end if del(ia) = deltg c ** Deltg=del(ia) is the average number of holes excited across c ** the gap. c ** Now calculate corrections to exfree and phfree when excitons c ** have been excited across a shell gap. exfree = exfree - deltg*dbig * fs(ia) phfree = phfree + (deltg-hlo) * fs(ia) 172 continue if (ish(1)+ish(2).le.0) go to 190 c ** Ish(1)+ish(2)=0 means no shell effects active in this nucleus hlotot = ihlotot exph = 0. if (phfree.gt.0.002) exph = exfree / phfree esm = 0. esm1 = 0. c -------------------------------------------------------- c 2. Generate correction to N_{Ei} = g0exf for gap crossings c -------------------------------------------------------- c ** g0exf = g_{i0} * (E-threshold) / (n-1); it is denoted N_{Ei} in c ** the code writeup where i is pi or nu for protons or neutrons. c ** Corrections are denoted Delta_{Na} and Delta_{Nb} in the c ** writeup but are esm1 and esm-esm1 here. Esm1 is stored c ** separately, because it is used alone in evaluating the "highly c ** constrained" cases. do 176 ia = 1,2 if (nq(ia).le.0) go to 176 ihlo = nhlo(ia) if (ihlo.le.0) go to 176 if (nph.le.1) go to 176 es1 = 0. esm2 = 0. hlo = ihlo ih = nhole(ia) h = ih ihhi = ih - ihlo hhi = ihhi hfree = hhi + del(ia) if (hfree.le.0.) go to 176 dbig = db(ia) dlim = (ex+ds(kp,ia)*0.5) / dbig c ** Dlim is the maximum number of excitons that could cross c ** the shell gap. scale = h/ph limesm = dlim * scale limesm = min0(limesm,ihlo) c ** Limesm is the maximum number scaled by this class' share c ** of total permutable excitation energy c ** and limited by initial number available. dlim = limesm x = 0. if (limesm.gt.0) then do 174 l=1,limesm y = l 174 x = x + 1./y x = x * dbig end if if (limesm.lt.ihlo) then y = (ex*scale - dlim*dbig) / (dlim+1.) x = x + y end if if (ihlotot+ihhi.lt.nph) then if (ihlotot+ihlo.gt.3) then difa = abs(idif(kp,ia)) y = 1. - difa*ds(kp,ia)*0.5/dbig scale = scale * y x = x * y end if end if x = x * gres(kc,kp,ia) / (hhi+1.) es1 = hlo * dbig / h es1 = amin1(es1,x) deltg = del(ia) if (hhi+deltg.gt.0.) then nd3 = mdfar(kp,ia) xnd3 = nd3 x = xnd3 - hhi - deltg - 1. if (x.gt.0.) esm2 = 0.5 * x * (hlo-deltg) / hlo if (nd3.le.ihhi+1) esm2 = - amin1(ex/(hlo*dbig),1.) end if esm = esm + (es1+esm2) * scale * fs(ia) esm1 = esm1 + es1*scale * fs(ia) 176 continue c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * c * Main loop to calculate xnum1 and xnum2 and the g's * c * * * * * * * * * * * * * * * * * * * * * * * * * * * * c c ----------------------------- c 1. Set preliminary quantities c ----------------------------- do 186 ia = 1,2 g0 = gres(kc,kp,ia) if (nq(ia).le.0) go to 186 dsmall = ds(kp,ia) ds2 = 0.5 * dsmall dbig = db(ia) db2 = 0.5 * dbig iq = nq(ia) q = iq ih = nhole(ia) ip = np(ia) + nh(ia) - ih p = ip h = ih iqp = iq - ip iqh = iq - ih idifa = abs(idif(kp,ia)) difa = idifa mnear = mdnear(kp,ia) mfar = mdfar(kp,ia) dnear = mnear dfar = mfar dod = dodd(kp,ia) c ** Xnum1 (xnum2) is the number of s.p. states accessible on the c ** same (opposite) side of the gap as the Fermi level in the S2ESM c ** Nd1,nd2,nd3 are numbers of active states in the corresponding c ** degenerate levels. c ** Xc1,xc2,xc3 are their contributions to xnum1 and xnum2. c ** Xc2 and xc3 both contribute to xnum2. c ** Eta1,eta2,eta3 are the average occupation numbers of the c ** degenerate levels. c ** Eta1t,eta2t,eta3t are the threshold occupation numbers, and c ** eta1l,eta2l,eta3l are the limiting values at high energy. c ** Delt is the maximum number of excitons that can be excited c ** into or out of a degenerate level. if (nph.gt.1) then c ---------------------------------------------------------------- c 2. Accessible single particle state with shells washed out. n>1 c ---------------------------------------------------------------- xeff = 1. / (ph-1.) g0ex = g0 * ex g0exf = g0ex * xeff xnum1e = p + g0exf - esm xnum2e = h + g0exf - esm if (ip.eq.1) xnum1e = xnum1e + esm*0.5 if (ih.eq.1) xnum2e = xnum2e + esm*0.5 if (xnum1e.lt.0.) xnum1e = 0. if (xnum2e.lt.0.) xnum2e = 0. if (ish(ia).le.0) then xnum1 = xnum1e xnum2 = xnum2e go to 182 end if c ------------------------------------------------------------------ c 3. Accessible S2ESM s.p. states on the Fermi level or particle c side of gap for n>1 -- xnum1 c ------------------------------------------------------------------ nd1 = max0(mnear-idifa-iqp,0) xnd1 = nd1 eta1t = amin0(ip,nd1) xnum1 = amax0(nd1,ip) + g0exf - esm xc1 = xnd1 if (nd1.gt.1) then exn = eta1t * exph g0exn = g0 * exn eta1l = amin1(p-1.,xnd1*0.5) delt = 0. if (ex.gt.dod-ds2) then dif = amax0(ip-nd1,0) b = dif + dod*g0 - 0.5 delt = 2.*g0exn + b*b delt = sqrt(delt) - b if (eta1l.lt.xnd1*0.5) then if (ip.lt.nph) delt = delt*0.75 end if end if eta1 = amax1(eta1t-delt,eta1l) if (eta1.ge.1.) then xcomb = fact(xnd1) / (fact(eta1) * fact(xnd1-eta1)) xcomb = xcomb**(1./eta1) xc1 = eta1*xcomb xnum1 = xnum1 - xnd1 + xc1 else if (ip.eq.1) then xnum1 = xnum1 - xnd1/(g0exf+ph) end if end if if (xnum1.lt.0.) xnum1=0. xnum1 = xnum1*fs(ia) + xnum1e*fe(ia) c ----------------------------------------------------------------- c 4. Accessible S2ESM s.p. states on the side of gap opposite the c Fermi level for n>1. -- xnum2 c ----------------------------------------------------------------- if (ih.le.0) go to 182 nd2 = 0 nd3 = 0 nd3t = 0 eta2t = 0. eta3t = 0. dbe = dbig ihlo = nhlo(ia) hlo = ihlo ihhi = ih - ihlo hhi = ihhi c ** hhi (hlo) is num of "holes" above (below) gap at thresh. c ------------------------------------------------------------- c 4a. Limiting values for the occupation of degenerate levels c and contributions to xnum2 from any ESM states below the gap. c ------------------------------------------------------------- if (iq.le.idifa-mnear) then c ** Last exciton at threshold is below 1st degenerate level xnum2 = h dbe = dod + (difa-dnear-1.-q) * dsmall dbe = amax1(dod,dbe) + dbig c ** dbe is an effective value of the shell gap including c ** effects from the odd spacing just outside it. etest1 = ex - dbe + dbig if (etest1.ge.-ds2) then xnum2 = idifa - iqh nd2 = mnear end if etest = etest1 - dbig if (etest.ge.-ds2) then xnum2 = xnum2 + dfar nd3 = mfar end if else if (iq.le.idifa) then c ** Last exciton at threshold is in 1st degenerate level xnum2 = idifa - iqh nd2 = min0(idifa-iqh,mnear) eta2t = ih - max0(idifa-iqh-mnear,0) etest = ex - dbig if (etest.ge.-ds2) then xnum2 = xnum2 + dfar nd3 = mfar end if else c ** Last exciton at threshold is across the shell gap xs = iq - idifa - mfar - 1 if (xs.ge.0.) dbe = dbe + dod + xs*dsmall nd2 = min0(idifa-iqh,mnear) nd2 = max0(nd2,0) eta2t = nd2 nd3 = mfar - max0(iqh-idifa,0) nd3t = nd3 xnd3t = nd3t eta3t = iq - max0(idifa,iqh) eta3t = amin1(eta3t,xnd3t) xnum2 = idifa - iqh + max0(iq-idifa,mfar) end if xnd2 = nd2 xnd3 = nd3 c ----------------------------------------------------------- c 4b. Average occupation numbers of degenerate levels 2 and 3 c ----------------------------------------------------------- limdel = (ex+ds2) / dbe limdel = min0(limdel,ihlo) deltc = limdel c ** Deltc is the maximum number of excitons of this class c ** that can be excited across the gap. (used 2 places below) exn = ex * hlo / ph deltg = amin1(deltc,exn/dbe) c ** Deltg calculated earlier assumed dbe = db(ia) if (ih.lt.nph) then if (ihhi.gt.0) then x = amin1(difa*dsmall/dbig,0.5) deltg = deltg * (1.-x) end if end if eta3l = amin1(xnd3*0.5,hhi+deltg-1.) eta3l = amax1(eta3l,0.) eta2l = amin1(h-1.-eta3l,xnd2*0.5) deltg = amin1(deltg,hlo-eta2l) xc3 = 0. eta3 = 0. if (nd3.gt.1) then eta3 = amin1(eta3t+deltg,xnd3) x = eta3t+deltg - eta3 if (eta3.gt.eta3l) then delt = 0. b = 0. if (ex.gt.dod-ds2) then dif = amax0(iq-idifa-mfar,0) exn = (eta3t+deltg) * exph c ** exph assumes that deltg is not reduced c ** for eta2l. exn = exn - dod*x - ds2*x*(dif+dif-1.+x) exn = amax1(exn,0.) b = dif+x + dod*g0 - 0.5 delt = 2. * g0 * exn + b*b delt = sqrt(delt) - b if (eta3l.lt.xnd3*0.5) then if (ih.lt.nph) delt = delt*0.75 end if end if eta3 = amax1(eta3-delt,eta3l) end if end if xc2 = 0. eta2 = 0. if (nd2.gt.1) then eta2 = eta2t - deltg if (eta2t.lt.eta2l) then c ** all "holes" are below the gap at threshold c ** add'l holes available below lower degen level delt = 0. if (ex.ge.dod-ds2) then g0exn = g0 * (h-eta2t) * exph b = g0*dod - 0.5 delt = sqrt(2.*g0exn+b*b) -b end if eta2 = amin1(eta2+delt,eta2l) else eta2 = amax1(eta2,eta2l) end if end if c ----------------------------------------------------- c 4c. Contributions from the degenerate levels to xnum2 c ----------------------------------------------------- if (nd3.gt.0) then xc3 = xnd3 if (eta3.ge.1.) then xcomb = fact(xnd3) / (fact(eta3)*fact(xnd3-eta3)) xcomb = xcomb**(1./eta3) xc3 = eta3*xcomb end if if (hhi+deltg.lt.1.) then xc3 = xnd3 * deltg if (ex.gt.dbe+ds2) xc3 = xc3 - xc3/(g0exf+ph) else if (hhi+deltg.lt.2) then if (ex.ge.(1.-hhi)*dbe+ds2) then xc3 = xc3 - xc3*(2.-hhi-deltg)/(g0exf+ph) end if end if xnum2 = xnum2 - xnd3 + xc3 end if if (nd2.gt.0) then if (eta2.gt.1.) then xcomb = fact(xnd2) / (fact(eta2)*fact(xnd2-eta2)) xcomb = xcomb**(1./eta2) xc2 = eta2*xcomb xnum2 = xnum2 - xnd2 + xc2 else xc2 = xnd2 end if c ------------------------------------------------ c 4d. Including effects of "Constraints" when most c excitons are trapped by gaps. c ------------------------------------------------ lcns1 = 0 lcns2 = 0 c ** lcns1 and lcns2 refer to two different types or c ** levels of constraint. limc = (ex-ds2)/dbe limdel = (ex+ds2)/dbe c ** Limc is limdel shifted up one s.p. spacing in energy if (ihhi+limc.eq.0) lcns1 = 1 c ** Lcns1=1 means there are no free excitons in this class if (ihtot.eq.nph) then if (ihhi.lt.2) then jhhi = nhole(3-ia) - nhlo(3-ia) if (jhhi.lt.2) then if (ihhi+jhhi.eq.2) limc=limdel jtest = 2 if (ih.lt.nph) jtest=4 jtest = min0(jtest,ih+1-jhhi) if (ihhi+limc.lt.jtest) lcns2 = 1 end if end if end if c ** Lcns2=1 means most excitons are trapped. lcnst = lcns1 + lcns2 c ** Lcnst=0 means there are no constraints if (lcnst.le.0) go to 180 if (ex.le.ds2) go to 180 if (limc.eq.limdel) then if (lcns1.eq.0) xnum2 = xnum2 - (g0exf-esm1)*0.5 xnum2 = xnum2 - xc2 - xc3 if (lcnst.eq.2) then xc2 = xc2 * 0.5 if (ihlotot.eq.nph) xc2 = 0. else ext = ex - (deltc-hhi)*dbe denom = dbe + hhi*dbe if (limc+ihhi.gt.1) then if (limc+ihhi.eq.jtest-1) then ext = ext + dbe denom = (denom+dbe)*2. end if end if x = 1. - 0.25*ext/denom xc2 = xc2 * x xc3 = xc3 * x end if xnum2 = xnum2 + xc2 + xc3 else if (ihtot.eq.nph) then cns = lcnst if(ihlotot.gt.1) cns= cns/(hlotot*(hlotot-1.)) xnum2 = xnum2 + (xc2+xc3)*cns if (lcnst.eq.2) then xnum2 = xnum2 - (g0exf-esm1)*0.5 else if (lcns2.eq.1) then if (ihhi+limdel.lt.jtest) then xnum2= xnum2- (g0exf-esm1)*0.5 end if end if end if 180 end if if (ihhi+limdel.gt.0) xnum2 = xnum2 + g0exf - esm if (xnum2.lt.0.) xnum2=0. xnum2 = xnum2*fs(ia) + xnum2e*fe(ia) c ----------------------------------------------------------- c 5. Convert xnum1 and xnum2 into effective g values; for n>1 c ----------------------------------------------------------- 182 g0exf = g0ex * xeff if (ish(ia).le.0) then if (ip.gt.0) gpart(ia) = g0 * xnum1 / (p+g0exf) if (ih.gt.0) ghole(ia) = g0 * xnum2 / (h+g0exf) else if (idif(kp,ia).gt.0) then c **i.e. if Z>Zmag or N>Nmag if (ip.gt.0) gpart(ia) = g0 * xnum1 / (p+g0exf) if (ih.gt.0) ghole(ia) = g0 * xnum2 / (h+g0exf) else if (ih.gt.0) gpart(ia) = g0 * xnum2 / (h+g0exf) if (ip.gt.0) ghole(ia) = g0 * xnum1 / (p+g0exf) end if else c ---------------------------------------------- c 6. Calculate xnum1 or xnum2 for case where n=1 c ---------------------------------------------- if (ish(ia).le.0) go to 186 if (idifa.lt.mnear) then if (ex.lt.dod-ds2) then xnum1 = dnear - difa else if (ex.lt.ds2) then xnum1 = dnear - difa + 1. else xnum1 = 1. end if else xnum1 = 1. end if xnum1 = xnum1*fs(ia) + fe(ia) if (idifa.le.0) then xnum2 = xnum1 else if (idifa.le.mnear) then if (ex.le.ds2) then xnum2 = difa else if (ex.lt.db2+db2-ds2) then xnum2 = 0. else if (ex.le.db2+db2+dod-ds2) then xnum2 = dfar else if (ex.le.db2+db2+ds2) then xnum2 = dfar + 1. else xnum2 = 1. end if else gs = difa * dsmall if (ex.lt.gs-dsmall-db2) then xnum2 = 1. else if (ex.le.gs-db2-dod) then xnum2 = dnear + 1. else if (ex.le.gs-db2) then xnum2 = dnear else if (ex.lt.gs-dsmall+db2) then xnum2 = 0. else if (ex.lt.gs-dsmall+db2+dod) then xnum2 = dfar else if (ex.le.gs+db2+dod) then xnum2 = dfar + 1. else xnum2 = 1. end if end if xnum2 = xnum2*fs(ia) + fe(ia) c ----------------------------------------------------- c 7. Convert xnum1, xnum2 to effective g value when n=1 c ----------------------------------------------------- if (idif(kp,ia).gt.0.) then if (ip.gt.0) gpart(ia) = g0*xnum1 if (ih.gt.0) ghole(ia) = g0*xnum2 else if (ih.gt.0) gpart(ia) = g0*xnum2 if (ip.gt.0) ghole(ia) = g0*xnum1 end if end if 186 continue 190 return end c c******************************************************************* c function tcor(kc,tppi,thpi,tpnu,thnu,ee,kp,i) c c written in 1990; revised June 1991 c c MAIN TERMS use type 2 blocking and state densities with their c natural particle and hole numbers c CORRECTION TERMS have no blocking and one state density c c called form: COMPOUND, OMEGA c calls to: ALPHA c common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin, 1 xncrit(7,2) common /state/ fac(52), zee common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3), 1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff flow = 1.e-18 spill = 1.e+18 tz2 = tze(kp)*2. e = ee + esym(kc,kp,i) - esym(kc,kp,i+1) c ---------------------------------------------------- c Set particle and hole numbers, c reversing the roles of protons and neutrons if Tz<0. c ---------------------------------------------------- if (itc(i).gt.0) then ppi = tppi hpi = thpi pnu = tpnu hnu = thnu gpi = gres(kc,kp,1) gnu = gres(kc,kp,2) else ppi = tpnu hpi = thnu pnu = tppi hnu = thpi gpi = gres(kc,kp,2) gnu = gres(kc,kp,1) end if ze = 0. denom = alpha(ppi,hpi,pnu,hnu,ee,ze,gpi,gnu) c ** Alpha returns a correct ESM value for ze=E-A(p,ppi) gpir = gpi/gnu gnur = gnu/gpi b = ppi*gnur + hnu*gpir cm = pnu + hpi c = cm + tz2 c ** The index i value transferred into TCOR is T-Tmin+1, c ** but TCOR needs T-abs(Tz)+1. Call this "it". c ** Here twot is 2Tmin delt = 0.5 * twot(kp) - abs(tz(kp)) it = delt+0.5 it = it + i c ---------------------------------------------------- c Calculate blocking fraction for type 2 isospin flips c ---------------------------------------------------- xi = it -1 ph = ppi + hpi + pnu + hnu sym = (tze(kp)+xi) * 4. / (gpi+gnu) zeb = ze - sym if (zeb.gt.0.) then f = (zeb/ze)**(ph-1.) if (f.lt.flow) f=0. if (f.gt.spill) f=0. f = 1.-f else f=1. end if c ----------------------------------------------------- c Calculate state densities according to T-Tz of states c ----------------------------------------------------- zet = ze + esym(kc,kp,i) - esym(kc,kp,i+1) if (zet.gt.0.) fac2 = (ph-1.) * (ph-2.) / (gpi*gnu * zet*zet) t2 = 0. if (it.eq.1) then t1 = denom if (zet.gt.0.) then bt = tz2+2. x1 = b / (c+3.) x2 = bt * ppi*hnu * fac2 / (c+2.) t2 = x1 + x2 end if else if (it.eq.2) then bt = tz2+2. bp0 = f*ppi bp1 = bp0 - f if (bp0.gt.bt) bp0=bt if (bp1.gt.bt) bp1=bt bh0 = f*hnu bh1 = bh0 - f if (bh0.gt.bt) bh0=bt if (bh1.gt.bt) bh1=bt b1 = bt - bp1 - bh0 + bp1*bh0/bt b2 = bt - bp0 - bh0 + bp0*bh0/bt b3 = bt - bp0 - bh1 + bp0*bh1/bt x1 = (pnu+1.) / (b1+cm+1.) x2 = 0. if (b2.gt.0.) x2 = b2 / (b2+cm) x3 = (hpi+1.) / (b3+cm+1.) t1 = x1 * alpha(ppi-1.,hpi,pnu+1.,hnu,ee,ze,gpi,gnu) t1 = t1 + x2 * denom t1 = t1 + x3 * alpha(ppi,hpi+1.,pnu,hnu-1.,ee,ze,gpi,gnu) if (zet.gt.0.) then x1 = x1 * ppi*gnur / (pnu+1.) x3 = x3 * hnu*gpir / (hpi+1.) bt = tz2+4. x11 = x1 * (b-gnur) / (c+6.) x12 = x1 * bt * (ppi-1.)*hnu * fac2 / (c+5.) x21 = x2 * b / (c+5.) x22 = x2 * bt * ppi*hnu * fac2 / (c+4.) x31 = x3 * (b-gpir) / (c+6.) x32 = x3 * bt * ppi*(hnu-1.) * fac2 / (c+5.) t2 = x11 + x12 + x21 + x22 + x31 + x32 end if else if (it.eq.3) then bt = tz2+4. bp0 = f*ppi bp1 = bp0 - f bp2 = bp1 - f if (bp0.gt.bt) bp0=bt if (bp1.gt.bt) bp1=bt if (bp2.gt.bt) bp2=bt bh0 = f*hnu bh1 = bh0 - f bh2 = bh1 - f if (bh0.gt.bt) bh0=bt if (bh1.gt.bt) bh1=bt if (bh2.gt.bt) bh2=bt b11 = bt - bp2 - bh0 + bp2*bh0/bt b12 = bt - bp1 - bh0 + bp1*bh0/bt b13 = bt - bp1 - bh1 + bp1*bh1/bt b22 = bt - bp0 - bh0 + bp0*bh0/bt b23 = bt - bp0 - bh1 + bp0*bh1/bt b33 = bt - bp0 - bh2 + bp0*bh2/bt den = b11+cm+2. x11 = (pnu+1.) * (pnu+2.) / (den * (den-1.)) den = b12+cm+1. x12 = (pnu+1.) * b12 / (den * (den-1.)) den = b13+cm+2. x13 = (pnu+1.) * (hpi+1.) / (den * (den-1.)) den = b23+cm+1. x23 = b23 * (hpi+1.) / (den * (den-1.)) den = b33+cm+2. x33 = (hpi+1.) * (hpi+2.) / (den * (den-1.)) den = b22+cm x22 = 0. if (den.gt.0.) x22 = b22 * (b22-1.) / (den * (den-1.)) t1 = x11 * alpha(ppi-2.,hpi,pnu+2.,hnu,ee,ze,gpi,gnu) y = alpha(ppi-1.,hpi,pnu+1.,hnu,ee,ze,gpi,gnu) t1 = t1 + (x12+x12) * y y = alpha(ppi-1.,hpi+1.,pnu+1.,hnu-1.,ee,ze,gpi,gnu) t1 = t1 + (x13+x13) * y t1 = t1 + x22 * denom y = alpha(ppi,hpi+1.,pnu,hnu-1.,ee,ze,gpi,gnu) t1 = t1 + (x23+x23) * y y = alpha(ppi,hpi+2.,pnu,hnu-2.,ee,ze,gpi,gnu) t1 = t1 + x33 * y if (zet.gt.0.) then x11 = x11 * ppi*(ppi-1.)*gnur*gnur / ((pnu+1.)*(pnu+2.)) x13 = x13 * ppi * hnu / ((pnu+1.)*(pnu+2.)) x33 = x33 * hnu*(hnu-1.)*gpir*gpir / ((hpi+1.)*(hpi+2.)) x12 = x12 * ppi*gnur / (pnu+1.) x23 = x23 * hnu*gpir / (hpi+1.) bt = tz2+6. x111 = x11 * (ppi-2.)*gnur / (c+9.) x113 = x11 * hnu*gpir / (c+9.) x133 = x13 * (hnu-1.)*gpir / (c+9.) x333 = x33 * (hnu-2.)*gpir / (c+9.) x = bt * fac2 / (c+8.) x112 = x11 * (ppi-2.)*hnu * x x132 = x13 * (ppi-1.)*(hnu-1.) * x x332 = x33 * ppi * (hnu-2.) * x x121 = x12 * (ppi-1.)*gnur / (c+8.) x123 = x12 * hnu*gpir / (c+8.) x233 = x23 * (hnu-1.)*gpir / (c+8.) x122 = x12 * bt * (ppi-1.)*hnu * fac2 / (c+7.) x232 = x23 * bt * ppi*(hnu-1.) * fac2 / (c+7.) x221 = x22 * b / (c+7.) x222 = x22 * bt * ppi*hnu * fac2 / (c+6.) t2 = x111 + x333 + x112 + x332 + x221 + x222 t2 = t2 + 2.*(x121 + x122 + x132 + x232 + x233) t2 = t2 + 3.*(x113 + x133) + 4.*x123 end if end if tcor = t1 if (zet.gt.0.) then tcor = tcor - t2 * alpha(ppi,hpi,pnu,hnu,e,ze,gpi,gnu) end if if (tcor.lt.0.) tcor=0. if (denom.gt.0.) then tcor = tcor / denom else tcor = 1. end if return end c c********************************************************************* c function well(n,nh,ee,v,vfull) c c Finite well depth correction functions for state densities c This includes averaging when an effective well depth due to c surface localization of the initial interaction is less than c the central well depth. c c called from: PRECO, NUTRA, EMMY c c * Therefore with joy shall ye draw waters out of the wells c * of salvation. (Isaiah 12.3) c c * But whosoever drinketh of the water that I shall give him c * shall never thirst, but the water that I shall give him shall c * be in him a well of water springing up into everlasting life. c * (John 4.14) c xx = 1. w = v * (vfull-v)/(vfull+vfull) if (n.eq.1) go to 8 c ------------------------------------------------------------ c For n>1 calculate the approximate average of the finite well c depth correction over a range of well depths if needed. c ------------------------------------------------------------ jmax = 5 if (v.ge.vfull-0.5) jmax=1 wtsum = 0. fwtsum = 0. ph1 = n-1 h = nh do 7 j=1,jmax xj = j-1 wt = (1.+exp(xj))*(1.+exp(-xj)) wt = 1./wt + 0.001*xj 2 vj = v + xj*w if (vj.gt.vfull) go to 6 if (vj.lt.0.) go to 6 f = 1. term = 1. xsign = 1. c ** Well is only called for h=1,2 in the exciton model, c ** but for 3 nucleon pickup or auxilliary configurations c ** in nutra, it can be called for higher n. c ** This programming is general. do 3 ih = 1,nh hi = ih ev = ee - hi*vj if (ev.le.0.01) go to 4 term = term * (h-hi+1.) / hi xsign = -xsign ev = alog10(term) + ph1*alog10(ev/ee) if (ev.lt.-5) go to 3 f = f + xsign * 10.**ev 3 continue 4 wtsum=wtsum+wt fwtsum=fwtsum+wt*f 6 if (xj.le.0.) go to 7 xj=-xj go to 2 7 continue if (wtsum.gt.0.) xx=fwtsum/wtsum go to 10 c ------------------------------------------ c Use a closed form average of fwd for n=h=1 c ------------------------------------------ 8 xx=0. if (ee.gt.1.16*vfull) go to 10 if (ee.lt.v) xx=1. if (v.ge.vfull) go to 10 xx=exp((ee-v)/w) xx=1./(1.+xx) 10 well=xx return end