C*******************************************************************************
C PROGRAM PLD (alias AVRIGEANU.FOR) *
C FORTRAN SUBROUTINES FOR CALCULATION OF PARTIAL LEVEL DENSITIES USED IN *
C PREEQUILIBRIUM NUCLEAR REACTION MODELS, FOR THE DEVELOPMENT OF THE *
C "REFERENCE INPUT PARAMETER LIBRARY" FOR NUCLEAR MODEL CALCULATIONS OF *
C NUCLEAR DATA (NUCLEAR DATA SECTION/IAEA-VIENNA PROJECT , 1994-1997). *
C *
C CARRIED OUT UNDER RESEARCH CONTRACT NO. 8886/R0-R1/RBF BETWEEN *
C INTERNATIONAL ATOMIC ENERGY AGENCY, VIENNA, AND *
C INSTITUTE FOR NUCLEAR PHYSICS AND ENGINEERING, BUCHAREST, ROMANIA. *
C CHIEF INVESTIGATOR: DR. M. AVRIGEANU, INPE-BUCHAREST. *
C SEPTEMBER 1997. *
C *
C OS/COMPILER: MS-DOS v.5/ MS-FORTRAN77 v.5.0 *
C*******************************************************************************
C
C I.DESCRIPTION OF THE SUBROUTINES/FUNCTIONS (* - 1997 new or major changes)
C ========================================
C
C 1. MAIN. * Reads the input data and organizes the calculation of the
C partial nuclear state densities (PSD) w(p,h,E) by using the specified
C formula [1-8], and partial level densities (PLD) D(p,h,E,J) by
C using the formalism of the spin cutoff factors of Fu [9,10]. It
C calculates (i) the nuclear state density w(E) as the PSD-sum over
C all allowed exciton numbers (p=h) for which PSD is higher than
C WMINACC=0.1 /MeV, and (ii) the total level density D(E) as the
C PLD-sum over the nuclear angular momentum J and allowed particle
C numbers p=h. At the same time are calculated the corresponding
C Wasym(E)- or Dasym(E)-values given by the asymptotical/'closed'
C formulae of the equidistant Fermi gas model (FGM) in order to make
C possible tests of the overall consistency. A two-fermion system
C correction [11] is involved optionally when the one-component
C Fermi gas formulae are used for calculation of PSDs.
C
C 2. PRINTIN. Prints the PSD/PLD formula and parameters used in calculation.
C
C 3. PRINTWN. Tabulates values of the calculated w(p,h,E) or D(p,h,E) (the
C latter being the sum over J of the D(p,h,E,J), i.e. the total
C level density for a given exciton configuration), for either
C (i) some given particle-hole configurations or (ii) all pairs of
C equal numbers of excited particles and holes. In the case (ii)
C there are printed also the w(E) or D(E)-values (within a first
C table including the PSD/PLD for p=h=1 to 7) and the corresponding
C Wasym(E)- or Dasym(E)-values given by the closed FGM formulae
C (within the second table including the PSD/PLD for p=h=8 to 16).
C In the case of the PLD calculation for a particular (p,h)
C configuration it is printed also a table of the D(p,h,E,J)-values
C but only for the last excitation energy E involved in the
C respective calculation; the total level density for this
C configuration is obtained and printed both as the D(p,h,E,J)-sum
C and 'D(p,h,E)-form' based on formula between this quantity and
C w(p,h,E).
C
C 4. WIL1 Calculates the partial nuclear state density w(p,h,E) of a given
C exciton configuration by means of the Williams formula [1] for
C one-component Fermi-gas.
C
C 5. WIL2 As WIL1 but using the two-fermion system formula.
C
C 6. WOB1 Calculates the partial nuclear state density of a given excited
C particle-hole configuration by means of the Betak-Dobes formula
C [2] with the nuclear potential finite-depth correction, for one-
C component Fermi-gas. If a value is specified for the nucleon
C binding energy it is carried out the calculation of the bound-
C state density according to Oblozinsky [3].
C
C 7. WOB2 As WOB1 but using the two-fermion system formula.
C
C 8. WFU1 Calculates the PSD of a given exciton configuration by means of
C the formula including the advanced pairing correction by Fu [4],
C for one-component Fermi-gas.
C
C 9. WFU2 * As WFU1 but using the two-fermion system formula.
C
C 10. WK1 Calculates the PSD of a given exciton configuration by using the
C improved implementation of the pairing correction by Kalbach [5,6]
C for one-component Fermi-gas.
C
C 11. WK2 * As WK1 but using the two-fermion system formula.
C
C 12. WK3 Calculates the PSD of a given exciton configuration by using the
C FGM energy-dependence of the single-excited particle and
C single-hole state densities, and/or the finite-depth correction
C including the nuclear-surface effects introduced by Kalbach [8]
C for one-component Fermi-gas.
C
C 13. WK4 * As WK3 but using the two-fermion system formula.
C
C 14. WM1 Calculates the PSD of a given exciton configuration by using the
C exact calculation of the Pauli-exclusion effect [7] and the
C pairing correction by Kalbach [5,6] for one-component Fermi-gas.
C
C 15. WM2 * As WM1 but using the two-fermion system formula.
C
C 16. WR1 * Calculates the PSD of a given exciton configuration by using
C (i) the improved implementation of the pairing correction by
C Kalbach [5,6], (ii) the FGM energy-dependence of the
C single-excited particle and single-hole state densities, and/or
C (iii) the finite-depth correction including the nuclear-surface
C effects introduced by Kalbach [8], for one-component Fermi-gas.
C
C 17. WR2 * As WR1 but using the two-fermion system formula.
C
C 18. PFU Calculates the advanced pairing correction by Fu [4], for the
C one-component Fermi-gas.
C
C 19. AK * Calculates the advanced pairing correction by Kalbach [5,6],
C for the one-component Fermi-gas.
C
C 20. FDC0 Calculates the nuclear potential finite-depth correction factor
C f(p,h,E,F) [8] for one-component Fermi-gas.
C
C 21. FDC * Calculates the nuclear potential finite-depth correction factor
C f(p+1,h,E,F) [8] in the case of bound states and the
C one-component Fermi-gas.
C
C 22. SUBPLD * Calculates the partial level density D(p,h,E,J) of a given
C 'p'-excited particle and 'h'-hole configuration for the excitation
C energy E and nuclear spin J, and the respective total level
C density D(p,h,E) as the sum over J of these PLDs, by using the
C partial state density w(p,h,E) and the formalism of Fu [9,10].
C
C 23. SIG2FU Calculates the spin cutoff factor for a given excited particle-
C hole configuration [9].
C
C 24. FCTR Calculates the factorial of natural numbers.
C
C
C II. REFERENCES
C ==========
C
C 1. F.C. Williams, Nucl. Phys. A166, 231 (1971)
C 2. E. Betak and J. Dobes, Z. Phys. A279, 319 (1976)
C 3. P.Oblozinsky, Nucl. Phys. A453, 127 (1986)
C 4. C.Y. Fu, Nucl.Sci.Eng. 86, 344 (1984)
C 5. C. Kalbach, Nucl.Sci.Eng. 95, 70 (1987)
C 6. C. Kalbach, Z.Phys. A 332, 157 (1989)
C 7. Mao Ming De and Guo Hua, J. Phys. G: 19, 421 (1993)
C 8. C. Kalbach, Phys. Rev. C 32, 1157 (1985)
C 9. C.Y. Fu, Nucl.Sci.Eng. 92, 440 (1986)
C 10. C.Y. Fu, Nucl.Sci.Eng. 109, 18 (1991)
C 11. J.M.Akkermans and H.Gruppelaar, Z.Phys. A 321, 605 (1985)
C
C
C III. INPUT DATA DESCRIPTION (Formatted read is used for all data in order to
C ====================== be possible the input of only few of them, while
C the rest are receiving zero-values)
C 1. NE,IOPTJ,TITLE
C ************** FORMAT(2I3,74A1)
C NE ..... Number of excitation energies.
C If NE>0 then the calculation is carried out for the energy values
C of E = 1, 2,.., NE MeV (maximum NEN=200), and record No. 2
C should be omitted
C IOPTJ .. Spin distribution (PLD calculation) option
C = 0 Calculate PSD-values, i.e. w(p,h,E) or w(pP,hP,pN,hN,E), and
C the nuclear state densities w(E) (obtained as sum over 'p' with
C the restriction p=h) and Wasym(E) given by the closed
C formulae
C = 1 Calculate the PLDs D(p,h,E,J) or D(pP,hP,pN,hN,E,J), the total
C level density given by the sum over J [ D(p,h,E) ] as well as
C over 'p' with the restriction p=h [ D(E) ], and Dform(E) given
C by the asymptotical/closed formula (Eq. (52) of [9])
C TITLE .. Title of the problem
C
C 2. E(I), I=1,|NE| (only if NE<0!) (it has to be omitted if NE>0!)
C ************** FORMAT(8F10.5)
C ... The excitation energies in MeV at which the PSD/PLDs are calculated
C
C 3. IMOD,ITFC,A,Z,UP
C **************** FORMAT(2I3,7F10.5)
C IMOD .. Select the partial state/level density formula (odd for one-fermion
C system and even for two-fermion system formulas):
C
C =-1, 0 Composite/recommended formula (present work)
C = 1, 2 F.C. Williams, Nucl. Phys. A166, 231 (1971)
C = 3, 4 P.Oblozinsky, Nucl. Phys. A453, 127 (1986), Eqs.(7,9)
C = 5, 6 C.Y. Fu, Nucl.Sci.Eng. 86, 344 (1984)
C = 7, 8 C. Kalbach, Nucl.Sci.Eng.95,70(1987),Z.Phys.A 332,157(1989)
C = 9,10 C. Kalbach, Phys.Rev. C 32, 1157 (1985)
C =11,12 Mao Ming De, J.Phys.G 19,421(1993)
C
C ITFC .. Option for two-fermion system correction
C = 0 No one
C = 1 J.M.Akkermans, H.Gruppelaar, Z.Phys. A 321,605(1985), Eq. (9)
C A ..... Mass number of the excited nucleus (may be omitted if GIN>0.)
C Z ..... Atomic number of the excited nucleus
C UP .... Pairing correction based on the odd-even mass differences
C If UP=-1. and A>0., P-values given by Eq.(9) of Dilg et al., Nucl.
C Phys. A217,269(1973), are adopted
C
C 4. NP0,NH0,GIN,FIN,BIN,F1IN
C ************************ FORMAT(2I3,7F10.5)
C NP0 ... Number 'p' of excited particles within the exciton configuration
C NH0 ... Number 'h' of holes within the exciton configuration
C If NP0=NH0=0 then the calculation is carried out for all pairs
C p=h=1, 2,... for which the PSD/PLD-value at the respective
C excitation energy E is higher than WMINACC=0.1 /MeV, followed by
C their sum to obtain the nuclear state density w(E), respectively
C the total level density D(E), and the values given by the
C asymptotic/closed formulae of the Fermi gas model
C GIN ... Input value of single-particle state density G in 1/MeV
C If GIN<0. the nuclear level density parameter DR(1) is read (record
C No. 6) and it is adopted the value G=(6/3.14**2)*DR(1)
C If GIN=0. and A=0. it is adopted the value G=1.0
C If GIN=0. and A>0. it is adopted the value G=A/13.
C FIN ... Input value of the Fermi energy F in MeV
C If FIN=0. it is adopted the value F=1.E+06
C BIN ... Input value of nucleon binding energy B in MeV
C If BIN=0. it is adopted the value B=1.E+06
C F1IN... Input value of the average effective Fermi energy F1 in MeV [8].
C If FIN=0.0 it is adopted the value F=1.E+06.
C If FIN<0.0 the PSD calculation by means of WR1 and WR2 functions
C is carried out for constant G.
C
C 5. NP0,NH0,NPN0,NHN0,GIN,GN,FIN,FN,BIN,BN,F1IN,F1N
C *********************************************** FORMAT(4I3,6F10.5)
C NP0 ... Number 'pP' of proton excited particles
C NH0 ... Number 'hP' of proton holes
C NPN0 .. Number 'pN' of neutron excited particles
C NHN0 .. Number 'hP' of neutron holes
C If NP0=NH0=NPN0=NHN0=0 the calculation is carried out for all
C configurations (pP=hP,pN=hN) - followed by increasing total exciton
C numbers N=pP+hP+pN+hN=2, 4,.. - for which the sum of PSD/PLDs for a
C given N at the respective excitation energy E is higher than
C WMINACC=0.1/MeV
C GIN ... Input value of single-proton state density G in 1/MeV
C GN .... Input value of single-neutron state density GN in 1/MeV
C If GIN<0 the nuclear level density parameter DR(1) is read (record
C No. 6) and there are adopted the values:
C G=Z/A*(6/3.14**2)*DR(1)
C GN=(A-Z)/A*(6/3.14**2)*DR(1)
C If GIN=0 and A=0. there are adopted the values G=GN=1.0
C If GIN=0 and A>0. are adopted the values G=Z/13 and GN=(A-Z)/13
C FIN ... Input value of the proton Fermi energy F in MeV
C FN .... Input value of the neutron Fermi energy FN in MeV
C If FIN=0. there are adopted the values F=FN=1.E+06
C BIN ... Input value of proton binding energy B in MeV
C BN .... Input value of neutron binding energy BN in MeV
C If BIN=0. there are adopted the values B=BN=1.E+06
C F1IN... Input value of the average effective proton Fermi energy F1 in MeV.
C If F1IN=0 it is adopted the value F1=1.E+06.
C If F1IN<0.0 the PSD calculation by means of the WR2 functions is
C carried out for constant G.
C F1N ... Input value of the average effective neutron Fermi energy F1N (MeV)
C If F1N=0 it is adopted the value F1N=1.E+06.
C If F1N<0.0 the PSD calculation by means of the WR2 functions is
C carried out for constant GN.
C
C 6. DR(K), K=1,3 (if GIN is negative, otherwise it has to be omitted!)
C ************ FORMAT(8F10.5)
C ... Parameters of the Back-Shifted Fermi Gas model for
C nuclear level density (W. Dilg et al., Nucl. Phys. A217,269(1973)]:
C DR(1) .. Nuclear level density parameter 'a'
C DR(2) .. Ratio of effective nuclear moment of inertia to rigid-body value
C calculated by using the reduced nuclear radius r=1.25 fm
C DR(3) .. Shift of the fictive nuclear ground state
C
C 7. ICONT,IEND
C ********** FORMAT(2I3,7F10.5)
C ICONT .. Output and recycle option
C =-1 Print of the first 2 tables only (see PRINTWN description)
C of calculated PSD/PLD, and resumption according to IEND
C = 0 Print of calculated PSD/PLD, and resumption according to IEND
C = 1 New calculation starting with input-data record 1, with the
C results printed together with those of the actual case
C = 2 New calculation starting with input-data record 3 (conserved
C energy grid)
C = 3 Calculation for another exciton configuration given by the
C input-data records either 4 or 5
C = 4 Calculation for another set of the BSFG model parameters given
C by the input-data record 6
C IEND ... Recycle option
C = 0 End
C = 1 New complete case starting with input-data record 1
C = 2 New calculation starting with input-data record 3
C = 3 New calculation for exciton configuration on records 4 or 5
C = 4 New calculation for BSFG model parameters given on record 6
C*******************************************************************************
c$debug
PARAMETER (NEN=200,NM=400,JM=50)
CHARACTER*1 TITLE(74)
REAL*8 WT,WTA,DT,SWK
DIMENSION D(NEN,NM),DT(NEN)
COMMON NPS(NM),NHS(NM),NPNS(NM),NHNS(NM),NMAX(NEN),
1 AVSIG2(NEN),DTA(NEN),E(NEN),TEF(NEN),
2 WT(NEN),WTA(NEN),W(NEN,NM),DPH(JM),PLDF,PSD,RSIG2,SIG2N
COMMON/CONST/A,Z,UP,G,F,F1,B,GN,FN,F1N,BN,GIN,FIN,F1IN,BIN,
1 IMOD,ICON,ITFC,IOPTJ,IPAR,TITLE,JMAX,DR(3),S
COMMON/UFU/DELTA2,ETH
EQUIVALENCE (D(1),W(1)),(DT(1),WT(1))
DATA R0/1.25/
C
OPEN(5,FILE=' ')
OPEN(6,FILE=' ')
WMINACC=.1
C
1 FORMAT(2I3,74A1)
2 FORMAT(8G10.3)
3 FORMAT(2I3,7F10.5)
4 FORMAT(4I3,8F10.5)
9 FORMAT(57X,'EXEC TIME=',F6.1,' S'/)
C
CALL TIME(T0)
ICONT=0
TEXT=0.
10 READ(5,1)NE,IOPTJ,TITLE
C ****
IF(NE.GT.0) THEN
DO 12 I=1,NE
12 E(I)=FLOAT(I)
ELSE IF(NE.EQ.0) THEN
GOTO 1000
ELSE
NE=IABS(NE)
READ(5,2)(E(I), I=1,NE)
C ****
ENDIF
L=1
C
20 READ(5,3)IMOD,ITFC,A,Z,UP
C ****
C
IF(Z.EQ.0) Z=FLOAT(IFIX(A)/2)
ICON=IABS(IMOD)-(IABS(IMOD)/2)*2
IPAR=NINT(A)-(NINT(A)/2)*2
IPARZ=NINT(Z)-(NINT(Z)/2)*2
UPN=0.
S=0.
ANC=0.
ANNC=0.
DR(1)=0.
SIG2EF=0.
IF(ICONT.NE.0) GOTO 30
DO 24 I0=1,NEN
WT(I0)=0.D+00
WTA(I0)=0.D+00
AVSIG2(I0)=0.
DO 24 L0=1,NM
24 W(I0,L0)=0.
C
30 IF(ICON.EQ.0) GOTO 35
READ(5,3)NP0,NH0,GIN,FIN,BIN,F1IN
C ****
G=GIN
GN=0.
F=FIN
B=BIN
F1=F1IN
TFC=1.
IF(A.EQ.0.AND.G.EQ.0) G=1.
IF(B.EQ.0.) B=1.E+06
IF(F.EQ.0.) F=1.E+06
IF(F1.EQ.0.) F1=1.E+06
IF(GIN.LT.0) GO TO 40
IF(G.EQ.0.) G=A/13.
GO TO 60
C
35 READ(5,4)NP0,NH0,NPN0,NHN0,GIN,GN,FIN,FN,BIN,BN,F1IN,F1N
C ****
G=GIN
F=FIN
B=BIN
IF(A.EQ.0.AND.G.EQ.0.) THEN
G=1.
GN=1.
ELSE
CONTINUE
ENDIF
F1=F1IN
IF(B.EQ.0.) B=1.E+06
IF(BN.EQ.0.) BN=1.E+06
IF(F.EQ.0.) F=1.E+06
IF(FN.EQ.0.) FN=1.E+06
IF(F1.EQ.0.) F1=1.E+06
IF(F1N.EQ.0.) F1N=1.E+06
IF(GIN.LT.0.) GO TO 40
IF(G.EQ.0.) G=Z/13.
IF(GN.EQ.0.) GN=(A-Z)/13.
IF(A.EQ.0.) THEN
IF(GN.EQ.0.) G=G/2.
IF(GN.EQ.0.) GN=G
ELSEIF(GIN.NE.0.) THEN
GN=G*(A-Z)/A
G=G*Z/A
ENDIF
GO TO 60
C
40 READ(5,2)(DR(K), K=1,3)
C ****
G=0.607927*DR(1)
IF(UP.EQ.-1.0.AND.A.GT.0.) THEN
IF(IPAR.EQ.1) THEN
UP=12.8/SQRT(A)
ELSE IF(IPARZ.EQ.1) THEN
UP=29.4/A
ELSE
UP=25.6/SQRT(A)
ENDIF
ELSE
CONTINUE
ENDIF
S=DR(3)-UP
C
60 CONTINUE
IF(ICONT.GE.2) L=L+1
IF(ICON.EQ.1) THEN
IF(NP0+NH0.GT.0) THEN
IOUT=1
ELSE
IOUT=2
ENDIF
ELSE IF(NP0+NH0+NPN0+NHN0.GT.0) THEN
IOUT=1
ELSE
IOUT=2
ENDIF
NMAX(1)=1
IF(GIN.GE.0.) DR(1)=(G+GN)/0.607927
UPP=UP/FLOAT(2-ICON) ! Mean-Gap Approximation (Fu'89)
IF(UP.GT.0.) THEN
ANC=1.584336*SQRT(G*UPP)
ANTC=1.584336*SQRT((G+GN)*UP)
SIGC2=0.274286*A**0.6666667*SQRT(G*UPP)
IF(ICON.EQ.1) THEN
ANNC=0.
SIGNC2=0.
ELSE
UPN=UP/2. ! Mean-Gap Approximation (Fu'89)
ANNC=1.584336*SQRT(GN*UPN)
SIGNC2=0.274286*A**0.6666667*SQRT(GN*UPN)
ENDIF
ELSE
CONTINUE
ENDIF
C
READ(5,3)ICONT,IEND
C ****
C End of the input data and exciton-configuration type setup
C Loop of the excitation energies
DO 600 I=1,NE
E1=E(I)
PEFFN=0.
GO TO (73,73,71,71,71,71,72,72,73,73,71,71,73,73), IMOD+2
71 E1EF=E1 ! IMOD=1-4,9,10
GO TO 80
72 E1EF=E1-UP-S ! IMOD=5-6
GO TO 80
73 IF(ANC.GT.0.) THEN ! IMOD=-1,0,7,8,11,12
IF(ANC.LE.4.48) THEN
c EK2=UPP*(1.+2.508/(ANC*ANC))
EK2=UPP
ELSE
EK2=UPP*(6.46/ANC-6.28/(ANC*ANC))
ENDIF
EK3=UPP/(1.+EXP(4.*(0.625-(E1/FLOAT(2-ICON))/UPP)))
PEFF=AMAX1(EK2,EK3)
IF(ANNC.GT.0.0) THEN
IF(ANNC.LE.4.48) THEN
c EK2=UPN*(1.+2.508/(ANNC*ANNC))
EK2=UPN
ELSE
EK2=UPN*(6.46/ANNC-6.28/(ANNC*ANNC))
ENDIF
EK3=UPN/(1.+EXP(4.*(0.625-(E1/2.)/UPN)))
PEFFN=AMAX1(EK2,EK3)
ELSE
CONTINUE
ENDIF
PEFF=(PEFF+PEFFN)/FLOAT(2-ICON)
IF(ICON.EQ.1)THEN
CONTINUE
ELSEIF(E1.LT.2.*UP)THEN
PEFF=PEFF*((E1-EK2)/(2.*UP-EK2)+1.)
ELSE
PEFF=PEFF*2.
ENDIF
E1EF=E1-PEFF-S
ELSE
E1EF=E1-S
ENDIF
C
C Initial/L-th exciton-configuration setup
80 NP=NP0+IOUT
NH=NH0+IOUT
IF(IOUT.EQ.2) L=1
NPS(L)=NP
NHS(L)=NH
IF(ICON.EQ.1) THEN
CONTINUE
ELSE
NPN=NPN0+1
NHN=NHN0+1
NPNS(L)=NPN
NHNS(L)=NHN
K=2
SWK=0.
ENDIF
EMAX=E1
WSSIG2=0.
SWI=0.
C
100 CONTINUE
P=NP-1
H=NH-1
IF(ICON.EQ.1) THEN
N=P+H
NN=0
ELSE
PN=NPN-1
HN=NHN-1
NN=PN+HN
N=P+H+PN+HN
ENDIF
C
IF(E1EF.GT.0.0) THEN
TEF(I)=0.5*(1.+SQRT(1.+ 4.*DR(1)*E1EF))/DR(1)
AMPN=1.3863*(G+GN)*TEF(I)
SIG2RB = 0.00957*A**1.666667*R0*R0*TEF(I)
SIG2EF = SIG2RB*DR(2)
IF(ICON.EQ.0.OR.ITFC.EQ.0) GOTO 102
TFC=1.02333*E1EF/(DR(1)**0.25*(E1EF+TEF(I))**1.25)
ELSE
CONTINUE
ENDIF
102 DO 105 JS=1,JM
105 DPH(JS)=0.
C
C Partial state density calculation
C
GO TO (297,298,110,120,130,140,150,160,170,180,190,200,210,220),
* IMOD+2
110 W(I,L)=TFC*WIL1(E1,NP,NH,G)
GO TO 300
120 W(I,L)=WIL2(E1,NP,NH,NPN,NHN,G,GN)
GO TO 300
130 W(I,L)=TFC*WOB1(E1,NP,NH,G,F,B)
Alpha=(P*P+H*H)/(2.*G)
EMAX=P*B+H*F-Alpha
GO TO 300
140 W(I,L)=WOB2(E1,NP,NH,NPN,NHN,G,GN,F,FN,B,BN)
Alphap=(P*P+H*H)/(2.*G)
Alphan=(PN*PN+HN*HN)/(2.*Gn)
Alpha=Alphap+Alphan
EMAX=P*B+H*F+PN*BN+HN*FN-Alpha
GO TO 300
150 W(I,L)=TFC*WFU1(E1,NP,NH,G,UP,S)
GO TO 300
160 W(I,L)=WFU2(E1,NP,NH,NPN,NHN,G,GN,UPP,UPN,S)
GO TO 300
170 W(I,L)=TFC*WK1(E1,NP,NH,G,UP)
GO TO 300
180 W(I,L)=WK2(E1,NP,NH,NPN,NHN,G,GN,UPP,UPN,S)
GO TO 300
190 W(I,L)=TFC*WK3(E1,NP,NH,G,F,F1)
GO TO 300
200 W(I,L)=WK4(E1,NP,NH,NPN,NHN,G,GN,F,FN,F1,F1N)
GO TO 300
210 W(I,L)=TFC*WM1(E1,NP,NH,G,F,B,UP)
EMAX=P*B+H*F-ETH
GO TO 300
220 W(I,L)=WM2(E1,NP,NH,NPN,NHN,G,GN,F,FN,B,BN,UPP,UPN,S)
EMAX=P*B+H*F+PN*BN+HN*FN-ETH
GO TO 300
297 W(I,L)=TFC*WR1(E1,NP,NH,G,F,F1,B,UP,S)
EMAX=P*B+H*F-ETH
GO TO 300
298 W(I,L)=WR2(E1,NP,NH,NPN,NHN,G,GN,F,FN,F1,F1N,B,BN,UPP,UPN,S)
EMAX=P*B+H*F+PN*BN+HN*FN-ETH
300 CONTINUE
C
C Partial level density calculation [C.Y.Fu, Nucl.Sci.Eng. 92,440(1986)]
IF(IOPTJ.EQ.0.OR.W(I,L).LT.WMINACC) THEN
PSD=0.
PLD=0.
PLDF=0.
SIG2N=0.
RSIG2=0.
GOTO 400
ELSE
PSD=W(I,L)
ENDIF
CALL SUBPLD(E1,E1EF,N,NN,ANC,ANNC,UPP,UPN,IPAR,SIGC2,SIGNC2,
* WMINACC,PSD,PLD,DPH,JMAX,SIG2,WSSIG2,SWI)
IF(SIG2.LE.0.) GOTO 400
D(I,L)=PLD
PLDF=PSD/(2.5066*SQRT(SIG2))
SIG2N=SIG2
IF(SIG2RB.GT.0.) RSIG2=SIG2N/SIG2RB
C
C Increase of exciton number for total exciton state/level density calculation
400 IF(IOUT.EQ.1) GO TO 500
IF(W(I,L).LT.WMINACC.AND.E1.LE.EMAX) THEN
W(I,L)=0.
IF(IMOD.NE.2.AND.E1.LT.N*F/2.) THEN
NMAX(I)=L-1
IF(N.GT.IFIX(1.25*AMPN)) GOTO 500
IF(ICON.EQ.1) GO TO 500
ELSE
CONTINUE
ENDIF
ELSE
WT(I)=WT(I)+W(I,L)
ENDIF
IF(L.EQ.NM) THEN
WRITE(6,*) ' EXCITON CONFIGURATIONS LIMITED AT N=',N
GOTO 500
ELSE
CONTINUE
ENDIF
L=L+1
IF(ICON.EQ.1) THEN
NP=NP+1
NH=NH+1
ELSE
SWK=SWK+W(I,L-1)
IF(NPN.LT.K) THEN
NP=NP-1
NH=NH-1
NPN=NPN+1
NHN=NHN+1
ELSE
IF(SWK.LT.WMINACC.AND.K.GT.3) GO TO 500
K=K+1
SWK=0.
NP=K
NH=K
NPN=NP0+1
NHN=NH0+1
ENDIF
NPNS(L)=NPN
NHNS(L)=NHN
ENDIF
NPS(L)=NP
NHS(L)=NH
GO TO 100
500 CONTINUE
C
C Calculation of closed forms (asymptotical) total state/level density
IF(IOUT.EQ.1.OR.E1EF.LE.0.) GOTO 600
IF(IOPTJ.GT.0) THEN
IF(SWI.GT.0.0.AND.WSSIG2.GT.0.) THEN
AVSIG2(I)=WSSIG2/SWI
FSIG=1./(2.5066*SQRT(AVSIG2(I)))
ELSE
FSIG=0.
ENDIF
ELSE
FSIG=1.
ENDIF
C
IF(ICON.EQ.1.AND.ITFC.EQ.0) THEN
IF(DR(1).LE.0.) THEN
WTA(I)=FSIG*EXP(3.141593*SQRT(2./3.*G*E1EF))/(SQRT(48.)*E1EF)
ELSE
WTA(I)=FSIG*EXP(2.*SQRT(DR(1)*E1EF))/(SQRT(48.)*E1EF)
ENDIF
ELSE
IF(DR(1).LE.0.) THEN
WTA(I)=(6.**0.25/12.)*EXP(3.141593*SQRT(2./3.*(G+GN)*E1EF))/
* ((G+GN)**0.25*E1EF**1.25)*FSIG
ELSE
WTA(I)=0.147704*EXP(2.*SQRT(DR(1)*E1EF))/
* (DR(1)**0.25*(E1EF+TEF(I))**1.25)*FSIG
ENDIF
ENDIF
600 CONTINUE
C End of excitation energy loop if any
C
GOTO (900,900,10,20,30,40), ICONT+2
C
C Print of the results for L-th exciton configurations at NE excitation energies
900 CALL PRINTIN
CALL PRINTWN(IOUT,L,NE,ICONT)
CALL TIME(T1)
TEX=(T1-T0)/100.
T0=T1
TEXT=TEXT+TEX
WRITE(6,9)TEX
ICONT=0
L=1
GOTO (1000,10,20,30,40), IEND+1
C
1000 WRITE(6,9)TEXT
STOP
END
C
SUBROUTINE PRINTIN
C*******************************************************************************
C PRINT OF THE NUCLEUS DATA *
C*******************************************************************************
CHARACTER*1 TITLE(74)
COMMON/CONST/A,Z,UP,G,F,F1,B,GN,FN,F1N,BN,GIN,FIN,F1IN,BIN,
1 IMOD,ICON,ITFC,IOPTJ,IPAR,TITLE,JMAX,DR(3),S
C
WRITE(6,1)TITLE
WRITE(6,2)Z,A,UP
IF(ICON.EQ.1) THEN
WRITE(6,3)
ELSE
WRITE(6,4)
ENDIF
IF(IMOD.EQ.-1.OR.IMOD.EQ.0) WRITE(6,10)
IF(IMOD.EQ.1.OR.IMOD.EQ.2) WRITE(6,11)
IF(IMOD.EQ.3.OR.IMOD.EQ.4) WRITE(6,13)
IF(IMOD.EQ.5.OR.IMOD.EQ.6) WRITE(6,15)
IF(IMOD.EQ.7.OR.IMOD.EQ.8) WRITE(6,17)
IF(IMOD.EQ.9.OR.IMOD.EQ.10) WRITE(6,19)
IF(IMOD.EQ.11.OR.IMOD.EQ.12) WRITE(6,21)
IF(ITFC.NE.0) WRITE(6,30)
IF(ICON.EQ.1) THEN
WRITE(6,40)G
IF(FIN.GT.0.0) WRITE(6,41)F
IF(F1IN.GT.0.) WRITE(6,42)F1
IF(BIN.GT.0.0) WRITE(6,43)B
IF(GIN.LT.0.0) WRITE(6,46) (DR(I), I=1,3),S
IF(GIN.EQ.0.AND.A.NE.0.) WRITE(6,48)
ELSE
WRITE(6,60)G,GN
IF(FIN.GT.0.0) WRITE(6,61)F,FN
IF(F1IN.GT.0.) WRITE(6,62)F1,F1N
IF(BIN.GT.0.0) WRITE(6,63)B,BN
IF(GIN.LT.0.0) WRITE(6,66)
IF(GIN.EQ.0.0.AND.A.NE.0.) WRITE(6,68)
ENDIF
C
1 FORMAT(1H ,'PARTIAL STATE/LEVEL DENSITIES CALCULATED FROM THE ',
1 'FOLLOWING PARAMETERS'/1H ,70('*')//1H ,74A1/)
2 FORMAT(2X,'NUCLEUS: CHARGE NO. Z=',F3.0,' MASS NO. A=',F4.0,
1 ' PAIRING Up=',F6.3,' MEV'/)
3 FORMAT(2X,'ONE-FERMION FORMULA:')
4 FORMAT(2X,'TWO-FERMION FORMULA:')
10 FORMAT(4X,'COMPOSITE FORMULA')
11 FORMAT(4X,'F.C. WILLIAMS, NUCL.PHYS. A166,231(1971)')
13 FORMAT(4X,'P.OBLOZINSKY, NUCL.PHYS. A453,127(1986), Eqs.(7,9)')
15 FORMAT(4X,'C.Y. FU, NUCL.SCI.ENG. 86,344(1984)')
17 FORMAT(4X,'C. KALBACH, NUCL.SCI.ENG. 95,70(1987),',
1 ' Z.PHYS. A 332,157(1989)')
19 FORMAT(4X,'C. KALBACH, PHYS.REV.C 32,1157(1985)')
21 FORMAT(4X,'MAO MING DE, GUO HUA, J.PHYS.G 19,421(1993)')
30 FORMAT(2X,'WITH TWO-FERMION SYSTEM CORRECTION:'
1 'AKKERMANS/GRUPPELAAR,Z.PHYS.A321,605(1985)')
40 FORMAT(/2X,'SINGLE-PARTICLE STATE DENSITY: G=',F6.3,' /MEV')
41 FORMAT(2X,'FERMI ENERGY: F=',F8.3,' MEV')
42 FORMAT(2X,'AV. EFF. FERMI ENERGY/1ST-2ND STEPS: F1=',F8.3,'MEV')
43 FORMAT(2X,'NUCLEON BINDING ENERGY: B=',F8.3,' MEV')
46 FORMAT(2X,'FROM: G=6/3.14**2*a; BSFG PARAMETERS: a=',F6.3,
1 '/MEV, I/Ir=',F4.2,', D=',F6.3,' MEV'/39X,'(S=',F6.3,' MEV)')
48 FORMAT(9X,'(PHENOMENOLOGICAL: G=A/13)')
60 FORMAT(/2X,'SINGLE-NUCLEON STATE DENSITIES: GP=',F8.3,
1 ' /MEV GN='F8.3' /MEV')
61 FORMAT(2X,'PROTON/NEUTRON FERMI ENERGIES: FP=',F8.3,' MEV',
1 ' FN=',F8.3,' MEV')
62 FORMAT(2X,'AV.EFF. FERMI EN./1ST-2ND STEPS: F1P=',F8.3,' MEV',
1 ' F1N=',F8.3,' MEV')
63 FORMAT(2X,'PROTON/NEUTRON BINDING ENERGIES: BP=',F8.3,' MEV',
1 ' BN=',F8.3,' MEV')
66 FORMAT(9X,'(FERMI-GAS MODEL: GP=3*Z/(2*FP), GN=3*(A-Z)/(2*FN)')
68 FORMAT(9X,'(PHENOMENOLOGICAL: GP=Z/13, GN=(A-Z)/13)')
RETURN
END
C
SUBROUTINE PRINTWN(IOUT,NL,NE,ICONT)
C*******************************************************************************
C PRINT OF THE TABLES OF CALCULATED PARTIAL-LEVEL DENSITIES *
C*******************************************************************************
PARAMETER (NEN=200,NM=400,JM=50)
CHARACTER*1 TITLE(74)
REAL*8 WT,WTA
COMMON NPS(NM),NHS(NM),NPNS(NM),NHNS(NM),NMAX(NEN),
1 AVSIG2(NEN),DTA(NEN),E(NEN),TEF(NEN),
2 WT(NEN),WTA(NEN),W(NEN,NM),DPH(JM),PLDF,PSD,RSIG2,SIG2N
COMMON/CONST/A,Z,UP,G,F,F1,B,GN,FN,F1N,BN,GIN,FIN,F1IN,BIN,
1 IMOD,ICON,ITFC,IOPTJ,IPAR,TITLE,JMAX,DR(3),S
C
S1=FLOAT(IPAR)/2.
JS1=S1+1.
WMINACC=1.
NEC=1
IF(E(NE).GT.E(1)) NEC=NE
DO 1 L=1,NMAX(NE)
IF(W(NEC,NMAX(NE)+1-L).GE.WMINACC) GO TO 2
1 CONTINUE
2 NMAX(NEC)=NMAX(NEC)+1-L
C
KF=1
NMIN=1
IF(IOUT.EQ.1) THEN
NMX=NL
NMX1=MIN0(8,NMX)
IF(NMX.GT.8) KF=(NMX-1)/8+1
ELSE
NMX=NMAX(NE)
NMX1=MIN0(7,NMX)
IF(NMX.GT.7) KF=(NMX-1)/7+1
IF(ICONT.LT.0) KF=2
ENDIF
C
DO 30 KE=1,KF
IF(KE.EQ.1) THEN
IF(IOUT.EQ.1) THEN
IF(ICON.EQ.1) THEN
IF(IOPTJ.LE.0) THEN
WRITE(6,130) (NPS(L)-1,NHS(L)-1, L=NMIN,NMX1)
ELSE
WRITE(6,132) (NPS(L)-1,NHS(L)-1, L=NMIN,NMX1)
ENDIF
ELSE
WRITE(6,140) (NPS(L)-1,NHS(L)-1,NPNS(L)-1,NHNS(L)-1,
* L=NMIN,NMX1)
ENDIF
ELSE
IF(ICON.EQ.1) THEN
IF(IOPTJ.LE.0) THEN
WRITE(6,150) (NPS(L)-1, L=NMIN,NMX1)
ELSE
WRITE(6,152) (NPS(L)-1, L=NMIN,NMX1)
ENDIF
ELSE
WRITE(6,170) (NPS(L)-1, NPNS(L)-1, L=NMIN,NMX1)
ENDIF
ENDIF
WRITE(6,188)
ELSE
IF(IOUT.EQ.1) THEN
NMIN=NMIN+8
NMX1=MIN0(NMX1+8,NMX)
ELSE
NMIN=NMIN+7
NMX1=MIN0(NMX1+7,NMX)
ENDIF
IF(KE.EQ.2) THEN
IF(IOUT.EQ.1) THEN
IF(ICON.EQ.1) THEN
IF(IOPTJ.LE.0) THEN
WRITE(6,130) (NPS(L)-1,NHS(L)-1, L=NMIN,NMX1)
ELSE
WRITE(6,132) (NPS(L)-1,NHS(L)-1, L=NMIN,NMX1)
ENDIF
ELSE
WRITE(6,140) (NPS(L)-1,NHS(L)-1,NPNS(L)-1,NHNS(L)-1,
* L=NMIN,NMX1)
ENDIF
ELSE
IF(ICON.EQ.1) THEN
IF(IOPTJ.LE.0) THEN
WRITE(6,180) (NPS(L)-1, L=NMIN,NMX1)
ELSE
WRITE(6,182) (NPS(L)-1, L=NMIN,NMX1)
ENDIF
ELSE
WRITE(6,185) (NPS(L)-1, NPNS(L)-1, L=NMIN,NMX1)
ENDIF
ENDIF
ELSE IF(IOUT.EQ.1) THEN
IF(ICON.EQ.1) THEN
IF(IOPTJ.LE.0) THEN
WRITE(6,130) (NPS(L)-1,NHS(L)-1, L=NMIN,NMX1)
ELSE
WRITE(6,132) (NPS(L)-1,NHS(L)-1, L=NMIN,NMX1)
ENDIF
ELSE
WRITE(6,140) (NPS(L)-1,NHS(L)-1,NPNS(L)-1,NHNS(L)-1,
* L=NMIN,NMX1)
ENDIF
ELSE IF(ICON.EQ.1) THEN
IF(IOPTJ.LE.0) THEN
WRITE(6,190) (NPS(L)-1, L=NMIN,NMX1)
ELSE
WRITE(6,192) (NPS(L)-1, L=NMIN,NMX1)
ENDIF
ELSE
WRITE(6,195) (NPS(L)-1, NPNS(L)-1, L=NMIN,NMX1)
ENDIF
WRITE(6,188)
ENDIF
C
DO 20 I=1,NE
IF(KE.EQ.1) THEN
IF(IOUT.EQ.1) THEN
WRITE(6,228) E(I), (W(I,L), L=NMIN,NMX1)
ELSE
WRITE(6,210) E(I),WT(I),(W(I,L), L=NMIN,NMX1)
ENDIF
ELSE IF(KE.EQ.2) THEN
IF(IOUT.EQ.1) THEN
WRITE(6,228) E(I), (W(I,L), L=NMIN,NMX1)
ELSE
WRITE(6,220) E(I),WTA(I),(W(I,L), L=NMIN,NMX1)
ENDIF
ELSE IF(IOUT.EQ.1) THEN
WRITE(6,228) E(I), (W(I,L), L=NMIN,NMX1)
ELSE
WRITE(6,230) E(I), (W(I,L), L=NMIN,NMX1)
ENDIF
20 CONTINUE
30 CONTINUE
C
IF(IOPTJ.GT.0) THEN
IF(IOUT.EQ.1) THEN
WRITE(6,300)
WRITE(6,200) E(NE),PSD,W(NE,NMX1),
* PLDF,SIG2N,RSIG2,(IJ+S1-1., DPH(IJ), IJ=JS1,JMAX)
ELSEIF(IOUT.EQ.2) THEN
WRITE(6,205) (E(I),AVSIG2(I), I=1,NE)
WRITE(6,206)
ENDIF
ENDIF
C
130 FORMAT(1H ,78('_')/1H ,6HENERGY,29X,9H w(p,h,E),/
1 1H ,6H (MeV),30X,8H (1/MeV),/1H ,78('_')
2 /1H ,6H(p,h)=,8(I2,1X,I2,4X))
132 FORMAT(1H ,78('_')/1H ,6HENERGY,29X,9H D(p,h,E),/
1 1H ,6H (MeV),30X,8H (1/MeV),/1H ,78('_')
2 /1H ,6H(p,h)=,8(I2,1X,I2,4X))
140 FORMAT(1H ,78('_')/1H ,6HENERGY,25X,17H w(pP,hP,pN,hN,E),/
1 1H ,6H (MeV),30X,8H (1/MeV),/1H ,6X,72('_')
2 /1H ,6H phph:,8(4I2,1X))
150 FORMAT(1H ,78('_')/1H ,6HENERGY,1X,8H w(E) ,27X,9H w(p,h,E),/
1 1H ,6H (MeV),1X,7H(1/MeV),29X,8H (1/MeV),/1H ,15X,63('_')
2 /1H ,15X,4Hp=h=,I1,3X,6(4X,I2,3X))
152 FORMAT(1H ,78('_')/1H ,6HENERGY,1X,8H D(E) ,27X,9H D(p,h,E),/
1 1H ,6H (MeV),1X,7H(1/MeV),29X,8H (1/MeV),/1H ,15X,63('_')
2 /1H ,15X,4Hp=h=,I1,3X,6(4X,I2,3X))
170 FORMAT(1H ,78('_')/1H ,6HENERGY,1X,8H w(E) ,23X,
1 17H w(P,hP,pN,hN,E),/
2 1H ,6H (MeV),1X,7H(1/MeV),29X,8H (1/MeV),/1H ,2X,76('_')
3 /1H ,2X,14H(pP=hP,pN=hN):,2I2,2X,6(3X,2I2,2X))
180 FORMAT(1H ,78('_')/1H ,6HENERGY,1X,8HWasym(E),27X,9H w(p,h,E),/
1 1H ,6H (MeV),1X,7H(1/MeV),29X,8H (1/MeV),/1H ,15X,63('_')
2 /1H ,15X,4Hp=h=,I1,3X,6(4X,I2,3X))
182 FORMAT(1H ,78('_')/1H ,6HENERGY,1X,8HDform(E),27X,9H D(p,h,E),/
1 1H ,6H (MeV),1X,7H(1/MeV),29X,8H (1/MeV),/1H ,15X,63('_')
2 /1H ,15X,4Hp=h=,I1,3X,6(4X,I2,3X))
185 FORMAT(1H ,78('_')/1H ,6HENERGY,1X,8HWasym(E),23X,
1 17H w(pP,hP,pN,hN,E),/
2 1H ,6H (MeV),1X,7H(1/MeV),29X,8H (1/MeV),/1H ,2X,76('_')
3 /1H ,2X,14H(pP=hP,pN=hN):,2I2,2X,6(3X,2I2,2X))
188 FORMAT(1H ,78('_'))
190 FORMAT(1H ,78('_')/1H ,6HENERGY,36X,9H w(p,h,E),/
1 1H ,6H (MeV),37X,8H (1/MeV),/1H ,15X,63('_')/
2 1H ,15X,4Hp=h=,I2,2X,6(4X,I2,3X))
192 FORMAT(1H ,78('_')/1H ,6HENERGY,36X,9H D(p,h,E),/
1 1H ,6H (MeV),37X,8H (1/MeV),/1H ,15X,63('_')/
2 1H ,15X,4Hp=h=,I2,2X,6(4X,I2,3X))
195 FORMAT(1H ,78('_')/1H ,6HENERGY,32X,17H w(pP,hP,pN,hN,E),/
1 1H ,6H (MeV),37X,8H (1/MeV),/1H ,2X,76('_')/
2 /1H ,2X,14H(pP=hP,pN=hN):,2I2,2X,6(3X,2I2,2X))
200 FORMAT(1H ,3H E=,F7.3,4H MEV,24X,17('_')/
1 2X,9Hw(p,h,E)=,6X,G11.4,6H /MEV;,5X,6H J ,11H D(p,h,E,J)/
2 2X,15HD(p,h,E,J)-sum=,G11.4,6H /MEV;,5X,6H(HBAR),2X,7H(1/MeV)/
3 2X,15HD(p,h,E)-form= ,G11.4,6H /MEV;,5X,17('_')/
4 2X,10HSIG2(E,n)=,G8.3,12HSIG2/SIG2RB=,G8.2,
5 F4.1,G11.4/(40X,F4.1,G11.4))
205 FORMAT(1H ,19('_')/1H ,6HENERGY,1X,12HAVERAGE-SIG2,/
1 1H ,6H (MeV),/1H ,19('_'),/(1X,F6.2,4X,G10.4))
206 FORMAT(1H ,19('_'))
210 FORMAT(1H ,F6.2,1X,G8.3,7(G9.3))
220 FORMAT(1H ,F6.2,1X,G8.3,7(G9.3))
228 FORMAT(1H ,F6.2, 8(G9.3))
230 FORMAT(1H ,F6.2,9X, 7(G9.3))
300 FORMAT(1H )
RETURN
END
C
SUBROUTINE SUBPLD(E1,E1EF,N,NN,ANC,ANNC,UP,UPN,IPAR,SIGC2,SIGNC2,
* WMINACC,PSD,PLD,DPH,JMAX,SIG2,WSSIG2,SWI)
C*******************************************************************************
C PARTIAL LEVEL DENSITY FORMULA: C.Y.FU,NUCL.SCI.ENG.92,440(1986) *
C*******************************************************************************
PARAMETER (JM=50)
DIMENSION DPH(JM)
IF(NN.EQ.0) THEN
SIG2=SIG2FU(E1,N,ANC,UP,SIGC2)
CONTINUE
ELSE
NP=N-NN
EP=E1*NP/N ! protons excitation energy
EN=E1*NN/N ! neutrons excitation energy
SIG2=SIG2FU(EP,NP,ANC,UP,SIGC2)
SIGN2=SIG2FU(EN,NN,ANNC,UPN,SIGNC2)
SIG2=SIG2+SIGN2
ENDIF
IF(SIG2.LE.0) GOTO 400
C97
IF(E1/UP.LT.1.5) SIG2NM=AMAX1(SIG2,0.345*SIGC2*(E1/UP)**0.927)
WSSIG2=WSSIG2+PSD*SIG2
SWI=SWI+PSD
PLD=0.
S1=FLOAT(IPAR)/2.
JS1=S1+1.
SJMAX=SQRT(2.*SIG2*E1EF)
SJMAX=AMAX1(10.,SJMAX)+1.
JMAX=MIN0(IFIX(SJMAX),JM)
SJ=S1
C
DO 320 JS=JS1,JMAX
DPH(JS)=PSD *
1 (2.*SJ+1.)/(5.013257*SIG2**1.5)*EXP(-(SJ+0.5)**2./(2.*SIG2))
IF(DPH(JS).LT.WMINACC) GOTO 330
PLD=PLD+DPH(JS)
SJ=SJ+1.
320 CONTINUE
C
330 JMAX=JS
400 RETURN
END
C
FUNCTION SIG2FU(E1,N,ANC,UP,SIGC2)
C*******************************************************************************
C SPIN CUTOFF FACTOR FOR STATES OF FIXED N:C.Y.FU,NUCL.SCI.ENG.92,440(1986) *
C*******************************************************************************
SIG2FU=0.
RNC=FLOAT(N)/ANC
IF(RNC.LE.0.446) THEN
UTH=UP*(3.23*RNC-1.57*RNC**2)
ELSE
UTH=UP*(1.+0.627*RNC**2)
ENDIF
IF(UTH.GE.E1) GOTO 10
X=-0.413+1.08*SQRT(RNC)-0.226*RNC
SIG2FU=SIGC2*1.38629*RNC*(1.-UTH/E1)**X
10 RETURN
END
C
REAL*8 FUNCTION FCTR(NN)
C***********************************************************************
C FACTORIAL CALCULATION *
C***********************************************************************
REAL*8 FACTOR
FACTOR=1.D+00
C
IF(NN.EQ.0) THEN
FCTR=1.
GO TO 10
ELSE
CONTINUE
ENDIF
DO 5 I=1,NN
FI=I
FACTOR=FACTOR*FI
5 CONTINUE
FCTR=FACTOR
10 RETURN
END
C
SUBROUTINE TIME(T)
C*******************************************************************************
C CALL DATE AND TIME *
C*******************************************************************************
CPC FOR MS-F77 v.5.0
CALL GETDAT(IYR,IMON,IDAY)
CALL GETTIM(IHR,IMIN,ISEC,I100TH)
WRITE(6,1)IYR,IMON,IDAY,IHR,IMIN,ISEC,I100TH
1 FORMAT(57X,I4,'-',I2,'-',I2,1X,I2,':',I2,':',I2,':',I2)
T=(IHR*3600.+IMIN*60.+FLOAT(ISEC))*100.+FLOAT(I100TH)
CPC
CVAX FOR VAX/VMS v.5.4-2 VAX FORTRAN77
CVAX CHARACTER*8 TSTR
CVAX CALL IDATE(IMON,IDAY,IYR)
CVAX CALL TIME(TSTR)
CVAX WRITE(6,1)IYR,IMON,IDAY,TSTR
CVAX1 FORMAT(56X,I4,'-',I2,'-',I2,1X,A8)
CVAX
RETURN
END
C
FUNCTION WIL1(E1,NP0,NH0,G)
C***********************************************************************
C ONE-FERMION FORMULA: F.C. WILLIAMS, NUCL.PHYS. A166,231(1971) *
C***********************************************************************
REAL*8 FCTR,T1,T2
NP=NP0-1
NH=NH0-1
P=FLOAT(NP)
H=FLOAT(NH)
Aph=(P*(P+1.)+H*(H-1.))/(4.*G)-H/(2.*G)
Ecor=E1-Aph
IF(Ecor.LE.0.)Ecor=0.
T1= DFLOAT(G**(P+H))/FCTR(NP+NH-1)
T2= FCTR(NP)*FCTR(NH)
WIL1=SNGL ( T1/T2*Ecor**(P+H-1.) )
RETURN
END
C
FUNCTION WIL2(E1,NPP0,NHP0,NPN0,NHN0,GP,GN)
C***********************************************************************
C TWO-FERMION FORMULA: F.C. WILLIAMS, NUCL.PHYS. A166,231(1971) *
C***********************************************************************
REAL*8 FCTR,T1,T2
NPP=NPP0-1
NHP=NHP0-1
NPN=NPN0-1
NHN=NHN0-1
NN=NPP+NHP+NPN+NHN
PP=FLOAT(NPP)
HP=FLOAT(NHP)
PN=FLOAT(NPN)
HN=FLOAT(NHN)
Ap=(PP*(PP+1.)+HP*(HP-3.))/(4.*GP)
An=(PN*(PN+1.)+HN*(HN-3.))/(4.*GN)
Apauli=Ap+An
Ecor=E1-Apauli
IF(Ecor.LE.0.)Ecor=0.
T1=DFLOAT(GP**(PP+HP))/(FCTR(NPP)*FCTR(NHP))
T2=DFLOAT(GN**(PN+HN))/(FCTR(NPN)*FCTR(NHN))
WIL2=SNGL( T1*T2* ( DFLOAT((Ecor**(NN-1)))/FCTR(NN-1)) )
RETURN
END
C
FUNCTION WOB1(E1,NP0,NH0,G,F,B)
C*******************************************************************************
C ONE-FERMION FORMULA: P.OBLOZINSKY, NUCL.PHYS. A453,127(1986), EQS.(7,9)*
C*******************************************************************************
REAL*8 FCTR,T1,T2,SUM
C
WOB1=0.
NP=NP0-1
NH=NH0-1
NN=NP+NH
P=FLOAT(NP)
H=FLOAT(NH)
T1=DFLOAT(G**(P+H))/FCTR(NN-1)
T2=FCTR(NP)*FCTR(NH)
CIN Alpha=(P*(P+1.)+H*(H-1.))/(2.*G)
CIN Aph=(P*(P+1.)+H*(H-3.))/(4.*G)
Alpha=(P*P+H*H)/(2.*G)
Aph=(P*(P-1.)+H*(H-1.))/(4.*G)
IF((E1+Alpha-P*B-H*F).GT.0.) GO TO 30
C
SUM=0.
DO 20 II=1,NP+1
I=II-1
DO 20 JJ=1,NH+1
J=JJ-1
Ecor=E1-Aph-I*B-J*F
Ecor1=E1-Alpha-I*B-J*F
IF(Ecor1.LE.0.)GO TO 20
D=(-1.)**(I+J)
CP=SNGL( FCTR(NP)/(FCTR(NP-I)*FCTR(I)) )
CH=SNGL( FCTR(NH)/(FCTR(NH-J)*FCTR(J)) )
SUM=SUM+D*CP*CH*SNGL(DFLOAT(Ecor**(NN-1))/T2)
20 CONTINUE
C
WOB1=SNGL(T1*SUM)
30 RETURN
END
C
FUNCTION WOB2(E1,NPP0,NHP0,NPN0,NHN0,GP,GN,FP,FN,BP,BN)
C*******************************************************************************
C TWO-FERMION FORMULA, CFM.: P.OBLOZINSKY, NUCL.PHYS. A453,127(1986) *
C*******************************************************************************
REAL*8 FCTR,T1,T2,SUM
C
WOB2=0.
NPP=NPP0-1
NHP=NHP0-1
NPN=NPN0-1
NHN=NHN0-1
NN=NPP+NHP+NPN+NHN
PP=FLOAT(NPP)
HP=FLOAT(NHP)
PN=FLOAT(NPN)
HN=FLOAT(NHN)
T1=DFLOAT(Gp**(PP+HP))/(FCTR(NPP)*FCTR(NHP))
T2=DFLOAT(Gn**(PN+HN))/(FCTR(NPN)*FCTR(NHN)*FCTR(NN-1))
Alphap=(PP*PP+HP*HP)/(2.*Gp)
Alphan=(PN*PN+HN*HN)/(2.*Gn)
Alpha=Alphap+Alphan
Ap=(PP*(PP-1.)+HP*(HP-1.))/(4.*Gp)
An=(PN*(PN-1.)+HN*(HN-1.))/(4.*Gn)
Apauli=Ap+An
IF((E1+Alpha-PP*BP-HP*FP-PN*BN-HN*FN).GT.0.) GO TO 30
C
SUM=0.
DO 20 IIp=1,NPP+1
Ip=IIp-1
DO 20 JJp=1,NHP+1
Jp=JJp-1
DO 20 IIn=1,NPN+1
In=IIn-1
DO 20 JJn=1,NHN+1
Jn=JJn-1
C
Ecor=E1-Apauli-Ip*Bp-In*Bn-Jp*Fp-Jn*Fn
Ecor1=E1-Alpha-Ip*Bp-In*Bn-Jp*Fp-Jn*Fn
IF(Ecor1.LE.0.)GO TO 20
D=(-1.)**(Ip+In+Jp+Jn)
CPP=SNGL( FCTR(NPP)/(FCTR(NPP-Ip)*FCTR(Ip)) )
CHP=SNGL( FCTR(NHP)/(FCTR(NHP-Jp)*FCTR(Jp)) )
CPN=SNGL( FCTR(NPN)/(FCTR(NPN-In)*FCTR(In)) )
CHN=SNGL( FCTR(NHN)/(FCTR(NHN-Jn)*FCTR(Jn)) )
SUM=SUM+D*CPP*CHP*CPN*CHN*DFLOAT((Ecor**(NN-1)))
20 CONTINUE
WOB2=SNGL(T1*T2*SUM)
30 RETURN
END
C
FUNCTION WFU1(E1,NP0,NH0,G,UP,S)
C*******************************************************************************
C ONE-FERMION FORMULA: C.Y. FU, NUCL.SCI.ENG. 86, 344 (1984) *
C*******************************************************************************
REAL*8 FCTR,T1,T2
COMMON/UFU/DELTA2,UTH
C
WFU1=0.
NP=NP0-1
NH=NH0-1
NN=NP+NH
P=FLOAT(NP)
H=FLOAT(NH)
T1=G**(P+H)/FCTR(NN-1)
T2=FCTR(NP)*FCTR(NH)
Aph=(P*(P+1.)+H*(H-3.))/(4.*G)
C
P1=PFU(E1,NN,G,UP)
IF(UTH.GT.E1) GO TO 10
Bph=Aph*SQRT(1.+4.*G*G*DELTA2/FLOAT(NN*NN))
Ecor=AMAX1(0., E1-P1-Bph-S )
WFU1=SNGL( T1*(DFLOAT(Ecor**(NN-1))/T2) )
10 RETURN
END
C
FUNCTION PFU(E1,NN,G,UP)
C*******************************************************************************
C PAIRING CORRECTION *
C ONE-FERMION FORMULA: C.Y. FU, NUCL.SCI.ENG. 86, 344 (1984) *
C*******************************************************************************
COMMON/UFU/DELTA2,UTH
IF(UP.LE.0.) THEN
PFU=0.
DELTA2=0.
UTH=0.
ELSE
ANC=1.584336*SQRT(G*UP)
RNC=FLOAT(NN)/ANC
IF(RNC.GE.0.446) THEN
UFUL=UP*(0.716+2.44*RNC**2.17)
UTH=UP*(1.+0.627*RNC**2)
ELSE
UFUL=0.
UTH=UP*(3.23*RNC-1.57*RNC**2)
ENDIF
IF(E1.GE.UFUL)THEN
R=0.996-1.76*RNC**1.60/(E1/UP)**0.68
ELSE
R=0.
ENDIF
PFU=UP*(1.-R*R)
DELTA2=4.*UP/G*R*R
ENDIF
10 RETURN
END
C
FUNCTION WFU2(E1,NPP0,NHP0,NPN0,NHN0,GP,GN,UP,UPN,S)
C*******************************************************************************
C TWO-FERMION FORMULA, CFM.: C. Y. FU, NOV.1989 *
C*******************************************************************************
REAL*8 FCTR,T1,T2
COMMON/UFU/DELTA2,UTH
C
WFU2=0.
NPP=NPP0-1
NHP=NHP0-1
NPN=NPN0-1
NHN=NHN0-1
NP=NPP+NHP ! total proton-exciton number
NN=NPN+NHN ! total neutron-exciton number
PP=FLOAT(NPP)
HP=FLOAT(NHP)
PN=FLOAT(NPN)
HN=FLOAT(NHN)
EP=E1*NP/(NN+NP) ! protons excitation energy
EN=E1*NN/(NN+NP) ! neutrons excitation energy
AP=(PP*(PP+1.)+HP*(HP-3.))/(4.*GP)
AN=(PN*(PN+1.)+HN*(HN-3.))/(4.*GN)
C
GPM=(GP+GN)/2.
IF(EP.GT.0.) THEN
PFUP=PFU(EP,NP,GPM,UP) ! Mean-Gap Approximation
DP2=DELTA2 ! pairing gap for protons
UTHP=UTH
ELSE
PFUP=0.
DP2=0.
UTHP=0.
ENDIF
IF(EN.GT.0.) THEN
PFUN=PFU(EN,NN,GPM,UPN) ! Mean-Gap Approximation
DN2=DELTA2 ! pairing gap for protons
UTHN=UTH
ELSE
PFUN=0.
DN2=0.
UTHN=0.
ENDIF
PFU2=PFUP+PFUN
UTH2=UTHP+UTHN
IF(UTH.GT.E1) GO TO 10
C
IF(NP.EQ.0.AND.NN.NE.0)THEN
BP=0.
BN=AN*SQRT(1.+4.*GN*GN*DN2/FLOAT(NN*NN)) ! new Pauli correction
ELSEIF(NP.NE.0.AND.NN.EQ.0.)THEN
BP=AP*SQRT(1.+4.*GP*GP*DP2/FLOAT(NP*NP)) ! new Pauli correction
BN=0.
ELSE
BP=AP*SQRT(1.+4.*GP*GP*DP2/FLOAT(NP*NP)) ! new Pauli correction
BN=AN*SQRT(1.+4.*GN*GN*DN2/FLOAT(NN*NN)) ! new Pauli correction
ENDIF
A2=BP+BN
Ecor=AMAX1(0.,E1-PFU2-A2-S)
T1=DFLOAT(GP**(PP+HP))/(FCTR(NPP)*FCTR(NHP))
T2=DFLOAT(GN**(PN+HN))/(FCTR(NPN)*FCTR(NHN))
WFU2=SNGL( T1*T2* ( DFLOAT((Ecor**(NP+NN-1)))/FCTR(NP+NN-1)) )
10 RETURN
END
C
FUNCTION WK1(E1,NP0,NH0,G,UP)
*******************************************************************************
C ONE-FERMION FORMULA: C. KALBACH, NUCL.SCI.ENG. 95,70(1987) *
C*******************************************************************************
REAL*8 FCTR,T1,T2
COMMON/UFU/DELTA2,ETH
C
WK1=0.
NP=NP0-1
NH=NH0-1
NN=NP+NH
P=FLOAT(NP)
H=FLOAT(NH)
T1=DFLOAT(G**(P+H))/FCTR(NN-1)
T2=FCTR(NP)*FCTR(NH)
C
AK1=AK(E1,NP,NH,NN,G,UP)
IF(ETH.GT.E1) GO TO 10
Ecor=AMAX1(0., E1-AK1 )
WK1=SNGL( T1*(DFLOAT(Ecor**(NN-1))/T2) )
10 RETURN
END
C
FUNCTION AK(E1,NP,NH,N,G,UP)
C*******************************************************************************
C PAIRING CORRECTION *
C ONE-FERMION FORMULA: C. KALBACH, NUCL.SCI.ENG. 95, 70 (1987) *
C*******************************************************************************
COMMON/UFU/DELTA2,ETH
AK=0.
ETH=0.
DELTA2=0.
P=FLOAT(NP)
H=FLOAT(NH)
PM=AMAX1(P,H)
IF(PM.EQ.0.) GOTO 10
FPH=0.
AKT=0.
IF(UP.GT.0.) THEN
ANC=1.584336*SQRT(G*UP)
RNC=FLOAT(NP+NH)/ANC
IF(RNC.GE.0.446) THEN
EPHA=UP*(0.716+2.44*RNC**2.17)
ELSE
EPHA=0.
ENDIF
IF(E1.GE.EPHA) THEN
R=AMAX1(0.,0.996-1.76*RNC**1.60/(E1/UP)**0.68)
ELSE
R=0.
ENDIF
DELTA2=4.*UP/G*R*R
ELSE
CONTINUE
ENDIF
cob ETH=(P*P+H*H)/(2.*G)
ETH=UP-G*DELTA2/4.+PM*SQRT((PM/G)**2+DELTA2)
IF(ETH.GT.E1) GO TO 10
AKS=(P*P+P+H*H+H)/(4.*G)
FPH=12.*(FLOAT(NP+NH)/FLOAT(N))
IF(PM.GT.0.0) FPH=FPH+(4.*G/PM)*(E1-ETH)
IF(FPH.GT.0.) AKT=(P*P+H*H-2.*(P+H)+2.)/(G*FPH)
cob AK=ETH-AKS
AK=ETH-AKS+AKT
10 RETURN
END
C
FUNCTION WK2(E1,NPP0,NHP0,NPN0,NHN0,GP,GN,UPP,UPN,S)
C*******************************************************************************
C TWO-FERMION FORMULA: C. KALBACH, NUCL.SCI.ENG. 95,70(1987) *
C*******************************************************************************
REAL*8 FCTR,T1,T2
COMMON/UFU/DELTA2,ETH
C
WK2=0.
NPP=NPP0-1
NHP=NHP0-1
NPN=NPN0-1
NHN=NHN0-1
NP=NPP+NHP
NN=NPN+NHN
N=NP+NN
PP=FLOAT(NPP)
HP=FLOAT(NHP)
PN=FLOAT(NPN)
HN=FLOAT(NHN)
PPM=AMAX1(PP,HP)
PNM=AMAX1(PN,HN)
EP=E1*FLOAT(NP)/FLOAT(N)
EN=E1*FLOAT(NN)/FLOAT(N)
C
IF(EP.GT.0.) THEN
AKP=AK(EP,NPP,NHP,N,GP,UPP)
IF(ETH.GE.EP) GOTO 10
ETHP=ETH
ELSE
AKP=0.
ETHP=0.
ENDIF
IF(EN.GT.0.) THEN
AKN=AK(EN,NPN,NHN,N,GN,UPN)
IF(ETH.GE.EN) GOTO 10
DN2=DELTA2
ETHN=ETH
ELSE
AKN=0.
DN2=0.
ETHN=0.
ENDIF
AK2=AKP+AKN
ETH=ETHP+ETHN
Ecor=AMAX1(0., E1-AK2-S )
T1=DFLOAT(GP**(PP+HP))/(FCTR(NPP)*FCTR(NHP))
T2=DFLOAT(GN**(PN+HN))/(FCTR(NPN)*FCTR(NHN))
WK2=SNGL( T1*T2/FCTR(N-1)*DFLOAT(Ecor**(N-1)) )
10 RETURN
END
C
FUNCTION WK3(E1,NP0,NH0,G,F,F10)
C*******************************************************************************
C ONE-FERMION FORMULA: C. KALBACH, PHYS.REV. C32,1157(1985) *
C*******************************************************************************
REAL*8 FCTR,T1,T2
C
WK3=0.
NP=NP0-1
NH=NH0-1
N=NP+NH
P=FLOAT(NP)
H=FLOAT(NH)
IF(H.LE.1.) THEN
F1=F10
ELSE
F1=F
ENDIF
C
Ep=0.
Eh=0.
IF(NP.EQ.0) Eh=E1/H
IF(NH.EQ.0) Ep=E1/P
IF(NP*NH.EQ.0) THEN
CONTINUE
ELSE
T11=FDC0(NP,NH,E1,F1)
IF(T11.LE.0.) GO TO 30
T12=FDC0(NP+1,NH,E1,F1)
RR=T12/T11
Ep=E1*RR/FLOAT(N)
Eh=(E1-P*Ep)/H
ENDIF
IF(Eh.GE.F) GO TO 30
GP=G*SQRT(1.+Ep/F)
GH=G*SQRT(1.-Eh/F)
GAV=((GP**P)*(GH**H))**(1./(P+H))
Alpha=(P*P+H*H)/(2.*GAV)
Aph=(P*(P-1.)+H*(H-1.))/(4.*GAV)
Ecor=E1-Aph
IF(Ecor.LE.0.) GO TO 30
C
T1=DFLOAT(GP**P)/FCTR(NP)
T2=DFLOAT(GH**H)/FCTR(NH)
WK3=SNGL(T1*T2/FCTR(N-1)*DFLOAT(Ecor**(N-1)*T11))
30 RETURN
END
C
FUNCTION FDC0(NP,NH,E1,F)
C*******************************************************************************
C FINITE-DEPTH CORRECTION FUNCTION f(p,h,E,F) *
C ONE-FERMION FORMULA: C. KALBACH, PHYS.REV. C32,1157(1985) *
C*******************************************************************************
REAL*8 FCTR
C
NN=NP+NH
SUM=0.
DO 20 JJ=1,NH+1
J=JJ-1
Ecorr=E1-J*F
IF(Ecorr.LE.0.) GO TO 20
D=(-1.)**J
CH=SNGL(FCTR(NH)/(FCTR(NH-J)*FCTR(J)))
SUM=SUM+D*CH*(Ecorr/E1)**(NN-1)
20 CONTINUE
FDC0=SUM
IF(FDC0.LT.0.) WRITE(6,*)' WARNING: FDC0.LT.0. '
RETURN
END
C
FUNCTION WK4(E1,NPP0,NHP0,NPN0,NHN0,GP,GN,FP,FN,F1P0,F1N0)
C*******************************************************************************
C TWO-FERMION FORMULA: C.KALBACH, PHYS.REV.C32,1157(1985); C33,818(1986) *
C*******************************************************************************
REAL*8 FCTR,T1P,T1N
WK4=0.
NPP=NPP0-1
NHP=NHP0-1
NPN=NPN0-1
NHN=NHN0-1
NP=NPP+NHP ! number of proton excitons
NN=NPN+NHN ! number of neutron excitons
N=NP+NN
PP=FLOAT(NPP) ! proton particle
HP=FLOAT(NHP) ! proton hole
PN=FLOAT(NPN) ! neutron particle
HN=FLOAT(NHN) ! neutron hole
IF(HP.LE.2.) THEN
F1P=F1P0
ELSE
F1P=FP
ENDIF
IF(HN.LE.2.) THEN
F1N=F1N0
ELSE
F1N=FN
ENDIF
EP=E1*FLOAT(NP)/FLOAT(N) ! proton-excitons excitation energy
EN=E1*FLOAT(NN)/FLOAT(N) ! neutron-excitons excitation energy
GPAV=0.
GNAV=0.
C
C protons
EPP=0.
EHP=0.
IF(NPP*NHP.EQ.0) THEN
IF(NPP.NE.0) EPP=EP/PP
IF(NHP.NE.0) EHP=EP/HP
ELSE
T11P=FDC0(NPP,NHP,EP,F1P)
IF(T11P.LE.0.) GO TO 1
T12P=FDC0(NPP+1,NHP,EP,F1P)
RP=T12P/T11P
EPP=EP*RP/FLOAT(NP)
EHP=(EP-PP*EPP)/HP
ENDIF
IF(EHP.GE.F1P) GO TO 30
1 CONTINUE
GPP=GP*SQRT(1.+EPP/FP)
GHP=GP*SQRT(1.-EHP/FP)
IF(NP.GT.0) GPAV=((GPP**PP)*(GHP**HP))**(1./NP)
C
C neutrons
EPN=0.
EHN=0.
IF(NPN*NHN.EQ.0) THEN
IF(NPN.NE.0) EPN=EN/PN
IF(NHN.NE.0) EHN=EN/HN
ELSE
T11N=FDC0(NPN,NHN,EN,F1N)
IF(T11N.LE.0.) GO TO 2
T12N=FDC0(NPN+1,NHN,EN,F1N)
RN=T12N/T11N
EPN=EN*RN/FLOAT(NN)
EHN=(EN-PN*EPN)/HN
ENDIF
IF(EHN.GE.F1N) GO TO 30
2 CONTINUE
GPN=GN*SQRT(1.+EPN/FN)
GHN=GN*SQRT(1.-EHN/FN)
IF(NN.GT.0) GNAV=((GPN**PN)*(GHN**HN))**(1./NN)
C
IF(GPAV.EQ.0.)THEN
AphP=0.
AlphaP=0.
ELSE
AphP=(PP*(PP-1.)+HP*(HP-1.))/(4.*GPAV) ! Pauli term for protons
AlphaP=(PP*PP+HP*HP)/(2.*GPAV) ! alpha term for protons
IF(AlphaP.GE.EP) GOTO 30
ENDIF
IF(GNAV.EQ.0.)THEN
AphN=0.
AlphaN=0.
ELSE
AphN=(PN*(PN-1.)+HN*(HN-1.))/(4.*GNAV) ! Pauli term for neutrons
AlphaN=(PN*PN+HN*HN)/(2.*GNAV) ! alpha term for neutrons
IF(AlphaN.GE.EN) GOTO 30
ENDIF
Aph=AphP+AphN
Alpha=AlphaP+AlphaN
Ecor=E1-Aph
C
T1P=DFLOAT(GPP**PP)/FCTR(NPP)*(DFLOAT(GHP**HP)/FCTR(NHP))
T1N=DFLOAT(GPN**PN)/FCTR(NPN)*(DFLOAT(GHN**HN)/FCTR(NHN))
WK4=SNGL(T1P*T1N/FCTR(N-1)*DFLOAT(Ecor**(N-1)))*
* FDC0(NPP+NPN,NHP+NHN,E1,(F1P+F1N)/2.)
30 RETURN
END
C
FUNCTION WM1(E1,NP0,NH0,G,F,BB,UP)
C*******************************************************************************
C ONE FERMION FORMULA: MAO MING DE+, J. PHYS. G 19,421(1993) *
C [ S. HILAIRE & A. KONNING (SMOLENICE '95) ] *
C*******************************************************************************
REAL*8 FCTR,CP,CH,T1,T2,SUM
DIMENSION B(0:31),C(0:31,0:31),BM(0:15,0:15,0:30)
COMMON/UFU/DELTA2,ETH
DATA B/ 1.,-.5,0.166666666,0.,-0.03333333,0.,0.023809523,
1 0.,-0.03333333,0.,0.075757576,0.,-0.253113553,0.,
2 1.166666667,0.,-7.092156863,0.,54.97117794,0.,-529.1242424,
3 0.,6192.123188,0.,-86580.25311,0.,1425517.167,0.,2729823.138,
4 0.,601580873.9,0./
C
C(0,0)=1.
DO 10 I=1,31
C(0,I)=0.
10 CONTINUE
DO 12 K=1,31
DO 12 LL=1,31
L=LL-1
CTOT=0.
DO 12 IKK=1,LL
IK=IKK-1
CTOT=CTOT+B(IK)*((-K/G)**IK)*C(K-1,L-IK)/FCTR(IK)
C(K,L)=CTOT
12 CONTINUE
C
DO 15 II1=1,16
I1=II1-1
DO 15 JJ1=1,16
J1=JJ1-1
N1=I1+J1
DO 15 LL=1,N1+1
L=LL-1
BM(I1,J1,L)=0.
DO 14 II2=1,L+1
I2=II2-1
DO 14 JJ2=1,L+1
J2=JJ2-1
IF(I2+J2.NE.L)GO TO 14
BM(I1,J1,L)=BM(I1,J1,L)+C(I1,I2)*C(J1,J2)
14 CONTINUE
15 CONTINUE
C
WM1=0.
NP=NP0-1
NH=NH0-1
NN=NP+NH
P=FLOAT(NP)
H=FLOAT(NH)
PM=AMAX1(P,H)
T1=DFLOAT(G**(P+H))/(FCTR(NP)*FCTR(NH))
AK1=AK(E1,NP,NH,NN,G,UP)
IF(E1+ETH-P*BB-H*F.GT.0.)GO TO 30
C
SUM=0.
DO 25 II=1,NP+1
I=II-1
DO 25 JJ=1,NH+1
J=JJ-1
D=(-1.)**(I+J)
CP=FCTR(NP)/(FCTR(NP-I)*FCTR(I))
CH=FCTR(NH)/(FCTR(NH-J)*FCTR(J))
Ecorr=E1-ETH-I*BB-J*F
IF(Ecorr.LE.0.)GO TO 25
DO 20 ILL=1,NN
IL=ILL-1
T2=CP*CH/FCTR(NN-1-IL)
SUM=SUM+D*T2*(Ecorr**(NN-1-IL))*BM(NP,NH,IL)
20 CONTINUE
25 CONTINUE
WM1=T1*DMAX1(0.,SUM)
IF(WM1.LT.0.) WRITE(6,*) ' WARNING: WM1.LT.0. ! '
30 RETURN
END
C
FUNCTION WM2(E1,NPP0,NHP0,NPN0,NHN0,GP,GN,FP,FN,BP,BN,UPP,UPN,S)
C*******************************************************************************
C TWO FERMION FORMULA: MAO MING DE, J. PHYS. G 19,421(1993), *
C A. KONNING (SMOLENICE '95) *
C*******************************************************************************
REAL*8 FCTR,T1,T2,SUM
DIMENSION A(0:30),B(0:31),C(0:31,0:31),BM(0:15,0:15,0:30)
COMMON/UFU/DELTA2,ETH
DATA B/ 1.,-.5,0.166666666,0.,-0.03333333,0.,0.023809523,
1 0.,-0.03333333,0.,0.075757576,0.,-0.253113553,0.,
2 1.166666667,0.,-7.092156863,0.,54.97117794,0.,-529.1242424,
3 0.,6192.123188,0.,-86580.25311,0.,1425517.167,0.,2729823.138,
4 0.,601580873.9,0./
C
WM2=0.
NPP=NPP0-1
NHP=NHP0-1
NPN=NPN0-1
NHN=NHN0-1
NP=NPP+NHP ! number of protron excitons
NN=NPN+NHN ! number of neutron excitons
N=NN+NP ! total number of excitons
IF(MAX1(NPP,NHP,NPN,NHN).GT.15.OR.N.GT.31) THEN
WRITE(6,*)' STOP: Vector-size limit for conf.',npp,nhp,npn,nhn
GOTO 40
ELSE
CONTINUE
ENDIF
PP=FLOAT(NPP)
HP=FLOAT(NHP)
PN=FLOAT(NPN)
HN=FLOAT(NHN)
GG=(GP**NP*GN**NN)**(1./FLOAT(N))
EP=E1*FLOAT(NP)/FLOAT(N)
EN=E1*FLOAT(NN)/FLOAT(N)
C
C(0,0)=1
DO 5 I=1,31
C(0,I)=0.
5 CONTINUE
DO 12 K=1,31
DO 12 LL=1,31
L=LL-1
C(K,L)=0.
DO 10 IKK=1,LL
IK=IKK-1
C(K,L)=C(K,L)+B(IK)*((-K/GG)**IK)*C(K-1,L-IK)/SNGL(FCTR(IK)) ! Eq.(23)
10 CONTINUE
12 CONTINUE
C
DO 15 II1=1,16
I1=II1-1
DO 15 JJ1=1,16
J1=JJ1-1
N1=I1+J1
DO 15 LL=1,N1+1
L=LL-1
BM(I1,J1,L)=0.
DO 14 II2=1,L+1
I2=II2-1
DO 14 JJ2=1,L+1
J2=JJ2-1
ITEST=I2+J2
IF(ITEST.NE.L)GO TO 14
BM(I1,J1,L)=BM(I1,J1,L)+C(I1,I2)*C(J1,J2) ! Eq.(22)
14 CONTINUE
15 CONTINUE
C
DO 20 IM=1,N+1
IL=IM-1
A(IL)=0.
DO 20 IN=1,IL+1
I=IN-1
20 A(IL)=A(IL)+BM(NPN,NHN,I)*BM(NPP,NHP,IL-I) ! Eq.(33)
C
IF(EP.GT.0.) THEN
AKP=AK(EP,NPP,NHP,N,GP,UPP)
IF(ETH.GE.EP) GOTO 40
ETHP=ETH
ELSE
ETHP=0.
ENDIF
IF(EN.GT.0.) THEN
AKN=AK(EN,NPN,NHN,N,GN,UPN)
IF(ETH.GE.EN) GOTO 40
ETHN=ETH
ELSE
ETHN=0.
ENDIF
ETH=ETHP+ETHN
C
IF((E1+ETH-PP*BP-HP*FP-PN*BN-HN*FN).GT.0.) GO TO 40
T1=SNGL((GP**NP)/(FCTR(NPP)*FCTR(NHP)))
T2=SNGL((GN**NN)/(FCTR(NPN)*FCTR(NHN)))
SUM=0.
DO 30 IIP=1,NPP+1
IP=IIP-1
DO 30 JJP=1,NHP+1
JP=JJP-1
DO 30 IIN=1,NPN+1
IN=IIN-1
DO 30 JJN=1,NHN+1
JN=JJN-1
DO 30 ILL=1,N
IL=ILL-1
Ecor=E1-ETH-IP*BP-IN*BN-JP*FP-JN*FN-S ! Eq.(36)
IF(Ecor.LE.0.)GO TO 30
D=(-1)**(IP+IN+JP+JN)
CPP=FCTR(NPP)/(FCTR(NPP-IP)*FCTR(IP))
CHP=FCTR(NHP)/(FCTR(NHP-JP)*FCTR(JP))
CPN=FCTR(NPN)/(FCTR(NPN-IN)*FCTR(IN))
CHN=FCTR(NHN)/(FCTR(NHN-JN)*FCTR(JN))
SUM=SUM+D*CPP*CHP*CPN*CHN/FCTR(N-1-IL)*(Ecor**(N-1-IL))*A(IL)
30 CONTINUE
WM2=T1*T2*DMAX1(0.,SUM) ! Eq.(35)
40 RETURN
END
C
FUNCTION WR1(E1,NP0,NH0,G,F,F10,B,UP,S)
C*******************************************************************************
C COMPOSITE ONE-FERMION FORMULA *
C*******************************************************************************
REAL*8 FCTR,T1,T2
COMMON/UFU/DELTA2,ETH
C
WR1=0.
NP=NP0-1
NH=NH0-1
N=NP+NH
P=FLOAT(NP)
H=FLOAT(NH)
IF(H.LE.2.0.AND.F10.GT.0.) THEN
F1=F10
ELSE
F1=F
ENDIF
C
Ep=0.
Eh=0.
IF(F10.LT.0.) GOTO 15
IF(NP*NH.EQ.0) THEN
IF(NP.NE.0) Ep=E1/P
IF(NH.NE.0) Eh=E1/H
ELSE
T11=FDC(E1,NP,NH,0,G,F1,B,UP,S)
IF(T11.LE.0.) GO TO 30
T12=FDC(E1,NP,NH,1,G,F1,B,UP,S)
RR=T12/T11
Ep=E1*RR/FLOAT(N)
Eh=(E1-P*Ep)/H
ENDIF
IF(Eh.GE.F) GO TO 30
15 GP=G*SQRT(1.+Ep/F)
GH=G*SQRT(1.-Eh/F)
GAV=((GP**P)*(GH**H))**(1./(P+H))
C
AK1=AK(E1,NP,NH,N,GAV,UP)
IF(E1+ETH-P*B-H*F.GT.0.)GO TO 30
T1=DFLOAT(GP**P)/FCTR(NP)
T2=DFLOAT(GH**H)/FCTR(NH)
WR1=SNGL(T1*T2/FCTR(N-1)*DFLOAT(E1**(N-1)))*
* FDC(E1,NP,NH,0,GAV,F1,B,UP,S)
30 RETURN
END
C
FUNCTION FDC(E1,NP,NH,IK,G,F,B,UP,S)
C*******************************************************************************
C EXTENDED FINITE-DEPTH CORRECTION FUNCTION f(p,h,E,F,AK,B) *
C BASED ON: C.KALBACH, PHYS.REV.C32,1157(1985), NUCL.SCI.ENG.95,70 (1987)*
C*******************************************************************************
REAL*8 FCTR,SUM
COMMON/UFU/DELTA2,ETH
C
N=NP+NH
FK1=1.
AK1=AK(E1,NP,NH,N,G,UP)
SUM=0.
DO 20 II=1,NP+1
I=II-1
DO 20 JJ=1,NH+1
J=JJ-1
Ecor=E1-AK1-I*B-J*F-S
Ecor1=E1-ETH-I*B-J*F-S
IF(Ecor1.LE.0.) GO TO 20
D=(-1.)**(I+J)
CP=SNGL(FCTR(NP)/(FCTR(NP-I)*FCTR(I)))
CH=SNGL(FCTR(NH)/(FCTR(NH-J)*FCTR(J)))
N1=N-1+IK
IF(IK.EQ.1) FK1=(1.+FLOAT(N)/FLOAT(NP)*(I*B/Ecor))**IK
SUM=SUM+D*CP*CH*(Ecor/E1)**N1*FK1
20 CONTINUE
FDC=SUM
IF(FDC.LT.0.) FDC=0.
RETURN
END
C
FUNCTION WR2(E1,NPP0,NHP0,NPN0,NHN0,GP,GN,FP,FN,F1P0,F1N0,B,BN,
* UPP,UPN,S)
C*******************************************************************************
C COMPOSITE TWO-FERMION FORMULA *
C*******************************************************************************
REAL*8 FCTR,T1P,T1N,sum
COMMON/UFU/DELTA2,ETH
WR2=0.
NPP=NPP0-1
NHP=NHP0-1
NPN=NPN0-1
NHN=NHN0-1
NP=NPP+NHP ! number of proton excitons
NN=NPN+NHN ! number of neutron excitons
N=NP+NN
PP=FLOAT(NPP) ! proton particle
HP=FLOAT(NHP) ! proton hole
PN=FLOAT(NPN) ! neutron particle
HN=FLOAT(NHN) ! neutron hole
IF(HP.LE.2.0.AND.F1P0.GT.0.) THEN
F1P=F1P0
ELSE
F1P=FP
ENDIF
IF(HN.LE.2.0.AND.F1N0.GT.0.) THEN
F1N=F1N0
ELSE
F1N=FN
ENDIF
F1=(F1P+F1N)/2.
EP=E1*FLOAT(NP)/FLOAT(N) ! proton-excitons excitation energy
EN=E1*FLOAT(NN)/FLOAT(N) ! neutron-excitons excitation energy
C
C protons
EPP=0.
EHP=0.
IF(F1P0.LT.0.) GOTO 1
IF(NPP*NHP.EQ.0) THEN
IF(NPP.NE.0) EPP=EP/PP
IF(NHP.NE.0) EHP=EP/HP
ELSE
T11P=FDC(EP,NPP,NHP,0,GP,F1P,B,UPP,S)
IF(T11P.LE.0.) GO TO 1
T12P=FDC(EP,NPP,NHP,1,GP,F1P,B,UPP,S)
RP=T12P/T11P
EPP=EP*RP/FLOAT(NP)
EHP=(EP-PP*EPP)/HP
ENDIF
IF(EHP.GE.FP) GO TO 40
1 CONTINUE
GPP=GP*SQRT(1.+EPP/FP)
GHP=GP*SQRT(1.-EHP/FP)
IF(NP.GT.0) GPAV=((GPP**PP)*(GHP**HP))**(1./NP)
C
C neutrons
EPN=0.
EHN=0.
IF(F1N0.LT.0.) GOTO 2
IF(NPN*NHN.EQ.0) THEN
IF(NPN.NE.0) EPN=EN/PN
IF(NHN.NE.0) EHN=EN/HN
CONTINUE
ELSE
T11N=FDC(EN,NPN,NHN,0,GN,F1N,BN,UPN,S)
IF(T11N.LE.0.) GO TO 2
T12N=FDC(EN,NPN,NHN,1,GN,F1N,BN,UPN,S)
RN=T12N/T11N
EPN=EN*RN/FLOAT(NN)
EHN=(EN-PN*EPN)/HN
ENDIF
IF(EHN.GE.FN) GO TO 40
2 CONTINUE
GPN=GN*SQRT(1.+EPN/FN)
GHN=GN*SQRT(1.-EHN/FN)
IF(NN.GT.0) GNAV=((GPN**PN)*(GHN**HN))**(1./NN)
GAV=(GPP**PP*GHP**HP*GPN**PN*GHN**HN)**(1./N)
C
IF(EP.GT.0.) THEN
csepf FDCP=FDC(EP,NPP,NHP,0,GPAV,F1P,B,UPP,S)
AKP=AK(EP,NPP,NHP,N,GPAV,UPP)
IF(ETH.GE.EP) GOTO 40
ETHP=ETH
ELSE
csepf FDCP=1.
AKP=0.
ETHP=0.
ENDIF
IF(EN.GT.0.) THEN
csepf FDCN=FDC(EN,NPN,NHN,0,GNAV,F1N,BN,UPN,S)
AKN=AK(EN,NPN,NHN,N,GNAV,UPN)
IF(ETH.GE.EN) GOTO 40
ETHN=ETH
ELSE
csepf FDCN=1.
AKN=0.
ETHN=0.
ENDIF
AK2=AKP+AKN
ETH=ETHP+ETHN
C
IF((E1+ETH-PP*B-HP*F1P-PN*BN-HN*F1N).GT.0.) GO TO 40
T1P=DFLOAT(GPP**PP)/FCTR(NPP)*(DFLOAT(GHP**HP)/FCTR(NHP))
T1N=DFLOAT(GPN**PN)/FCTR(NPN)*(DFLOAT(GHN**HN)/FCTR(NHN))
csepf WR2=SNGL(T1P*T1N/FCTR(N-1)*DFLOAT(E1**(N-1)))*FDCP*FDCN
csepf GOTO 35
C
SUM=0.
DO 30 IIP=1,NPP+1
IP=IIP-1
DO 30 JJP=1,NHP+1
JP=JJP-1
DO 30 IIN=1,NPN+1
IN=IIN-1
DO 30 JJN=1,NHN+1
JN=JJN-1
Ecor=E1-AK2-IP*B-IN*BN-JP*F1P-JN*F1N-S
IF(Ecor.LE.0.)GO TO 30
D=(-1)**(IP+IN+JP+JN)
CPP=FCTR(NPP)/(FCTR(NPP-IP)*FCTR(IP))
CHP=FCTR(NHP)/(FCTR(NHP-JP)*FCTR(JP))
CPN=FCTR(NPN)/(FCTR(NPN-IN)*FCTR(IN))
CHN=FCTR(NHN)/(FCTR(NHN-JN)*FCTR(JN))
SUM=SUM+D*CPP*CHP*CPN*CHN/FCTR(N-1)*(Ecor**(N-1))
30 CONTINUE
WR2=SNGL(T1P*T1N*SUM)
35 IF(WR2.LT.0.) WRITE(6,*) ' WARNING: WR2.LT.0. ! '
40 RETURN
END