PROGRAM HFSKY COMMON /DER/ FCT(100),DF(100),H,NNN COMMON RM(100),DU(100),VC(100),U(100),VS(100), + TAUN(100),TAUP(100),TAUT(100),DJN1(100),DJP1(100), +DJT1(100),VN(100),VP(100),VSON(100),VSOP(100), + ECOULB(100),VQRE(100),RM2(100),RM3(100),HMEN(100),HMEP(100), +DHMEN(100),DHMEP(100),D2HMEN(100),D2HMEP(100),DN(100),DP(100), +DT(100),DN1(100),DP1(100),DT1(100),DN2(100),DP2(100),DT2(100), +DJN(100),DJP(100),DJT(100),DC(100) + ,AVP(100 ,45),DUNL(100 ,45),UNL(100 ,45) DIMENSION NN(45),LL(45),LJ(45),LT(45),DEG(45),EHF(45) + ,EVSO(45) DIMENSION REC(8) , UNLOSC(100,45) DIMENSION NX(29),LX(29),JX(29) CHARACTER*6 SX(29),ST(45) REAL MA,MESMP(100),MESMN(100) DATA NNN,NEOC/100,45/ DATA CP/.79577471E-01/,QP/12.566371/ DATA SX/' 1S1/2',' 1P3/2',' 1P1/2',' 1D5/2',' 2S1/2',' 1D3/2', 1' 1F7/2',' 2P3/2',' 1F5/2',' 2P1/2',' 1G9/2',' 2D5/2',' 1G7/2', 2' 3S1/2',' 2D3/2','1H11/2',' 2F7/2',' 1H9/2','1I13/2',' 3P3/2', 3' 2F5/2',' 3P1/2','1I11/2','2G 9/2','1J15/2',' 3D5/2',' 4S1/2', 4' 2G7/2',' 3D3/2'/ DATA NX/1,1,1,1,2,1,1,2,1,2,1,2,1,3,2,1,2,1,1,3,2,3,1,2,1,3,4,2,3/ DATA LX/0,1,1,2,0,2,3,1,3,1,4,2,4,0,2,5,3,5,6,1,3,1,6,4,7,2,0,4,2/ DATA JX/1,3,1,5,1,3,7,3,5,1,9,5,7,1,3,11,7,9,13,3,5,1,11,9, + 15,5,1,7,3/ 200 FORMAT( //12H ITERATION I2,10H XMU=F5.3//) 201 FORMAT (2X,A6,I2,5X2HE=E11.4,5X4HDEG=E11.4,5X3H D=E8.2) 203 FORMAT(/) 204 FORMAT(5X,A6,I2,5X7HNON LIE) 205 FORMAT(/18H NOYAU DE MASSE A=F5.1,3X15HET DE CHARGE Z=I2/) 206 FORMAT(1H1) 209 FORMAT(8E10.3) 211 FORMAT(10E12.5) 212 FORMAT(10E12.5) 213 FORMAT(//18H DENSITE NEUTRONS/) 214 FORMAT(//18H DENSITE PROTONS /) 215 FORMAT(5E15.8) 216 FORMAT(//18H DENSITE CHARGE /) 217 FORMAT(/6H EC/A=E15.8,7H EHF/A=E15.8,8H EREA/A=E15.8,8H ETOT/A=E 1 15.8//) 226 FORMAT(/10H FICHIER I4//) 221 FORMAT(//4H RN=E13.5,4H RP=E13.5,4H RC=E13.5//) 400 FORMAT(20I3) 401 FORMAT(8F10.5) 403 FORMAT(8(F12.2,3X)) C CCC PARAMETRES DU CALCUL DE HARTREE-FOCK CCC C LECTURE DES PARAMETRES CCC C READ 401,T0,T1,T2,T3,X0,X1,X2,X3,GAM,W2 ACCEPT *,T0,T1,T2,T3,X0,X1,X2,X3,GAM,W2 PRINT 211,T0,T1,T2,T3,X0,X1,X2,X3,GAM,W2 W1=W2 *2. C 1 READ 400,NOS,NPO,NI,IOC1,IDEG1,IOC2,IDEG2,IOC3,IDEG3 C + ,IOC4,IDEG4,IOC5,IDEG5 1 ACCEPT *,NOS,NPO,NI,IOC1,IDEG1,IOC2,IDEG2,IOC3,IDEG3 + ,IOC4,IDEG4,IOC5,IDEG5 IF(NOS.EQ.0) GO TO 1000 PRINT 400,NOS,NPO,NI,IOC1,IDEG1,IOC2,IDEG2,IOC3,IDEG3 + ,IOC4,IDEG4,IOC5,IDEG5 C READ 400,KOP1,IDCC,IPUNCH ACCEPT *,KOP1,IDCC,IPUNCH PRINT 400,KOP1,IDCC,IPUNCH DO 2 J=1,NNN 2 U(J)=0. RMAX = 20. H = RMAX / NNN IPRE = 2 DO 95 J=1,NNN RM(J) = H*FLOAT(J) RM2(J) = RM(J) * RM(J) RM3(J) = RM(J)*RM2(J) 95 CONTINUE UNS4PI = CP NN1 = NNN- 1 NN2 = NNN- 2 NN3 = NNN- 3 NN4 = NNN - 4 W3 = W2 W0 = W1 CRO = .25*(T1*(1.+.5*X1)+T2*(1.+.5*X2)) CROQ = .125*(T2*(1.+2.*X2) -T1*(1.+2.*X1)) TXRO = T0*(1. + .5*X0) TXROQ = -T0*(X0 +.5) TXDRO = .5*(T2*(1.+.5*X2)-T1*(1.+.5*X1)) TXD2RO = .125*(T2*(1.+.5*X2) - 3.*T1*(1.+.5*X1)) TXDROQ = .25*(T1*(1.+2.*X1) + T2*(1.+2.*X2)) TXD2RQ = .0625*(3.*T1*(1. +2.*X1) + T2*(1.+2.*X2)) TXTAU = CRO TXTAUQ = CROQ DO 344 I = 1,NOS IF (I-NPO) 376,376,377 376 LT(I) = 1 J = I GO TO 378 377 LT(I) = 0 J = I - NPO 378 NN(I) = NX(J) LL(I) = LX(J) LJ(I) = JX(J) 344 ST(I) = SX(J) CCC CCC DEFINITIN DES FACTEURS D OCCUPATION DO 402 I = 1,NOS 402 DEG(I) = LJ(I) + 1 IF(IOC1.NE.0) DEG(IOC1)=FLOAT(IDEG1) IF(IOC2.NE.0) DEG(IOC2)=FLOAT(IDEG2) IF(IOC3.NE.0) DEG(IOC3)=FLOAT(IDEG3) IF(IOC4.NE.0) DEG(IOC4)=FLOAT(IDEG4) IF(IOC5.NE.0) DEG(IOC5)=FLOAT(IDEG5) CCC JZ = 0. MA = 0. DO 100 I=1,NOS JZ = JZ + LT(I)*DEG(I) 100 MA = MA + DEG(I) PRINT 205,MA,JZ DMSHB = .04823*MA/(MA-1) C HB = 1./DMSHB CCC CCC PUITS DE SAXON WOODS CCC R0 = 1.21 A1 = 55.5 A2 = 33.2 A3 = .36 A5 = .68 CWPI = 1.41 CW2 = CWPI*CWPI AL = A5 R = R0*MA**(1./3.) CCC CCC STARTING VALUE CCC ITER = 0 C VAX GO TO(137,341,601), KOP1 IF (KOP1.EQ.1) GOTO 137 IF (KOP1.EQ.2) GOTO 341 IF (KOP1.EQ.3) GOTO 601 TYPE *,' ====> VALORE NON AMMESSO DI KOP1 : ',KOP1 STOP 137 READ 209, ((AVP(J,I),J=1,NNN),I=1,NOS),(EHF(I),I=1,NOS) GO TO 107 341 CONTINUE DO 94 J=1,NNN HMEN(J) = HB HMEP(J) = HB DHMEN(J) = 0. DHMEP(J) = 0. D2HMEN (J) = 0. D2HMEP (J) = 0. MESMN(J) = 1. MESMP(J) = 1. 94 CONTINUE KX = 1 E = 0. E22 = 1.43986* H F = 0. RC = R*1.09/R0 SYM = (MA-2*JZ)/MA Z = 0. DO 96 K = 1,NNN X = RM(K) DP(K)=1./(1.+EXP((X-RC)/.55)) 96 Z = Z + X*X*DP(K) Y = FLOAT(JZ)/(QP*Z*H) DO 97 K = 1,NNN 97 DP(K) = Y*DP(K) DO 98 K=1,NNN Y = QP*RM(K)*RM(K) E = E + Y*DP(K) F = F + Y*DP(K)/RM(K) 98 VC(K) = E/RM(K) - F FF = F DO 101 J=1,NNN VC(J) = E22*(VC(J)+FF) X = RM(J) PE=EXP((X-R)/AL) FPE = -1./(1.+PE) VN(J) = FPE VS(J) = - PE*FPE*FPE/AL DO 101 I=1,NOS N = NN(I) L = LL(I) JM = LJ(I) FL = (JM*(JM+2) - 4*L*(L+1) - 3)/8. ISIG = 2*LT(I) - 1 V = A1 + ISIG*A2*SYM AVP(J,I) = V*VN(J) - A3*CW2*VS(J)*2.*FL IF (I.LE.NPO) AVP(J,I) = AVP(J,I) + VC(J) 101 CONTINUE DO 102 I=1,NOS EHF(I) = .3 * AVP(1,I) 102 IF(I.LE.NPO) EHF(I) = EHF(I) +.3*VC(1) 601 CONTINUE 107 CONTINUE ITMAX = NI CCC CCC DEBUT D ITERETION XMU=0.3 CCC 104 ITER = ITER + 1 IF(ITER.EQ.10) XMU=0.5 IF(ITER.EQ.15) XMU=0.7 IF(ITER.EQ.ITMAX) PRINT 200,ITER,XMU CCC CCC RESOLUTION DES EQUATIONS DE SCHORODINGER CCC EPOTE=0. DO 110 I=1,NOS N = NN(I) L = LL(I) LL1 = L*(L+1) EI = EHF(I) K = 0 S2 = 0. DO 300 J=1,NNN IF( LT(I).EQ.0) VQRE(J) = MESMN(J) * (AVP(J,I) - .25 * DHMEN(J) * 1 DHMEN(J) / HMEN(J) + .5*D2HMEN(J) ) + (1.-MESMN(J))*EI IF( LT(I).EQ.1) VQRE(J) = MESMP(J) * (AVP(J,I) - .25 * DHMEP (J)* 1 DHMEP(J) / HMEP(J) + .5*D2HMEP(J) ) + (1.-MESMP(J))*EI VQRE(J) = VQRE(J)/HB + LL1/(RM(J)*RM(J)) 300 CONTINUE GO TO 112 113 K=K+2 IF ( K- 200) 114,114,115 115 PRINT 204,ST(I),LT(I) GO TO 111 114 EI = -.5*FLOAT(K) 112 E = EI * DMSHB NO = N - 1 C CALL NUM1L(NNN,H,E,S2,VQRE,U,NO,.1E-03) IF (NO.LT.0) GO TO 113 111 CONTINUE CCC DO 12 J=1,NNN Y = SQRT( MESMN(J)*(1-LT(I)) + MESMP(J)*LT(I)) U(J) = U(J) * Y 12 FCT(J) = U(J) CALL DERIV S = 0. DO 303 J= 1,NNN DU(J) = DF(J) 303 S = S + U(J)*U(J) S = SQRT(S*H) SIG = U(1)/ABS(U(1)) DO 304 J=1,NNN DU(J) = SIG * DU(J)/S 304 U(J) = SIG*U(J)/S E = E/DMSHB EHF(I) = E DO 122 J=1,NNN DUNL(J,I) = DU(J) 122 UNL(J,I) = U(J) X = (E-EI)/E X = ABS(X) IF(ITER.EQ.ITMAX) PRINT 201,ST(I),LT(I),E,DEG(I),X IF(ITER.EQ.ITMAX.and.lt(i).eq.1.and.deg(i).ne.0.) write(3,*)e IF(ITER.EQ.ITMAX.and.lt(i).eq.0.and.deg(i).ne.0.) write(2,*)e EPOTE=EPOTE+E*DEG(I) 110 CONTINUE CCC CCC DENSITES 352 DO 143 J=1,NNN X = RM(J) R2 = RM2(J) R3 = RM3(J) E = 0. F = 0. E2= 0. F2 = 0. TE = 0. TF= 0. DO 142 I= 1,NOS L = LL(I) N = NN(I) JM = LJ(I) LL1 = L*(L+1) FL = (JM*(JM+2)-4 *L*(L+1) - 3.)/8. EVSO (I) = FL Y = UNL(J,I)*UNL(J,I)*DEG(I) E = E + Y*(1-LT(I)) F = F + Y*LT(I) EE = 2*FL*Y E2 = E2 + EE*(1-LT(I)) F2 = F2 + EE*LT(I) Y = (-UNL(J,I)/X + DUNL(J,I))/X Y1 = UNL(J,I)/R2 Y2 = DEG(I) *(Y*Y + LL1*Y1*Y1) TE = TE + Y2*(1-LT(I)) TF = TF + Y2 *LT(I) 142 CONTINUE DP(J) = F*UNS4PI/R2 DN(J) = E* UNS4PI/R2 DT(J) = DN(J) + DP(J) DJN(J) = E2*UNS4PI/R3 DJP(J) = F2*UNS4PI/R3 DJT(J) = DJN(J) + DJP(J) TAUN(J) = TE*UNS4PI TAUP(J) = TF * UNS4PI TAUT(J) = TAUN(J) + TAUP(J) 143 CONTINUE C CC DERIVEES DES DENSITES C DO 13 J=1,NNN 13 FCT(J) = DN(J) CALL DERIV DO 14 J=1,NNN DN1(J) = DF(J) 14 FCT(J) = DP(J) CALL DERIV DO 15 J=1,NNN DP1(J) = DF(J) DT1(J) = DN1(J) + DP1(J) 15 FCT(J) = DN1(J) CALL DERIV DO 16 J=1,NNN DN2(J) = DF(J) 16 FCT(J) = DP1(J) CALL DERIV DO 17 J=1,NNN DP2(J) = DF(J) DT2(J) = DP2(J) + DN2(J) 17 FCT(J) = DJN(J) CALL DERIV DO 18 J=1,NNN DJN1(J) = DF(J) 18 FCT(J) = DJP(J) CALL DERIV DO 19 J=1,NNN DJP1(J) = DF(J) 19 DJT1(J) = DJP1(J) + DJN1(J) C C CALCUL DES CHAMPS CCC TERME COULOMBIEN E=0. F = 0. DO 144 J=1,NNN X = RM(J) Y = QP*X*X E = E + Y*DP(J) F = F + Y*DP(J)/X 144 VC(J) = E/X - F DO 11 J=1,NNN 11 VC(J) = E22*(VC(J)+F) C DO 167 J = 1,NNN HMEN(J) =(HB + CRO*DT(J) + CROQ*DN(J) )*(1.-XMU) + HMEN(J) * XMU HMEP(J) =(HB + CRO*DT(J) + CROQ*DP(J) )*(1.-XMU)+HMEP(J)*XMU MESMN(J) = HB/HMEN(J) MESMP(J) = HB/HMEP(J) DHMEN(J) =(CRO*DT1(J) + CROQ*DN1(J) )*(1.-XMU) + DHMEN(J)*XMU DHMEP(J) =(CRO*DT1(J) + CROQ*DP1(J) )*(1.-XMU)+DHMEP(J)*XMU D2HMEN(J) =(CRO*DT2(J) + CROQ*DN2(J) )*(1.-XMU)+ D2HMEN(J) *XMU D2HMEP(J) =(CRO*DT2(J) + CROQ*DP2(J) )*(1.-XMU) +D2HMEP(J)*XMU 167 CONTINUE DO 1500 J=1,NNN VN(J) = TXRO*DT(J) + TXDRO *DT1(J)/RM(J) + TXD2RO*DT2(J) + 1 TXTAU * TAUT(J) - W3*(DJT(J)/RM(J) + DJT1(J)/2.) 1 +T3*0.0833333*((1.+0.5*X3)*DT(J)**(GAM+1.)*(2.+GAM) 2 -GAM*(X3+0.5)*(DN(J)**2+DP(J)**2)*DT(J)**(GAM-1.)) VSON(J) = W0* DT1(J)/(2.* RM(J)) - .125*( T1*X1 + T2*X2)*DJT(J) 1 /RM(J) VP (J) = VN(J) + TXROQ*DP(J) + TXDROQ* DP1(J)/RM(J) +TXD2RQ * 1 DP2(J) +TXTAUQ *TAUP(J) - W3*(DJP(J)/RM(J) +.5*DJP1(J)) 1 -T3*0.16666666*(X3+0.5)*DP(J)*DT(J)**GAM VSOP(J) = VSON(J) + W0*DP1(J)/(2.*RM(J)) DO 1502 I=1,NPO 1502 AVP(J,I) = (VP(J) +VSOP(J)*EVSO(I) +VC(J)+0.0*RM(J)**2)*(1. 1 -XMU)+XMU*AVP(J,I) VN(J) = VN(J) + TXROQ*DN(J) + TXDROQ*DN1(J)/RM(J) + TXD2RQ* DN2(J 1) + TXTAUQ * TAUN(J) - W3*( DJN(J)/RM(J) + .5*DJN1(J)) 1 -T3*0.16666666*(X3+0.5)*DN(J)*DT(J)**GAM VSON(J) = VSON(J) + W0*DN1(J)/(2.*RM(J)) NNEUT = NOS - NPO DO 1504 I2=1,NNEUT I = I2 + NPO 1504 AVP(J,I) = (VN(J) + VSON(J)*EVSO(I)-0.0*RM(J)**2)*(1.-XMU) 1 + XMU*AVP(J,I) 1500 CONTINUE IF (ITER. GE.ITMAX) GO TO 116 GO TO 104 CCC CCC FIN D ITERATION 116 CONTINUE do i=1,nnn write(88,*)rm(i),vn(i) write(89,*)rm(i),vp(i)+vc(i) enddo C CALCUL DES ENERGIES EREA=0. DO 1501 J=1,NNN 1501 EREA=EREA+GAM*((1.+0.5*X3)*DT(J)**2-(X3+0.5)*(DN(J)**2+ 1 DP(J)**2))*DT(J)**GAM*T3/24.*RM2(J) EREA=EREA*H*QP EREA=-EREA/MA ECTOT=0. DO 1510 J=1,NNN 1510 ECTOT=ECTOT+TAUT(J)*RM2(J) ECTOT=ECTOT*QP*H*HB EPOTE=EPOTE/MA ECTOT=ECTOT/MA H0=0.5*(EPOTE+ECTOT) ETOT=H0+EREA PRINT 217,ECTOT,H0,EREA,ETOT C CCCCC CALCUL DE LA DENSITE DE CHARGE DO 1503 J=1,NNN X=RM(J) 1503 DC(J)=DSC(X) CCCCCCCCC CALCUL DES RAYONS RN=0. RT = 0. RP=0. RC=0. DO 105 J=1,NNN X4=RM2(J)*RM2(J) RN=RN+DN(J)*X4 RT = RT + DT(J)*X4 RP=RP+DP(J)*X4 RC=RC+DC(J)*X4 105 CONTINUE RN=SQRT(QP*(H*RN/(MA-JZ))) RT = SQRT(QP*(H*RT/MA)) RP=SQRT(QP*(H*RP/JZ)) RC=SQRT(QP*(H*RC/JZ)) PRINT 221,RN,RP,RC PRINT 213 PRINT 212,(DN(J),J=1,NNN) PRINT 214 PRINT 212,(DP(J),J=1,NNN) PRINT 216 PRINT 212,(DC(J),J=1,NNN) C IF(IDCC.EQ.1) PUNCH1209,(DC(J),J=1,NNN) 1209 FORMAT(5E16.9) 156 CONTINUE DO 157 I=1,NOS PRINT 239,ST(I),LT(I) 239 FORMAT (5X,A6,I2/) PRINT 211,(UNL(J,I),J=1,NNN) PRINT 203 157 CONTINUE non=npo+1 do 158 j=1,nnn write(90,*)rm(j),(unl(j,jn),jn=non,nos) 158 continue PRINT 211,((UNL(J,3)/RM(J))**2/SQRT(QP),J=1,NNN) PRINT 203 GO TO 1 1000 CONTINUE STOP END SUBROUTINE NUM1L(N,H,E,S2,U,S,NO,EPS) CC VERSION CORRIGEE LE 21 NOV 72 C*****INTEGRATION DE L'EQUATION DE SCHROEDINGER PAR LA METHODE DE NUMERO C*****POUR E NEGATIF C*****RECHERCHE DE L'ENERGIE PROPRE PAR LA METHODE DE RAPHSON-NEWTON DIMENSION U(1),S(1) DATA RAP1,RAP2/0.,0./ H12=H*H/12. C*****CONTROLE DES CONDITIONS ASYMPTOTIQUES IF(E.GT..0) E=.0 DEI=.0 EPSS=.1E-10 IF(U(N-1).GT.EPSS) GO TO 10 DEI=U(N-1)-EPSS DO 8 K=1,N 8 U(K)=U(K)-DEI 10 U(N)=U(N-1) C*****CALCUL DU NOMBRE D'ETATS LIES PAR INTEGRATION A ENERGIE NULLE S(1)=1.E-10 B0=0. AA=H12*U(1) IF (S2) 16,18,16 16 B0=-S(1)*AA 18 B1=S(1)*(1.-AA) DO 38 K=2,N B2=12.*S(K-1)-10.*B1-B0 IF (ABS(B2).LT.1.E+10) GO TO 22 B2=B2*1.E-20 B1=B1*1.E-20 22 AA=H12*U(K) S(K)=B2/(1.-AA) B0=B1 38 B1=B2 DO 42 K=5,N N0=K IF(U(K).LT.0.) GO TO 44 42 CONTINUE 44 NEL=0 DO 52 K=N0,N IF (S(K-1)*S(K)) 46,50,52 46 NEL=NEL+2 GO TO 52 50 NEL=NEL+1 52 CONTINUE NEL=NEL/2 IF(NEL.GT.NO) GO TO 64 IF(NEL.EQ.NO) GO TO 60 62 NO=-1 RETURN 60 RAP1=S(N-1)/S(N) RAP2=EXP(H*SQRT(U(N-1)-E)) IF(RAP1.LT.RAP2) GO TO 62 C*****CALCUL DE EMIN ET EMAX ENTRE LESQUELLES SE TROUVE L'ENERGIE PROPRE 64 UMIN=U(1) DO 70 K=2,N IF(U(K).LT.UMIN) UMIN=U(K) 70 CONTINUE EMIN=UMIN EMAX=0. C*****DEBUT DE LA RECHERCHE DE L'ENERGIE PROPRE DANS L'INTERVALLE MAXIMU TE=EMAX-EMIN C*****REJET DE L'ENERGIE D'ESSAI E PROPOSEE SI ELLE EST A L'EXTERIEUR DE C*****BORNES (EMIN,EMAX) IF((E.LT.EMIN).OR.(E.GT.EMAX)) E=EMIN+.5*TE E1=EMIN E2=EMAX J=2 I=1 GO TO 102 C*****REDUCTION DES BORNES EMIN ET EMAX 90 EMIN=E1 EMAX=E2 TE=EMAX-EMIN J=2 98 I=1 100 E=EMIN+TE*FLOAT(I)/FLOAT(J) 102 DE=0. 104 E=E+DE IF(E.GT.0.) GO TO 204 S(N)=1.E-10 N1=N-1 RAP2=EXP(H*SQRT((U(N-1)+U(N))/2.-E)) S(N1)=S(N)*RAP2 AA=H12*(U(N1)-E) B0=S(N)*(1.-AA) B1=S(N1)*(1.-AA) N1=N-2 DO 138 KAUX=1,N1 K=N1-KAUX+1 B2=12.*S(K+1)-10.*B1-B0 AA=H12*(U(K)-E) S(K)=B2/(1.-AA) B0=B1 B1=B2 IF(U(K).LT.E) GO TO 140 138 CONTINUE 140 N1=K C*****NORMALISATION DE LA FONCTION D'ONDE A S(N1) DO 146 KAUX=N1,N K=N-KAUX+N1 146 S(K)=S(K)/S(N1) C*****DEBUT DE L'INTEGRATION VERS L'EXTERIEUR JUSQU'A N1 S(1)=1.E-10 B0=0. AA=H12*(U(1)-E) IF(S2) 156,158,156 156 B0=-S(1)*AA 158 B1=S(1)*(1.-AA) DO 170 K=2,N1 B2=12.*S(K-1)-10.*B1-B0 AA=H12*(U(K)-E) S(K)=B2/(1.-AA) B0=B1 170 B1=B2 C*****NORMALISATION DE LA FONCTION A S(N1) DO 174 K=1,N1 174 S(K)=S(K)/S(N1) C*****CALCUL DE LA CORECTION D'ENERGIE SOM=0. DO 180 K=1,N 180 SOM=SOM+S(K)*S(K) DE=((-S(N1-1)+2.-S(N1+1))/(H*H)+U(N1)-E)/SOM IF(ABS(DE).GT.EPS) GO TO 104 C*****CALCUL DU NOMBRE DE NOEUDS DE L'ETAT PROPRE TROUVE DO 182 K=5,N IF(U(K).LT.E) GO TO 184 182 CONTINUE 184 N0=K NEL=0 DO 192 K=N0,N1 IF(S(K-1)*S(K)) 186,190,192 186 NEL=NEL+2 GO TO 192 190 NEL=NEL+1 192 CONTINUE NEL=NEL/2 C*****L'ETAT PROPRE TROUVE EST-IL LE BON IF(NEL-NO) 198,214,202 198 IF(E.GT.E1) E1=E GO TO 204 202 IF(E.LT.E2) E2=E 204 I=I+2 IF (I.LE.J) GO TO 100 J=2*J IF(ABS(E1-EMIN).GT.EPS.OR.ABS(EMAX-E2).GT.EPS) GO TO 90 GO TO 98 C*****NORMALISATION DE LA FONCTION PROPRE 214 SOM=1./SQRT(SOM*H) DO 218 K=1,N 218 S(K)=S(K)*SOM E=E+DEI RETURN C*****DEBUT FORMATS 2000 FORMAT(/,35X,56HL@ETAT DEMANDE N@EST PAS LIE. RETOUR DE NUM1L AVEC 1 NO=-1,/) C*****FIN FORMATS END FUNCTION DSC(Z1) CCC CCC DENSITE DE CHARGE CCC FACTEUR DE FORME GAUSSIEN CCC MU=0.65 FM CCC DIMENSION YK(10),PI(10),SIG(2) DATA (YK(I),I=1,10)/0.24534071,0.73747373,1.2340762,1.7385377,2.25 149740,2.7888061,3.3478546,3.9447640,4.6036824,5.3874809/ DATA(PI(I),I=1,10)/0.49092150,0.49384339,0.49992087,0.50967903, 10.52408035,0.54485174,0.57526244,0.62227870,0.70433296,0.89859196/ DATA SIG/1.,-1./ X=Z1 XMU=0.65 ZW=1./(0.4431125*XMU**3) SG=0. XSM=X/XMU DO 100 I=1,10 DO 100 K=1,2 Y=X+XMU*YK(I)*SIG(K) IF(Y)100,100,101 101 F=ZW*GUGUS(Y)*Y*Y YSM=Y/XMU SG=SG+XMU*PI(I)*F*VKG(0,XSM,YSM) 100 CONTINUE DSC=SG RETURN END FUNCTION VKG(I1,Z1,Z2) DIMENSION FJ(50) L=I1 X=Z1 Y=Z2 DELTA=Y-X A=2.*X*Y IF(A-15.)101,101,113 101 IF(A-0.01)102,102,104 102 ARG=EXP(-(X*X+Y*Y)) FJ0=ARG*(1.+A*A/6.) IF(L-1)103,103,105 103 FJ(1)=FJ0 FJ(2)=ARG*A*(10.+A*A)/30. GO TO 107 104 U=DELTA*DELTA V=(X+Y)*(X+Y) U=EXP(-U) IF(V.GT.670.) V=669. V=EXP(-V) FJ0=(U-V)/(2.*A) 105 L2=L+5 L2=L+10 FJ(L2)=1.E-10*FLOAT(2*L2+1)/A FJ(L2+1)=1.E-10 L3=L2-1 DO 106 LL=1,L3 L1=L2-LL FJ(L1)=FLOAT(2*L1+1)*FJ(L1+1)/A+FJ(L1+2) IF(FJ(L1)-1.E+30)106,106,111 111 DO 112 L4=L1,L2 112 FJ(L4)=1.E-10*FJ(L4) 106 CONTINUE ZZ=FJ0/FJ(1) L2=L2-9 DO 109 L1=1,L2 109 FJ(L1)=ZZ*FJ(L1) GO TO 107 113 U=DELTA*DELTA V=(X+Y)*(X+Y) U=EXP(-U) IF(V.GT.670.) V=669. V=EXP(-V) FJ(1)=(U-V)/(2.*A) FJ(2)=((A-1.)*U+(A+1.)*V)/(2.*A*A) IF(L-1)107,107,114 114 DO 115 L1=2,L 115 FJ(L1+1)=-FLOAT(2*L1-1)*FJ(L1)/A+FJ(L1-1) 107 VKG=FLOAT(L+L+1)*FJ(L+1) RETURN END FUNCTION HWF(I1,I2,Z1,Z2) N=I1 L=I2 B=Z1 R=Z2 X=R/B XSQ=-2.*X*X Y=-X*X/2. B3=B*B*B RHO=1. VNL=1. IF(N)101,101,102 102 XA=1. K=2*L+1 NS=N+1 NI=0 DO 103 I=1,N K=K+2 NI=NI+1 NS=NS-1 XA=XA*XSQ*NS/(K*NI) VNL=VNL+XA 103 RHO=(K*RHO)/(2.*NI) 101 DPL=1. XSL=1. IF(L)104,104,105 105 K=1 DO 106 I=1,L K=K+2 XK=K DPL=2.*DPL/XK 106 XSL=XSL*X 104 ARG=RHO*DPL/(B3*1.772454) XA=2.*XSL*R*SQRT(ARG) HWF=XA*VNL*EXP(Y) RETURN END SUBROUTINE DERIV C FIVE POINT DIFFERENTATION FORMULA (ABRAMOWITZ, PAG 914) COMMON /DER/ FCT(100),DF(100),H,NNN DIMENSION A(5,5) DATA A(1,1),A(1,2),A(1,3),A(1,4),A(1,5) + /-50.,96.,-72.,32.,-6./ DATA A(2,1),A(2,2),A(2,3),A(2,4),A(2,5) + /-6.,-20.,36.,-12.,2./ DATA A(3,1),A(3,2),A(3,3),A(3,4),A(3,5) + /2.,-16.,0.,16.,-2./ DATA A(4,1),A(4,2),A(4,3),A(4,4),A(4,5) + /-2.,12.,-36.,20.,6./ DATA A(5,1),A(5,2),A(5,3),A(5,4),A(5,5) + /6.,-32.,72.,-96.,50./ DATA EMFACT/24./ NMX=NNN-2 DO 1 J=1,NNN K=3 IF(J.LT.3)K=J IF(J.GT.NMX)K=J-NNN+5 SUM=0. DO 2 I=1,5 JJ=J+I-K 2 SUM=SUM+A(K,I)*FCT(JJ) 1 DF(J)=SUM/(H*EMFACT) RETURN END FUNCTION GUGUS(X) C---- INTERPOLATION QUADRATIQUE DE FONCTIONS PAIRES SUR 0-RMAX COMMON /DER/ FCT(100),DF(100),H,NNN COMMON RM(100),DU(100),VC(100),U(100),VS(100), + TAUN(100),TAUP(100),TAUT(100),DJN1(100),DJP1(100), +DJT1(100),VN(100),VP(100),VSON(100),VSOP(100), + ECOULB(100),VQRE(100),RM2(100),RM3(100),HMEN(100),HMEP(100), +DHMEN(100),DHMEP(100),D2HMEN(100),D2HMEP(100),DN(100),DP(100), +DT(100),DN1(100),DP1(100),DT1(100),DN2(100),DP2(100),DT2(100), +DJN(100),DJP(100),DJT(100),DC(100) + ,AVP(100 ,45),DUNL(100 ,45),UNL(100 ,45) X2=2.*H IF(X.GT.X2) GO TO 2 GUGUS=((4.*DP(1)-DP(2))+(DP(2)-DP(1))*X*X)/3. RETURN 2 I=X/H IF(I.GE.NNN) GO TO 4 Y=(X-I*H)/H S=DP(I)+.5*Y*(DP(I+1)-DP(I-1)) GUGUS=S+.5*Y*Y*(DP(I+1)+DP(I-1)-2.*DP(I)) RETURN 4 GUGUS=.0 RETURN END