rtang/PtolemyGUI
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Merge remote-tracking branch 'origin'

Browse Source
This commit is contained in:
Ryan Tang 2025-02-17 18:23:07 -05:00
parent 8d1b390d17 299cefb9ce
commit d0b359859f
7 changed files with 329 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
dwuck4/ADWCK4.FOR
View File

@ -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

190
dwuck4/BDWCK4.FOR
View File

@ -269,96 +269,163 @@ 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 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
DR2(N)=DR(N)**2/12.0
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)
DO 29 ISS=1,JS
I=I+1
LM(I)=0
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
30 CONTINUE
DO 40 IQ=1,NP
F1(IQ)=0.0
F2(IQ)=0.0
c------------------------------- end of J_state loop
30 CONTINUE ! end of loop N
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
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
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, 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
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) ! 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) ! 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)
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 ! 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)
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 radial grids
LX=1
drrc = 0.1
DO 120 N=1,2
R2=FK(N)*R(N)
R1=R2-DR(N)*FK(N)
@ -383,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)
@ -392,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
@ -407,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

112
dwuck4/CDWCK4.FOR
View File

@ -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

16
dwuck4/DWtest2.DAT Normal file
View File

@ -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

66
dwuck4/LGNDR.py Executable file
View File

@ -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)

30
dwuck4/culib8/LGNDR.FOR
View File

@ -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

3
dwuck4/culib8/POLFCT.FOR
View File

@ -125,6 +125,9 @@ c Calculate Tensor polarization = <Sij(final)>
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 )