program duflo28p common/cpair/epair open (unit=10,file='mass_duflo.out') do 10 kz=20,100 zz=kz aa=zz+zz do 20 m=1,5 20 aa=zz*(1.983588+0.016412*aa**0.6666667) ia=aa-2.5 kn=ia-kz kn=kn+mod(kn,2) ka=kn+kz c read(*,*) kz,ka c if (kz.eq.0) goto 1000 xz=float(kz) xn=float(ka-kz) call emassdz(xz,xn,y) emasse=8.07138*xn+7.28903*xz-y write(10,2) kz,kn,ka,y,epair write(*,2) kz,kn,ka,y,epair c 2 format('kz=',i2,' ka=',i3,' M=',f9.2) 2 format(3i4,2f9.2) 10 continue end c SUBROUTINE EMASSDZ (Z,X,Y) ! Z protons, X neutrons , Y bind.Energy common/cpair/epair DIMENSIONA(33),DYDA(33),FYDA(33),op(2,3,2),fyd0(33),onps(2),n4(2), 1ONP(0:8,2,2),OT(0:8,2,2),OEI(2),DEI(2),NN(2),NOC(18,2),OP2(2), 2JUP(2),JUD(2),YM(2),OP1(2),n2(2),shell(2),sshell(2) CC***************************************** data a/ 9.0914, 6.3355, 4.5791, 19.8946, 1.7325, 7.5247, & -7.1953, -39.9787, -0.3976, 0.8131, -0.7435, -3.7291, & -0.1305, -0.6387, 0.4534, 2.0605, 0.3449, 1.4727, & -1.0433, 0.0000, 5.2495, 0.0000, -32.1007,-151.1164, & -4.6103, -32.4238, -37.3226, -52.1673, 0.9597, 3.0024, & 0.6977, 6.0390, 17.7960/ c Data=1751 RMS= 0.330 (mars 95) cFM+* 9.09 fm+* 6.34 FS+* 4.58 fs+* 19.89 FS-* 1.73 fs-* 7.52 cFC+* -7.20 fc+* -39.98 PM+* -0.40 pm+* 0.81 PS+* -0.74 ps+* -3.73 cPS-* -0.13 ps-* -0.64 S3 * 0.45 s3 * 2.06 SQ-* 0.34 sq-* 1.47 cD3 * -1.04 d3 0.00 QQ+* 5.25 qq+ 0.00 D0* -32.10 d0*-151.12 cQQ-* -4.61 qq-* -32.42 TT* -37.32 tt* -52.17 SS* 0.96 ss* 3.00 c C * 0.70 P0 * 6.04 P1 * 17.80 c ---------------------- CC DEBUT DU SS PROG. c c a(27)=0. c a(28)=0. c a(29)=0. c a(30)=0. c a(31)=0. c a(32)=0. c a(33)=0. IMAX=18 MO=2 MAXP=8 NN(1)=X NN(2)=Z nx=x nz=z V=X+Z T=ABS(X-Z) r=v**(1./3.) s=r*r Rc=R*(1.-.25*(t/v)**2) Ra=(rc*rc)/r C****************** DO NDEF=1,2 ! NDEF=1 sph. NDEF=2 def. ym(ndef)=0. C NDEF= 1->SPH// 2->DEF// JU=0 JUP(1)=0 JUP(2)=0 JUD(1)=0 JUD(2)=0 C----------- IF(NDEF.EQ.2.and.nz.gt.50)JU=4 DO KK=1,33 FYDA(KK)=0. FYD0(KK)=0. dyda(KK)=0. END DO DO J=1,2 DO I=1,IMAX NOC(I,J)=0 ENDDO DO K=0,MAXP DO M=1,MO ONP(K,M,J)=0. ENDDO ENDDO ENDDO DO J=1,2 !LOOP OVER FLUIDS (1=N,2=P) NCUM=0 I=0 20 I=I+1 i2=(I/2)*2 if(i2.ne.i)then ID=I+1 if(ncum.lt.nn(j))sshell(J)=1. !sscouche J else ID=I*(I-2)/4 if(ncum.lt.nn(j))sshell(J)=2. ! SSC R endif NCUM=NCUM+ID IF(NCUM.LE.NN(J))THEN NOC(I,J)=ID GO TO 20 ENDIF shell(J)=I !N0 de sscouche (SSc 2-> 0 nucl) IP=(I-1)/2 !N0 couche HO MOC=NN(J)-NCUM+ID IF(NDEF.EQ.2) THEN if(i2.ne.i)then JUD(J)=MAX (0,JU-MOC) JUP(J)=0 Else JUP(J)=MIN (JU,MOC) JUP(J)=JU JUD(J)=0 ENDIF ENDIF NOC(I,J)=MOC-JUP(J)+JUD(J) NOC(I+1,J)=JUP(j) NOC(I-1,J)=NOC(I-1,J)-JUD(J) if(i2.ne.i)then OEI(J)=MOC+IP*(IP-1)-JU DEI(J)=IP*(IP+1)+2 ELSE OEI(J)=MOC-JU DEI(J)=(IP+1)*(IP+2)+2 ENDIF CC HERE,DEGENERACIES AND NUMBER OF ACTIVE PARTICLES FOR EI. IPL=0 vmr=0. vmj=0. DO II=1,IMAX onps(j)=0. IP=(II-1)/2 DEGI=(IP+1)*(IP+2) FAC=1./SQRT(DEGI) IF(IP.NE.IPL)IPL=IPL+1 if((2*IP+1).eq.II)then VM2=(.5*IP)/(IP+1) degr=ip*(ip-1) if(ip.gt.2)then vmr=(.5*(ip-1))/ip vmj=-1./ip if(noc(ii,j).le.degr)onps(j)=noc(ii,j)*vmr if(noc(ii,j).gt.degr)onps(j)=degr*vmr+(noc(ii,j)-degr)*vmj endif endif if((2*IP+1).ne.II)VM2=-1./(IP+1) !SSc. j ONP(IPL,2,J)=ONP(IPL,2,J)+NOC(II,J)*VM2 ONP(IPL,1,J)=ONP(IPL,1,J)+NOC(II,J)*FAC fyd0(29)=fyd0(29)+onps(j)*(onp(ipl,1,j)+onp(ipl,2,j)) ENDDO ENDDO !END OF LOOP OVER FLUIDS c************** IF(NDEF.EQ.2) THEN ALFA=0. ELSE ALFA=0.5 ENDIF FACN=DEI(1)**ALFA FACZ=DEI(2)**ALFA DNNB=OEI(1)*(DEI(1)-OEI(1))/DEI(1) DZZB=OEI(2)*(DEI(2)-OEI(2))/DEI(2) QN=DNNB*FACN/SQRT(DEI(1)) QZ=DZZB*FACZ/SQRT(DEI(2)) DI1N=DNNB*(2*OEI(1)-DEI(1))*FACN*FACN/(DEI(1)) DI1Z=DZZB*(2*OEI(2)-DEI(2))*FACZ*FACZ/(DEI(2)) S3=DI1Z+DI1N QQ0=(QN+QZ)*(QN+QZ) QQ1=(QN-QZ)*(QN-QZ) qqp=qq0+qq1 qqm=qq0-qq1 c------------ DO M=1,MO DO I=0,MAXP OT(I,M,1)=ONP(I,M,1)+ONP(I,M,2) OT(I,M,2)=ONP(I,M,1)-ONP(I,M,2) ENDDO ENDDO c------------ DO L=1,2 DO M=1,3 DO J=1,2 OP(L,M,J)=0. ENDDO ENDDO ENDDO DO I=0,MAXP DEGI=(I+1)*(I+2) FAC=SQRT(DEGI) DO N=1,2 DO M=1,MO OP(1,M,N)=OP(1,M,N)+OT(I,M,N) OTX=OT(I,M,N)*OT(I,M,N)*FAC OP(2,M,N)=OP(2,M,N)+OTX ENDDO ENDDO ENDDO DO N=1,2 OP1(N)=0. OP2(N)=0. DO M=1,MO OPXX=OP(1,M,N)*OP(1,M,N) OP(1,M,N)=OPXX ENDDO ENDDO DO N=1,2 DO I=0,MAXP DEGI=(I+1)*(I+2) FAC=sqrt(degi) OP1(N)=OP1(N)+OT(I,1,N)/FAC OP2(N)=OP2(N)+OT(I,2,N)/FAC ENDDO ENDDO DO N=1,2 OP(1,3,N)=OP1(N)*OP2(N) ENDDO K=-1 DO L=1,2 DO M=1,3 DO N=1,2 K=K+2 FYD0(K)=OP(L,M,N) ENDDO ENDDO ENDDO C----------- DO JF=1,17,4 FYD0(JF)=FYD0(JF)+FYD0(JF+2) FYD0(JF+2)=FYD0(JF)-2.*FYD0(JF+2) ENDDO C----------- fyda(1)=fyd0(1) fyda(3)=fyd0(5) fyda(5)=fyd0(7) fyda(7)=fyd0(9) fyda(9)=fyd0(13) fyda(11)=fyd0(17) fyda(13)=fyd0(19) FYDA(27)= T*(T+2)/s fyda(29)=fyd0(29) C----------- IF(NDEF.EQ.1)then FYDA(15)=S3 FYDA(17)=QQm ELSE FYDA(19)=S3 FYDA(21)=QQp FYDA(23)=16.-qqm FYDA(25)=QQm ENDIF C----------- DO MSS=1,29,2 DYDA(MSS)=FYDA(MSS)/Ra DYDA(MSS+1)=-DYDA(MSS)/Ra ym(ndef)=ym(ndef)+dyda(mss)*a(mss)+dyda(mss+1)*a(mss+1) ENDDO c--------- Z2=Z*(Z-1) DYDA(31)=(-Z2+.76*Z2**(2./3.))/Rc !Coulomb c--------- rxz=1./Ra !Pairing vxz=1./v txz=rxz*(t/v) uxz=rxz-txz do l=1,2 n2(l)=2*(nn(l)/2) dyda(32)=-rxz+dyda(32) if(sshell(l).eq.2.) dyda(33)=vxz+dyda(33) !effet de couche en 1/A if(n2(l).eq.nn(l))dyda(32)=uxz+dyda(32) enddo j=2 ! Z>N if(nn(1).ge.nn(2))j=1 ! N>ou=Z k=3-j if(n2(j).eq.nn(j).and.n2(k).ne.nn(k))dyda(32)=-txz+dyda(32) c******** DO MSS=31,33 ym(ndef)=ym(ndef)+dyda(mss)*a(mss) ENDDO ENDDO !End sph and def calculations Y=ym(1) if(z.gt.50.and.ym(2).gt.ym(1)) Y=ym(2) av=1.717 v4=v/(1-(t/2./v)**2.)**2*(1.+2./9.*(t/v)**2) c v4=v c Y=Y-v4*(a(1)+a(9))*av+v4/Ra*(a(2)+a(10))*av c Y=Y-dyda(1)*a(1)+dyda(2)*a(2)-dyda(9)*a(9)+dyda(10)*a(10) t22=t*(t+2) Y=Y-15.07*v**(4./3.)/Ra+15.00*v**(4./3.)/Ra**2 & +33.88*t22/v**(2./3.)/Ra-48.51*t22/v**(2./3.)/Ra**2 epair=dyda(32)*a(32)+dyda(33)*a(33) END