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