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

snapshot for DWUCK4

Browse Source
This commit is contained in:
Ryan@Home 2025-02-16 17:24:22 -05:00
parent 359b7d0def
commit 299cefb9ce
5 changed files with 194 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

112
dwuck4/CDWCK4.FOR
View File

@ -537,12 +537,13 @@ c
igam2=.false. igam2=.false.
endif endif
c c
MPLUS=JX/2+1 MPLUS=JX/2+1 ! j/2 + 1, spin transfer + 1
IFACT=MPLUS*JR*JS 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 J2K=(1.0+PHASEF(NS(2)))/2.0
M2K=JX-MPLUS-MPLUS+2 M2K=JX-MPLUS-MPLUS+2
LX=LTR+LTR LX=LTR+LTR
TEMP2=SQRT(FLOAT((JX+1)*(IS(1)+1)))*flfact 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 IF(FN.EQ.0.0) THEN
c clear amplitude storage unless for coherent sum c clear amplitude storage unless for coherent sum
IND=2*LPLUS*IFACT IND=2*LPLUS*IFACT
@ -550,9 +551,12 @@ c clear amplitude storage unless for coherent sum
D(M)=0.0 D(M)=0.0
10 CONTINUE 10 CONTINUE
ENDIF ENDIF
IS1=-IS(1)
write(6,*) 'Debug', IND
IS1=-IS(1) ! m-state of sa ? step +2
DO 95 I=1,JR DO 95 I=1,JR
IS2=-IS(2) IS2=-IS(2) ! m-state of sb ? step +2
DO 90 J=1,JS DO 90 J=1,JS
IF(NLTR.NE.1) GO TO 14 IF(NLTR.NE.1) GO TO 14
IF(JR*JS.EQ.1) GO TO 15 IF(JR*JS.EQ.1) GO TO 15
@ -565,10 +569,10 @@ C
READ (2)(FLL(INDEX),INDEX=1,INCR) READ (2)(FLL(INDEX),INDEX=1,INCR)
15 continue 15 continue
c final L loop c final L loop
DO 80 LL=1,LPLUS DO 80 LL=1,LPLUS ! loop on Lb
lf=LL-1 lf=LL-1 ! Lb = momentum transfer
LLX=lf+lf LLX=lf+lf ! 2Lb
JLX=LLX+IS2 JLX=LLX+IS2 ! 2 Jb = 2Lb + 2sb
IF(JLX.LT.0) GO TO 75 IF(JLX.LT.0) GO TO 75
if(i_sym(2)) then if(i_sym(2)) then
if(phasef(lf).gt.0.0) then if(phasef(lf).gt.0.0) then
@ -579,20 +583,20 @@ c final L loop
else else
temp4=temp2 temp4=temp2
endif endif
TEMP4=temp4*SQRT(FLOAT(JLX+1))*float(llx+1) 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 if(igam2) then ! skip by (d,p)
temp4=temp4*sqrt(float(ll)/(float(lf)+1.e-6)) 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 endif
LSTOR=lf*IFACT LSTOR=lf*IFACT ! num. of Lb * (2sa+1) * (2sb+1)
LP1=IABS(LL-LTR-1)+1 LP1=IABS(LL-LTR-1)+1
LP2=MIN0(LL+LTR,LPLUS) LP2=MIN0(LL+LTR,LPLUS)
LK=0 LK=0
c initial L loop c initial L loop
DO 60 LP=LP1,LP2,2 DO 60 LP=LP1,LP2,2 ! loop for La
li=lp-1 li=lp-1
LPX=LP+LP-2 LPX=LP+LP-2 ! 2La
JPX=LPX+IS1 JPX=LPX+IS1 ! 2Ja
IF(JPX.GE.0) then IF(JPX.GE.0) then
if(i_sym(1)) then if(i_sym(1)) then
if(phasef(li).gt.0.0) 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)) temp3=temp3*sqrt(float(lp)/(float(li)+1.e-6))
1 *(z(2)+za(1)*(-fm(2)/fma(2))**li) 1 *(z(2)+za(1)*(-fm(2)/fma(2))**li)
endif endif
c write(6,*) LLX, isf, JLX, LX, ISB, JX, LPX, ISI, JPX
TEMP=temp3*SQRT(FLOAT(LPX+1))*PHASEF((LP-LL-LTR)/2) 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 INDEX=LK+LL
KT=0 KT=0
c Initial state spins c Initial state spins
@ -617,7 +631,7 @@ c Initial state spins
c Final state spins c Final state spins
MS =-IS(2) MS =-IS(2)
DO 55 MS2=1,JS 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 c
DO 50 M=1,MPLUS DO 50 M=1,MPLUS
MK=M+M-1 MK=M+M-1
@ -625,11 +639,48 @@ c
ML2=MSP-MX-MS ML2=MSP-MX-MS
ML=IABS(ML2/2) ML=IABS(ML2/2)
IF(ML.LE.lf) then 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) FACT=VCP*VCC(JLX,JX,JPX,MSP-MX,MX)*VCC(LLX,IS(2),JLX,ML2,MS)
1 *SQRT(YXFCT(lf+ML,lf-ML))*TEMP 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 D(IND)=D(IND)+FLL(INDEX)*FACT
endif 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 50 CONTINUE
KT = KT+MPLUS KT = KT+MPLUS
MS =MS +2 MS =MS +2
@ -645,6 +696,13 @@ c
90 CONTINUE 90 CONTINUE
IS1=IS1+2 IS1=IS1+2
95 CONTINUE 95 CONTINUE
c write(6,*) "=============", IND
c DO 888 QQ=1,192
c write(6,*) QQ, D(QQ)
c 888 CONTINUE
RETURN RETURN
END END
@ -702,7 +760,7 @@ c initial state average factor for Gamma ray
FACTA=sqrt(FACTR) FACTA=sqrt(FACTR)
endif endif
c 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 NPLUS=(JTR+IS(1)+IS(2))/2+1
MPLUS=JTR/2+1 MPLUS=JTR/2+1
IFACT = MPLUS*JR*JS IFACT = MPLUS*JR*JS
@ -737,6 +795,12 @@ c
DO 40 LL=1,LPLUS DO 40 LL=1,LPLUS
ML1 =ML1 +1 ML1 =ML1 +1
SUM1 = SUM1+D(IND)*PLM(ML1) 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
C CALCULATE TOTAL INELASTIC SIGMA C CALCULATE TOTAL INELASTIC SIGMA
C C
@ -744,7 +808,7 @@ C
L=LL-1 L=LL-1
ML = iabs(ML) ML = iabs(ML)
if(ML.le.L) then 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 IF(M2 .NE. 0) FACT=FACT*2.0
TotSig=TotSig+FACT TotSig=TotSig+FACT
endif endif
@ -757,6 +821,10 @@ C
index2 = index2-1 index2 = index2-1
SUM(index2) = SUM1*PHAS2 *FACTA SUM(index2) = SUM1*PHAS2 *FACTA
endif 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 if(nth.eq.2) write(*,'(a,4i3, 1p4e12.4)')
c 1 ' Is2,Is1 M, ML :',is2,is1,M,ML,SUM(Index1),SUM(Index2) c 1 ' Is2,Is1 M, ML :',is2,is1,M,ML,SUM(Index1),SUM(Index2)
KT = KT+MPLUS 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) IMPLICIT REAL*8(A-H,O-Z)
DIMENSION PLM(459) DIMENSION PLM(459)
c
c
c write(6, *) "========== ",MPLUS, LPLUS, THET
c c
THETA=THET /57.295779 THETA=THET /57.295779
Y=COS(THETA) Y=COS(THETA)
Z=SIN(THETA) Z=SIN(THETA)
IX=0 IX=0
DO 100 M=1,MPLUS DO 100 M=1,MPLUS ! for MPLUS = 1, LPLUS = 16
LX=M-1 LX=M-1 ! LX = 0
L2=0 L2=0 ! L2 = 0
P3=1.0 P3=1.0 ! P3 = 1.0
FL1=LX FL1=LX ! FL1 = 0
IF(LX.EQ.0) GO TO 41 IF(LX.EQ.0) GO TO 41
DO 40 LT=1,LX DO 40 LT=1,LX
FL1=FL1+1.0 FL1=FL1+1.0
P3=P3*FL1*Z/2.0 P3=P3*FL1*Z/2.0
40 CONTINUE 40 CONTINUE
41 P2=0.0 41 P2=0.0 ! P2 = 0.0
FL2=FL1+1.0 FL2=FL1+1.0 !FL2 = 1.0
FL3=1.0 FL3=1.0 ! FL3 = 1.0
DO 90 LT=1,LPLUS c================================= loop Lb
IX1=IX+LT DO 90 LT=1,LPLUS ! loop Lb
IF(L2-LX)50,70,60 IX1=IX+LT
IF(L2-LX)50,70,60 ! if L2 < Lx -> 50; L2 == Lx -> 70; L2 > LX -> 60
50 PLM(IX1)=0.0 50 PLM(IX1)=0.0
GO TO 75 GO TO 75
60 P3=(FL2*Y*P2-FL1*P1)/FL3 60 P3=(FL2*Y*P2-FL1*P1)/FL3
@ -40,10 +44,14 @@ c
FL2=FL2+2.0 FL2=FL2+2.0
FL3=FL3+1.0 FL3=FL3+1.0
70 PLM(IX1)=P3 70 PLM(IX1)=P3
c write(6, *) 'PLM, ',THETA*57.295779, IX1, PLM(IX1)
P1=P2 P1=P2
P2=P3 P2=P3
75 L2=L2+1 75 L2=L2+1
90 CONTINUE 90 CONTINUE
c================================== end of Loop Lb
IX=IX+LPLUS IX=IX+LPLUS
100 CONTINUE 100 CONTINUE
RETURN RETURN

3
dwuck4/culib8/POLFCT.FOR
View File

@ -125,6 +125,9 @@ c Calculate Tensor polarization = <Sij(final)>
endif endif
140 continue 140 continue
c 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(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(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 ) Pol(3) =Pol(3) + conjg(sr(mp,m ,my,mx)) * b(mp,m )