From 4dd9e71ed5ec358f288c5e7233fd109918a722df Mon Sep 17 00:00:00 2001 From: "Ryan@Home" Date: Wed, 12 Feb 2025 23:34:08 -0500 Subject: [PATCH 1/3] add comments in ADWUCK4.FOR and BDWUCK.FOR --- dwuck4/ADWCK4.FOR | 4 +-- dwuck4/BDWCK4.FOR | 89 +++++++++++++++++++++++++++++------------------ 2 files changed, 57 insertions(+), 36 deletions(-) diff --git a/dwuck4/ADWCK4.FOR b/dwuck4/ADWCK4.FOR index 1ecb4e3..dc863f5 100644 --- a/dwuck4/ADWCK4.FOR +++ b/dwuck4/ADWCK4.FOR @@ -380,8 +380,8 @@ C Q=E E=(ECM(1)+Q)*(FM(2)+FMA(2))/FMA(2) ENDIF - IS(N)=FS(N) - NS(N)=IS(N)+1 + IS(N)=FS(N) ! 2*spin, FS(N) will divided by 2 later + NS(N)=IS(N)+1 ! 2*spin+1, number of m-state IF(AMASS.EQ.0.0) AMASS=FMA(1) IF(IK.EQ.0) DR(N)=DRF*AMASS/FMA(N) KMXX=KMAX diff --git a/dwuck4/BDWCK4.FOR b/dwuck4/BDWCK4.FOR index 92978a2..64592c2 100644 --- a/dwuck4/BDWCK4.FOR +++ b/dwuck4/BDWCK4.FOR @@ -269,48 +269,56 @@ c*********************************************************************** 5 ,(DTEMP(1201),GP),(DTEMP(1601),S ),(DTEMP(2001),C ) DATA ETA3/10.E+00/ C + IWORD=0 - JT=NS(1)+NS(2) - NP=LPL2*JT + JT=NS(1)+NS(2) ! NS(1) is number of m-state in incoming channel, similar for NS(2), for (d,p), JT = 5 + NP=LPL2*JT ! LPLUS = LMAX + 1, LPL2 = 2*LPLUS, NP = 2*(LMAX+1)*JT, for LMAX = 15 (d,p) NP = 160, NP= number of partial wave, real + imag I=0 - DO 30 N=1,2 - DR2(N)=DR(N)**2/12.0 + DO 30 N=1,2 !N = 1 : incoming, 2 : outgoing + DR2(N)=DR(N)**2/12.0 ! this is for the Numerov method R(N)=0.0 - JS=NS(N) - DO 29 ISS=1,JS - I=I+1 - LM(I)=0 + JS=NS(N) ! number of m-state + DO 29 ISS=1,JS ! loop all m-state + I=I+1 + LM(I)=0 ! all LM(I) = 0 29 CONTINUE - 30 CONTINUE - DO 40 IQ=1,NP - F1(IQ)=0.0 - F2(IQ)=0.0 + 30 CONTINUE ! end of loop N + DO 40 IQ=1,NP ! initial distorted wave, have NP points + F1(IQ)=0.0 ! + F2(IQ)=0.0 ! F1 is n-1 of F2 Q1(IQ)=0.0 - Q2(IQ)=0.0 + Q2(IQ)=0.0 ! Q1 is n-1 of Q2 40 CONTINUE C - DO 100 M=1,K - MK=M+M-1 - IX=0 - I=0 - DO 90 N=1,2 + + WRITE(6,*) 'Debug: K=', K, ', NP=', NP + DO 100 M=1,K ! K = ABS(RMAX)/DRF + 1.0E-08 + MK=M+M-1 ! 2*M-1, odd number from 1 to K, odd for real, even for imag + IX=0 ! loop from 1 to 128, step 32 = 2*LPLUS, why? + I=0 ! loop from 1 to 5, total m-state in incoming and outgoing channel? + + DO 90 N=1,2 ! looping incoming and outgoing channel R(N)=R(N)+DR(N) - DRR2(N)=DR2(N)/R(N)**2 - Q(1)=1.0+DR2(N)*U(MK ,N) - Q(2)= DR2(N)*U(MK+1,N) - LTEMP=2.0*FK(N)*R(N)+ETA3 - LTEMP=MIN0(LTEMP,LPLUS) - FI=-FS(N) - JS=NS(N) - SFACT=FS(N)**2+FS(N) - DO 89 ISS=1,JS - I=I+1 - FL=0.0 - DO 80 LL=1,LPLUS - FJ=FL+FI - IX1=IX+LL+LL-1 - FLFACT=FL**2+FL - FACT=DR2(N)*(FJ**2+FJ-FLFACT-SFACT)*0.5 + DRR2(N)=DR2(N)/R(N)**2 ! for L(L+1)/r^2 term in the potential + Q(1)=1.0+DR2(N)*U(MK ,N) ! seems to be the Numerov k_n + Q(2)= DR2(N)*U(MK+1,N) + LTEMP=2.0*FK(N)*R(N)+ETA3 ! ETA3 = 10, + LTEMP=MIN0(LTEMP,LPLUS) ! set the maximum acceptable L + + FI=-FS(N) ! m-state of S of channle-N, N=1 for incoming, =2 for outgoing + JS=NS(N) ! number of m-state of S + SFACT=FS(N)**2+FS(N) ! s * (s+1) + + DO 89 ISS=1,JS ! loop the m-state + I=I+1 ! I is the id of m-state + FL=0.0 ! FL runs from 0 to LMAX + + DO 80 LL=1,LPLUS ! loop all L, fortan start index is 1, so need to run from 1 to LMAX + 1 + FJ=FL+FI ! J = L + S, looping possible J-state + IX1=IX+LL+LL-1 ! loop from 1 to 159 odd, odd for real? even for imag? + FLFACT=FL**2+FL + + FACT=DR2(N)*(FJ**2+FJ-FLFACT-SFACT)*0.5 ! for L(L+1)/r^2 Q(3 )=Q(1)+FACT*V(MK ,N)-DRR2(N)*FLFACT Q(4 )=Q(2)+FACT*V(MK+1,N) IF(LL.LE.LM(I)) GO TO 70 @@ -318,10 +326,12 @@ C LM(I)=LM(I)+1 IF(FJ-ABS(FL-FS(N)).LT.0.0) GO TO 72 c calculate approximate starting value + f2(ix1 )=1.0 do 50 ii=1,ll f2(ix1 )=f2(ix1 )*(fk(n)*r(n))/float(2*ii-1) 50 continue + F2(IX1+1)=0.0 Q2(IX1 )=Q(3)*f2(ix1 ) Q2(IX1+1)=Q(4)*f2(ix1 ) @@ -345,6 +355,12 @@ c Q1(IX1+1)=Q2(IX1+1) Q2(IX1 )=CTEMP(1) Q2(IX1+1)=CTEMP(2) + + IF(N.EQ.1 .AND. FL.EQ.1 .AND. FJ.EQ.1) THEN + WRITE(6,'(A, F7.3, F10.6, F10.6, F10.6, F10.6)') + 1'Debug:',R(N),Q1(IX1),Q1(IX1+1),Q2(IX1), Q2(IX1+1) + ENDIF + 72 CONTINUE FL=FL+1.0 80 CONTINUE @@ -357,8 +373,13 @@ C WRITE RADIAL WAVE FUNCTIONS ON TAPE 4 C WRITE(4)(F2(J),J=1,NP) 100 CONTINUE + +C ================= end of loop +C + LX=1 drrc = 0.1 + DO 120 N=1,2 R2=FK(N)*R(N) R1=R2-DR(N)*FK(N) From 359b7d0defe0f09ed8c33c38ae0a194ee4c249ec Mon Sep 17 00:00:00 2001 From: "Ryan@Home" Date: Thu, 13 Feb 2025 21:46:55 -0500 Subject: [PATCH 2/3] add comments to BDWUCK4.FOR --- dwuck4/BDWCK4.FOR | 157 +++++++++++++++++++++++++++++++--------------- 1 file changed, 106 insertions(+), 51 deletions(-) diff --git a/dwuck4/BDWCK4.FOR b/dwuck4/BDWCK4.FOR index 64592c2..2f0d787 100644 --- a/dwuck4/BDWCK4.FOR +++ b/dwuck4/BDWCK4.FOR @@ -271,19 +271,29 @@ c*********************************************************************** C IWORD=0 - JT=NS(1)+NS(2) ! NS(1) is number of m-state in incoming channel, similar for NS(2), for (d,p), JT = 5 + JT=NS(1)+NS(2) ! NS(1) is number of J-state in incoming channel, similar for NS(2), for (d,p), JT = 5 NP=LPL2*JT ! LPLUS = LMAX + 1, LPL2 = 2*LPLUS, NP = 2*(LMAX+1)*JT, for LMAX = 15 (d,p) NP = 160, NP= number of partial wave, real + imag I=0 - DO 30 N=1,2 !N = 1 : incoming, 2 : outgoing +c-------------------------------------- loop incoming and outgoing channel + DO 30 N=1,2 !N = 1 : incoming channel, 2 : outgoing DR2(N)=DR(N)**2/12.0 ! this is for the Numerov method R(N)=0.0 - JS=NS(N) ! number of m-state - DO 29 ISS=1,JS ! loop all m-state + JS=NS(N) ! number of J-state in N-channel +c------------------------------- loop all J-state + DO 29 ISS=1,JS ! loop all J-state I=I+1 LM(I)=0 ! all LM(I) = 0 29 CONTINUE +c------------------------------- end of J_state loop 30 CONTINUE ! end of loop N - DO 40 IQ=1,NP ! initial distorted wave, have NP points +c---------------------------------------- end of channel loop +c Nemerov method +c y''(r) = g(r) y(r) +c k(n+1)y(n+1) = (12 - 10*k(n))*y(n+1) - k(n-1)*y(n-1) +c k(n) = 1 + dr^2/12 * g(n) +c F2 = distorted wave = y(n) +c Q2 = k(n)*y(n) + DO 40 IQ=1,NP ! initial distorted wave, have NP partial wave F1(IQ)=0.0 ! F2(IQ)=0.0 ! F1 is n-1 of F2 Q1(IQ)=0.0 @@ -292,90 +302,126 @@ C C WRITE(6,*) 'Debug: K=', K, ', NP=', NP +c============================================== loop radial grids DO 100 M=1,K ! K = ABS(RMAX)/DRF + 1.0E-08 MK=M+M-1 ! 2*M-1, odd number from 1 to K, odd for real, even for imag - IX=0 ! loop from 1 to 128, step 32 = 2*LPLUS, why? - I=0 ! loop from 1 to 5, total m-state in incoming and outgoing channel? - + IX=0 ! loop from 1 to 128, step 32 = 2*LPLUS, each group of 32 is all L-states for a J-state + I=0 ! loop from 1 to 5, total J-state in incoming and outgoing channel +c------------------------------------ loop channels DO 90 N=1,2 ! looping incoming and outgoing channel R(N)=R(N)+DR(N) DRR2(N)=DR2(N)/R(N)**2 ! for L(L+1)/r^2 term in the potential - Q(1)=1.0+DR2(N)*U(MK ,N) ! seems to be the Numerov k_n - Q(2)= DR2(N)*U(MK+1,N) - LTEMP=2.0*FK(N)*R(N)+ETA3 ! ETA3 = 10, + Q(1)=1.0+DR2(N)*U(MK ,N) ! for real part , seem to be the Numerov k_n + Q(2)= DR2(N)*U(MK+1,N) ! for imag part + LTEMP=2.0*FK(N)*R(N)+ETA3 ! ETA3 = 10, this is the theoritic maximum L LTEMP=MIN0(LTEMP,LPLUS) ! set the maximum acceptable L - - FI=-FS(N) ! m-state of S of channle-N, N=1 for incoming, =2 for outgoing - JS=NS(N) ! number of m-state of S + FI=-FS(N) ! s-state of S of channle-N, N=1 for incoming, =2 for outgoing + JS=NS(N) ! number of J-state of S SFACT=FS(N)**2+FS(N) ! s * (s+1) - - DO 89 ISS=1,JS ! loop the m-state - I=I+1 ! I is the id of m-state +c-------------------------- loop J-state + DO 89 ISS=1,JS ! loop the J-state + I=I+1 ! I is the id of J-state FL=0.0 ! FL runs from 0 to LMAX - +c---------------- loop the L-state DO 80 LL=1,LPLUS ! loop all L, fortan start index is 1, so need to run from 1 to LMAX + 1 FJ=FL+FI ! J = L + S, looping possible J-state - IX1=IX+LL+LL-1 ! loop from 1 to 159 odd, odd for real? even for imag? + IX1=IX+LL+LL-1 ! index in memomry, loop from 1 to 159 odd, odd for real? even for imag? FLFACT=FL**2+FL - FACT=DR2(N)*(FJ**2+FJ-FLFACT-SFACT)*0.5 ! for L(L+1)/r^2 Q(3 )=Q(1)+FACT*V(MK ,N)-DRR2(N)*FLFACT Q(4 )=Q(2)+FACT*V(MK+1,N) IF(LL.LE.LM(I)) GO TO 70 IF(LTEMP.LT.LL) GO TO 72 - LM(I)=LM(I)+1 - IF(FJ-ABS(FL-FS(N)).LT.0.0) GO TO 72 -c calculate approximate starting value - + LM(I)=LM(I)+1 !this control the calculateing of start value, weird but work. + IF(FJ-ABS(FL-FS(N)).LT.0.0) GO TO 72 ! j < ||l-s|, increase FL by 1 +c........... calculate approximate starting value +c for FL < 9, R=0.1 is set +c for FL = 10, R=0.5 is set +c FL = 11, R=0.89 is set +c FL = 12, R=1.3 +c FL = 13, R=1.7 +c FL = 14, R=2.1 +c FL = 15, R=2.5 f2(ix1 )=1.0 do 50 ii=1,ll f2(ix1 )=f2(ix1 )*(fk(n)*r(n))/float(2*ii-1) 50 continue +c IF(N.EQ.1 ) THEN +c WRITE(6,*) R(N), FL, FJ, IX, IX1, f2(ix1), M +c ENDIF + +c IF(N.EQ.1 .AND. FL.EQ.1 .AND. FJ.EQ.1) THEN +c WRITE(6,5678) +c 1'Debug:',R(N),Q1(IX1),Q2(IX1),F2(IX1), +c 2 Q1(IX1+1),Q2(IX1+1),F2(IX1+1), +c 3 Q(3), Q(4) +c ENDIF + F2(IX1+1)=0.0 Q2(IX1 )=Q(3)*f2(ix1 ) Q2(IX1+1)=Q(4)*f2(ix1 ) C C EVALUATE Q AT ORIGIN FOR L=1 C - IF(LL.EQ.2) Q1(IX+3)=-f2(ix1 )/6.0 + IF(LL.EQ.2) Q1(IX+3)=-f2(ix1 )/6.0 ! when LL is 2 or FL = 1 GO TO 72 +c........... end of starting value 70 CONTINUE c c Step equations forward by dr(n) via Numerov-Fox-Goodwin-Milne method c - CTEMP(1)=12.*F2(IX1 )-10.*Q2(IX1 )-Q1(IX1 ) + +c IF(N.EQ.1 .AND. FL.EQ.1 .AND. FJ.EQ.1) THEN +c WRITE(6,5678) +c 1'Debug:',R(N),Q1(IX1),Q2(IX1),F2(IX1), +c 2 Q1(IX1+1),Q2(IX1+1),F2(IX1+1), +c 3 Q(3), Q(4) +c 5678 FORMAT(A, F7.3, F10.6, F10.6, F10.6, F10.6, +c 1 F10.6, F10.6, F10.6, F10.6) +c ENDIF +c +c Q2 (n+1) = 12 * y (n) - 10 * Q2 (n) - Q2 (n-1) for real +c Q2'(n+1) = 12 * y'(n) - 10 * Q2'(n) - Q2'(n-1) for imag +c a(n) = Q(3) = k(n) for real ? +c b(n) = Q(4) = k(n) for imag ? +c y (n+1) = [ Q2(n+1)*a(n) + Q2'(n+1)*b(n) ] / (a(n)^2 + b(n)^2) +c y'(n+1) = [ -Q2(n+1)*b(n) + Q2'(n+1)*a(n) ] / (a(n)^2 + b(n)^2) +c + CTEMP(1)=12.*F2(IX1 )-10.*Q2(IX1 )-Q1(IX1 ) ! k(n+1)y(n+1) = (12 - 10*k(n))*y(n+1) - k(n-1)*y(n-1) CTEMP(2)=12.*F2(IX1+1)-10.*Q2(IX1+1)-Q1(IX1+1) - F1(IX1 )=F2(IX1 ) + F1(IX1 )=F2(IX1 ) ! save the old f2 to f1 F1(IX1+1)=F2(IX1+1) - DET=Q(3)**2+Q(4)**2 - F2(IX1 )=(CTEMP(1)*Q(3 )+CTEMP(2)*Q(4 ))/DET + DET=Q(3)**2+Q(4)**2 ! real^2 + imag^2 + F2(IX1 )=(CTEMP(1)*Q(3 )+CTEMP(2)*Q(4 ))/DET ! new f2 = F2(IX1+1)=(CTEMP(2)*Q(3 )-CTEMP(1)*Q(4 ))/DET Q1(IX1 )=Q2(IX1 ) Q1(IX1+1)=Q2(IX1+1) Q2(IX1 )=CTEMP(1) Q2(IX1+1)=CTEMP(2) - IF(N.EQ.1 .AND. FL.EQ.1 .AND. FJ.EQ.1) THEN - WRITE(6,'(A, F7.3, F10.6, F10.6, F10.6, F10.6)') - 1'Debug:',R(N),Q1(IX1),Q1(IX1+1),Q2(IX1), Q2(IX1+1) - ENDIF +c IF(N.EQ.1 .AND. FL.EQ.1 .AND. FJ.EQ.1) THEN +c WRITE(6,'(A, F7.3, F10.6, F10.6, F10.6, F10.6, F10.6, F10.6)') +c 1'Debug:',R(N),Q1(IX1),Q2(IX1),F2(IX1), +c 2 Q1(IX1+1),Q2(IX1+1),F2(IX1+1) +c ENDIF 72 CONTINUE FL=FL+1.0 80 CONTINUE +c---------------- end of loop the L-state FI=FI+1.0 - IX=IX+LPL2 + IX=IX+LPL2 ! after L-state loop, go to next index 89 CONTINUE +c-------------------------- end of loop J-state 90 CONTINUE +c------------------------------------ end of loop channels C C WRITE RADIAL WAVE FUNCTIONS ON TAPE 4 C WRITE(4)(F2(J),J=1,NP) 100 CONTINUE - -C ================= end of loop -C +c============================================== end of loop radial grids LX=1 drrc = 0.1 @@ -404,8 +450,8 @@ c DO 200 N=1,2 JS=NS(N) FI=-FS(N) - ARG=S(LX)-S(LX-LL+1) - Q(1)=COS(ARG) + ARG=S(LX)-S(LX-LL+1) ! Coulomb phase shift ? + Q(1)=COS(ARG) Q(2)=SIN(ARG) Q(3)=Q(1)**2-Q(2)**2 Q(4)=2.0*Q(1)*Q(2) @@ -413,14 +459,14 @@ c FJ=FL+FI I=I+1 DET=F(LX)*GP(LX)-FP(LX)*G(LX) - A(1)=(F1(IX1 )*GP(LX)-F2(IX1 )*G (LX))/DET - A(2)=(F1(IX1+1)*GP(LX)-F2(IX1+1)*G (LX))/DET - B(1)=(F2(IX1 )*F (LX)-F1(IX1 )*FP(LX))/DET - B(2)=(F2(IX1+1)*F (LX)-F1(IX1+1)*FP(LX))/DET + A(1)=(F1(IX1 )*GP(LX)-F2(IX1 )*G (LX))/DET ! real part of F + A(2)=(F1(IX1+1)*GP(LX)-F2(IX1+1)*G (LX))/DET ! imag part of F + B(1)=(F2(IX1 )*F (LX)-F1(IX1 )*FP(LX))/DET ! real part of G + B(2)=(F2(IX1+1)*F (LX)-F1(IX1+1)*FP(LX))/DET ! imag part of G IF(LL.LE.LM(I).and.FJ-ABS(FL-FS(N)).ge.0.0) then - DET=(A(1)+B(2))**2+(A(2)-B(1))**2 - CTEMP(1)=(A(1)+B(2))/DET - CTEMP(2)=(B(1)-A(2))/DET + DET=(A(1)+B(2))**2+(A(2)-B(1))**2 ! this is the normalization + CTEMP(1)=(A(1)+B(2))/DET ! CTEMP(1) = cos(sigma)/sqrt(DET) + CTEMP(2)=(B(1)-A(2))/DET ! = - sin(sigma)/sqrt(DET) else CTEMP(1)=0.0 CTEMP(2)=0.0 @@ -428,13 +474,22 @@ c C C C=NORMALIZATION CONSTANTS C - C(IX1 )=Q(1)*CTEMP(1)-Q(2)*CTEMP(2) - C(IX1+1)=Q(1)*CTEMP(2)+Q(2)*CTEMP(1) + C(IX1 )=Q(1)*CTEMP(1)-Q(2)*CTEMP(2) ! cos(sigma)^2 + sin(sigma)^2 = 1/ sqrt(DET) + C(IX1+1)=Q(1)*CTEMP(2)+Q(2)*CTEMP(1) ! 0 ? + C -C E=PARTIAL WAVE SCATTERING AMPLITUDES +C E=PARTIAL WAVE SCATTERING AMPLITUDES, ScatAmp = (S-1)/2i C - E(2*I-1)=B(1)*CTEMP(1)-B(2)*CTEMP(2) - E(2*I )=B(1)*CTEMP(2)+B(2)*CTEMP(1) + E(2*I-1)=B(1)*CTEMP(1)-B(2)*CTEMP(2) ! real part of ScatAmp = s2/2 + E(2*I )=B(1)*CTEMP(2)+B(2)*CTEMP(1) ! image = (1-s1)/2 + +c IF(N.EQ.1) THEN +c WRITE(6,*) 'Norm ', FL, FJ, A(1), A(2), B(1), B(2), +c 1 CTEMP(1), CTEMP(2), +c 2 ARG, C(IX1), C(IX1+1), 1./DET, +c 3 E(2*I-1), E(2*I) +c ENDIF + T1 = E(2*I-1) T2 = E(2*I ) if(isym(N) .and. is(N).eq.0 ) then From 299cefb9cec193361e6e4c3146d6a7b816dd001b Mon Sep 17 00:00:00 2001 From: "Ryan@Home" Date: Sun, 16 Feb 2025 17:24:22 -0500 Subject: [PATCH 3/3] snapshot for DWUCK4 --- dwuck4/CDWCK4.FOR | 112 +++++++++++++++++++++++++++++++-------- dwuck4/DWtest2.DAT | 16 ++++++ dwuck4/LGNDR.py | 66 +++++++++++++++++++++++ dwuck4/culib8/LGNDR.FOR | 30 +++++++---- dwuck4/culib8/POLFCT.FOR | 3 ++ 5 files changed, 194 insertions(+), 33 deletions(-) create mode 100644 dwuck4/DWtest2.DAT create mode 100755 dwuck4/LGNDR.py diff --git a/dwuck4/CDWCK4.FOR b/dwuck4/CDWCK4.FOR index be67e0a..b281d97 100644 --- a/dwuck4/CDWCK4.FOR +++ b/dwuck4/CDWCK4.FOR @@ -537,12 +537,13 @@ c igam2=.false. endif c - MPLUS=JX/2+1 - IFACT=MPLUS*JR*JS + MPLUS=JX/2+1 ! j/2 + 1, spin transfer + 1 + IFACT=MPLUS*JR*JS ! (j+1)*num of m-state of sa * num of m-state of sb J2K=(1.0+PHASEF(NS(2)))/2.0 M2K=JX-MPLUS-MPLUS+2 - LX=LTR+LTR - TEMP2=SQRT(FLOAT((JX+1)*(IS(1)+1)))*flfact + LX=LTR+LTR + TEMP2=SQRT(FLOAT((JX+1)*(IS(1)+1)))*flfact ! IS(1) = 2*sa, JX = 2*j, sqrt((2j+1)*(2sa+1)) * flfact, for (d,p), flfact = 100 * sqrt((2l+1)/(2j+1)) +c for (d,p), TEMP2 = 100 * sqrt((2sa+1)(2l+1)) IF(FN.EQ.0.0) THEN c clear amplitude storage unless for coherent sum IND=2*LPLUS*IFACT @@ -550,9 +551,12 @@ c clear amplitude storage unless for coherent sum D(M)=0.0 10 CONTINUE ENDIF - IS1=-IS(1) + + write(6,*) 'Debug', IND + + IS1=-IS(1) ! m-state of sa ? step +2 DO 95 I=1,JR - IS2=-IS(2) + IS2=-IS(2) ! m-state of sb ? step +2 DO 90 J=1,JS IF(NLTR.NE.1) GO TO 14 IF(JR*JS.EQ.1) GO TO 15 @@ -565,10 +569,10 @@ C READ (2)(FLL(INDEX),INDEX=1,INCR) 15 continue c final L loop - DO 80 LL=1,LPLUS - lf=LL-1 - LLX=lf+lf - JLX=LLX+IS2 + DO 80 LL=1,LPLUS ! loop on Lb + lf=LL-1 ! Lb = momentum transfer + LLX=lf+lf ! 2Lb + JLX=LLX+IS2 ! 2 Jb = 2Lb + 2sb IF(JLX.LT.0) GO TO 75 if(i_sym(2)) then if(phasef(lf).gt.0.0) then @@ -579,20 +583,20 @@ c final L loop else temp4=temp2 endif - TEMP4=temp4*SQRT(FLOAT(JLX+1))*float(llx+1) - if(igam2) then + TEMP4=temp4*SQRT(FLOAT(JLX+1))*float(llx+1) ! sqrt(2Jb+1)*(2Lb+1), (d,p) : temp4 = 100 * sqrt((2sa+1)(2l+1)(2Jb+1)) (2Lb+1) + if(igam2) then ! skip by (d,p) temp4=temp4*sqrt(float(ll)/(float(lf)+1.e-6)) - 1 *(z(1)+za(1)*(-fm(1)/fma(1))**lf) + 1 *(z(1)+za(1)*(-fm(1)/fma(1))**lf) endif - LSTOR=lf*IFACT + LSTOR=lf*IFACT ! num. of Lb * (2sa+1) * (2sb+1) LP1=IABS(LL-LTR-1)+1 LP2=MIN0(LL+LTR,LPLUS) LK=0 c initial L loop - DO 60 LP=LP1,LP2,2 + DO 60 LP=LP1,LP2,2 ! loop for La li=lp-1 - LPX=LP+LP-2 - JPX=LPX+IS1 + LPX=LP+LP-2 ! 2La + JPX=LPX+IS1 ! 2Ja IF(JPX.GE.0) then if(i_sym(1)) then if(phasef(li).gt.0.0) then @@ -607,8 +611,18 @@ c initial L loop temp3=temp3*sqrt(float(lp)/(float(li)+1.e-6)) 1 *(z(2)+za(1)*(-fm(2)/fma(2))**li) endif + +c write(6,*) LLX, isf, JLX, LX, ISB, JX, LPX, ISI, JPX + TEMP=temp3*SQRT(FLOAT(LPX+1))*PHASEF((LP-LL-LTR)/2) - 1 *VCC(LLX,LX,LPX,0,0)*WINEJ(LLX,isf,JLX,LX,ISB,JX,LPX,isi,JPX) + 1 *VCC(LLX,LX,LPX,0,0)*WINEJ(LLX,isf,JLX,LX,ISB,JX,LPX,isi,JPX) ! WINEJ = 9j, VCC = CG-coeff +c (d,p) TEMP = temp4 * sqrt(2La+1) (-1)^{La - Lb - l} * CG[{Lb, 0}, {l, 0}, {La, 0}] * NineJ[everything is double] +c VCC = take 2*l, 2*s, 2*j, but compute CG[{l, m}, {s, ms}, {k, m+ms}] +c LLX,isf,JLX,LX,ISB,JX,LPX,isi,JPX +c -> 2Lb,2sb,2Jb,2l,2s ,2j,2La,2sa,2Ja +c TEMP = 100 * sqrt((2sa+1)(2l+1)(2Jb+1)) (2Lb+1) sqrt(2La+1) (-1)^{La - Lb - l} * CG[{Lb, 0}, {l, 0}, {La, 0}] * NineJ[everything is double] +c + INDEX=LK+LL KT=0 c Initial state spins @@ -617,7 +631,7 @@ c Initial state spins c Final state spins MS =-IS(2) DO 55 MS2=1,JS - VCP=VCC(LPX,IS(1),JPX,0,MSP) + VCP=VCC(LPX,IS(1),JPX,0,MSP) ! CG[{La, 0}, {sa, ma}, {Ja, ma}] c DO 50 M=1,MPLUS MK=M+M-1 @@ -625,11 +639,48 @@ c ML2=MSP-MX-MS ML=IABS(ML2/2) IF(ML.LE.lf) then - IND=LSTOR+KT+M + IND=LSTOR+KT+M ! +c +c LSTOR = l * num of total Ja, Jb states +c KT = MPLUS state +c +c write(6,*) LSTOR, KT, M, IND + FACT=VCP*VCC(JLX,JX,JPX,MSP-MX,MX)*VCC(LLX,IS(2),JLX,ML2,MS) 1 *SQRT(YXFCT(lf+ML,lf-ML))*TEMP +c MSP = 2ma +c MX = 2m +c MS = 2mb +c ML2 = 2(ma - m - mb) +c VCC(JLX,JX,JPX,MSP-MX,MX) = CG[{Jb, ma - m}, {j, m}, {Ja, ma}] +c VCC(LLX,IS(2),JLX,ML2,MS) = CG[{Lb, ma - m -mb}, {sb, mb}, {Jb, ma - m}] +c +c FACT = CG[{Jb, ma - m}, {j, m}, {Ja, ma}] * CG[{Lb, ma - m -mb}, {sb, mb}, {Jb, ma - m}] * sqrt((Lb+abs(ma-m-mb)!/(Lb-abs(ma-m-mb)!) * TEMP +c temp = 100 * sqrt((2sa+1)(2l+1)(2Jb+1)) (2Lb+1) sqrt(2La+1) (-1)^{La - Lb - l} * CG[{Lb, 0}, {l, 0}, {La, 0}] * NineJ[everything is double] +c + +c write(6,5432) LPX/2.,JPX/2., LLX/2.,JLX/2., +c 1 MSP, MX, MS, TEMP, VCP, VCC(JLX,JX,JPX,MSP-MX,MX), +c 2 VCC(LLX,IS(2),JLX,ML2,MS), SQRT(YXFCT(lf+ML,lf-ML)), +c 3 FACT, FLL(INDEX) +c 5432 FORMAT(F5.1, F5.1, F5.1, F5.1, +c 1 I4, I4, I4, F15.6, F15.6, F15.6, F15.6, +c 2 F15.6, F15.6, F15.6, F15.6) + D(IND)=D(IND)+FLL(INDEX)*FACT endif + +c IND = (2sa+1)*(2sb+1) * (L+1) + (2sb+1)(ma+sa+1) + (mb+sb+1), loop from (-sa,-sb), (-sa, -sb+1), ...(-sa, sb), (-sa+1,-sb),.... +c D(Lb, ma, mb, m), m > 0 + +c IF(IND.EQ.3) then +c write(6,8765) lf, LPX, JPX, LLX, JLX, MSP, ! l, La, Ja, Lb, Jb, ma, mb, m +c 1 MS, MX, D(IND), FACT*FLL(INDEX), FACT, FLL(INDEX) +c 8765 FORMAT(8I4, 7F15.7) +c ENDIF +c +c + 50 CONTINUE KT = KT+MPLUS MS =MS +2 @@ -645,6 +696,13 @@ c 90 CONTINUE IS1=IS1+2 95 CONTINUE + +c write(6,*) "=============", IND +c DO 888 QQ=1,192 +c write(6,*) QQ, D(QQ) +c 888 CONTINUE + + RETURN END @@ -702,7 +760,7 @@ c initial state average factor for Gamma ray FACTA=sqrt(FACTR) endif c - M2K=(1.0-PHASEF(IS(3)))/2.0 + M2K=(1.0-PHASEF(IS(3)))/2.0 ! PHASEF = (-1)^N NPLUS=(JTR+IS(1)+IS(2))/2+1 MPLUS=JTR/2+1 IFACT = MPLUS*JR*JS @@ -737,6 +795,12 @@ c DO 40 LL=1,LPLUS ML1 =ML1 +1 SUM1 = SUM1+D(IND)*PLM(ML1) + +c sum only for j > m > 0 +c SUM1 = sum( D(Lb, ma, mb, m) P(Lb, ma-m+mb, Cos(theta)), {ma, -sa, sa}, {mb, -sb, sb}) + +c write(6,*) MPLUS, IND, D(IND), ML1, PLM(ML1), SUM1 + C C CALCULATE TOTAL INELASTIC SIGMA C @@ -744,7 +808,7 @@ C L=LL-1 ML = iabs(ML) if(ML.le.L) then - FACT = conjg(D(IND))*D(IND)*YXFCT(L-ML,L+ML)/FLOAT(2*L+1) + FACT = conjg(D(IND))*D(IND)*YXFCT(L-ML,L+ML)/FLOAT(2*L+1) !YXFCT = N!/M! IF(M2 .NE. 0) FACT=FACT*2.0 TotSig=TotSig+FACT endif @@ -757,6 +821,10 @@ C index2 = index2-1 SUM(index2) = SUM1*PHAS2 *FACTA endif + +c write(6,*) IS1, IS2, M, ML, SUM1, FACTA, +c 1 index1, SUM(index1), index2, SUM(index2) + c if(nth.eq.2) write(*,'(a,4i3, 1p4e12.4)') c 1 ' Is2,Is1 M, ML :',is2,is1,M,ML,SUM(Index1),SUM(Index2) KT = KT+MPLUS diff --git a/dwuck4/DWtest2.DAT b/dwuck4/DWtest2.DAT new file mode 100644 index 0000000..52cb5d6 --- /dev/null +++ b/dwuck4/DWtest2.DAT @@ -0,0 +1,16 @@ +10001310500100000 16O(D,P)17O G.S. d5/2 orbital ++181. +00. +01.0 ++15+01+02+05 ++00.10 +12. ++20.00 +02. +01. +16. +08. +01.30 +02. ++01. -88.955 +01.149 +00.675 -02.348 +01.345 +00.603 ++02. +01.394 +00.687 +40.872 +01.394 +00.687 +-04. -14.228 +00.972 +01.011 +01.562 +00.477 ++01.92 +01. +01. +17. +08. +01.42 +01. ++01. -49.544 +01.146 +00.675 -02.061 +01.146 +00.675 ++02. +30.680 +01.302 +00.528 +-04. -21.184 +00.934 +00.590 +00.424 +00.934 +00.590 +-04.14 +01. +00. +16. +08. +01.30 +01. +-01. -01. +01.10 +00.65 +24. ++00. +02. +05. +01. +58. +9 END OF DATA DWUCK4 test cases \ No newline at end of file diff --git a/dwuck4/LGNDR.py b/dwuck4/LGNDR.py new file mode 100755 index 0000000..6bb09f7 --- /dev/null +++ b/dwuck4/LGNDR.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 + +import numpy as np +from scipy.special import lpmv + +def lgndr(mplus, lplus, thet): + """ + Calculates Legendre polynomials Plm + + Parameters: + mplus : int + Number of m's > 0 + lplus : int + Number of l's > 0 + thet : float + Angle in degrees + + Returns: + plm : list + List containing Legendre polynomials + """ + + theta = np.radians(thet) + y = np.cos(theta) + z = np.sin(theta) + plm = np.zeros(459, dtype=np.float64) + + ix = 0 + for m in range(1, mplus + 1): # For MPLUS = 1, LPLUS = 16 + lx = m - 1 # LX = 0 + l2 = 0 # L2 = 0 + p3 = 1.0 # P3 = 1.0 + fl1 = float(lx) # FL1 = 0 + + if lx != 0: + for lt in range(1, lx + 1): + fl1 += 1.0 + p3 *= fl1 * z / 2.0 + + p2 = 0.0 # P2 = 0.0 + fl2 = fl1 + 1.0 # FL2 = 1.0 + fl3 = 1.0 # FL3 = 1.0 + + for lt in range(1, lplus + 1): # Loop Lb + ix1 = ix + lt + + if l2 < lx: + plm[ix1] = 0.0 + else: + if l2 > lx: + p3 = (fl2 * y * p2 - fl1 * p1) / fl3 + fl1 += 1.0 + fl2 += 2.0 + fl3 += 1.0 + plm[ix1] = p3 + print(f'PLM, {lx:3d}, {l2:3d}, {ix1:3d}, {plm[ix1]:15.10f}') + p1, p2 = p2, p3 + + l2 += 1 + + ix += lplus + + return plm + + +plm = lgndr(3, 16, 1) \ No newline at end of file diff --git a/dwuck4/culib8/LGNDR.FOR b/dwuck4/culib8/LGNDR.FOR index 00e8b46..2886832 100644 --- a/dwuck4/culib8/LGNDR.FOR +++ b/dwuck4/culib8/LGNDR.FOR @@ -12,27 +12,31 @@ c IMPLICIT REAL*8(A-H,O-Z) DIMENSION PLM(459) +c +c +c write(6, *) "========== ",MPLUS, LPLUS, THET c THETA=THET /57.295779 Y=COS(THETA) Z=SIN(THETA) IX=0 - DO 100 M=1,MPLUS - LX=M-1 - L2=0 - P3=1.0 - FL1=LX + DO 100 M=1,MPLUS ! for MPLUS = 1, LPLUS = 16 + LX=M-1 ! LX = 0 + L2=0 ! L2 = 0 + P3=1.0 ! P3 = 1.0 + FL1=LX ! FL1 = 0 IF(LX.EQ.0) GO TO 41 DO 40 LT=1,LX FL1=FL1+1.0 P3=P3*FL1*Z/2.0 40 CONTINUE - 41 P2=0.0 - FL2=FL1+1.0 - FL3=1.0 - DO 90 LT=1,LPLUS - IX1=IX+LT - IF(L2-LX)50,70,60 + 41 P2=0.0 ! P2 = 0.0 + FL2=FL1+1.0 !FL2 = 1.0 + FL3=1.0 ! FL3 = 1.0 +c================================= loop Lb + DO 90 LT=1,LPLUS ! loop Lb + IX1=IX+LT + IF(L2-LX)50,70,60 ! if L2 < Lx -> 50; L2 == Lx -> 70; L2 > LX -> 60 50 PLM(IX1)=0.0 GO TO 75 60 P3=(FL2*Y*P2-FL1*P1)/FL3 @@ -40,10 +44,14 @@ c FL2=FL2+2.0 FL3=FL3+1.0 70 PLM(IX1)=P3 + +c write(6, *) 'PLM, ',THETA*57.295779, IX1, PLM(IX1) + P1=P2 P2=P3 75 L2=L2+1 90 CONTINUE +c================================== end of Loop Lb IX=IX+LPLUS 100 CONTINUE RETURN diff --git a/dwuck4/culib8/POLFCT.FOR b/dwuck4/culib8/POLFCT.FOR index fbd3790..1f4875c 100644 --- a/dwuck4/culib8/POLFCT.FOR +++ b/dwuck4/culib8/POLFCT.FOR @@ -125,6 +125,9 @@ c Calculate Tensor polarization = endif 140 continue c + +c write(6,*) mp, m, my, mx, sr(mp, m, my, mx) + Pol(1) =Pol(1) + conjg(sr(mp,m ,my,mx)) * sr(mp,m ,my,mx) Pol(2) =Pol(2) + conjg(sr(mp,m ,my,mx)) * a(mp,m ) Pol(3) =Pol(3) + conjg(sr(mp,m ,my,mx)) * b(mp,m )