OSBRACKETS tested
This commit is contained in:
parent
26d5a97d47
commit
ca33915f88
|
|
@ -1,6 +1,9 @@
|
|||
# probably not recommended
|
||||
Project.toml
|
||||
|
||||
# Compiled FORTRAN libraries
|
||||
*.so
|
||||
|
||||
# Misc
|
||||
temp/*
|
||||
scratch/*
|
||||
|
|
|
|||
|
|
@ -0,0 +1,71 @@
|
|||
The programs are written in Fortran-90. The double precision is adopted.
|
||||
|
||||
The following six files make up the package:
|
||||
|
||||
allosbrac.f90
|
||||
testallosbrac.f90
|
||||
allosoutput
|
||||
|
||||
osbrac.f90
|
||||
testosbrac.f90
|
||||
osoutput
|
||||
|
||||
These are two independent groups of files that correspond to the two versions of the present
|
||||
program.
|
||||
|
||||
In one version, the subroutine producing the brackets is named ALLOSBRAC. It is contained in
|
||||
the file allosbrac.f90 along with the routines this subroutine uses.
|
||||
|
||||
In the other version, the subroutine producing the brackets is named OSBRAC. It is contained
|
||||
in the file osbrac.f90 along with the routines this subroutine uses.
|
||||
|
||||
To test the set of programs contained in the file allosbrac.f90, one may say, e.g.,
|
||||
gfortran -O2 testallosbrac.f90 allosbrac.f90. To test the set of programs contained in the
|
||||
file osbrac.f90, one may say, e.g., gfortran -O2 testosbrac.f90 osbrac.f90. The files
|
||||
allosoutput and osoutput are respective output files that contain results of the tests.
|
||||
|
||||
The files testallosbrac.f90 and testosbrac.f90 contain, respectively, only the programs
|
||||
TESTALLOSBRAC and TESTOSBRAC. The tests performed are described in comment lines
|
||||
in these programs and in the accompanying CPC paper.
|
||||
|
||||
The parameters that are set at the moment in the programs TESTALLOSBRAC and TESTOSBRAC
|
||||
are just those for which the results in the files allosoutput and osoutput are listed. With these
|
||||
parameters, the programs run less than a second on a notebook.
|
||||
-------------------------------------------------------------------------------------------------------------------------------------------
|
||||
As said above, the brackets are produced either by the subroutine ALLOSBRAC or by the
|
||||
subroutine OSBRAC. The subroutines contain generous comments.
|
||||
-------------------------------------------------------------------------------------------------------------------------------------------
|
||||
In the case when ALLOSBRAC is employed, all the routines used to calculate the brackets,
|
||||
which are contained in the file allosbrac.f90, are as follows: the subroutine ALLOSBRAC (with
|
||||
comments), the subroutines ARR, FLPHI (with comments), COEFREL (with comments), and the
|
||||
function WIGMOD.
|
||||
|
||||
The subroutine ALLOSBRAC calls for the ARR, COE, FLPHI, and COEFREL subroutines. FLPHI calls
|
||||
for the function WIGMOD. The meaning of the parameters of the ALLOSBRAC subroutine is
|
||||
explained in comments at its beginning.
|
||||
------------------------------------------------------------------------------------------------------------------------------------------
|
||||
In the case when OSBRAC is employed, all the routines used to calculate the brackets,
|
||||
which are contained in the file osbrac.f90, are as follows: the subroutine OSBRAC (with
|
||||
comments), the subroutines ARR, FLPHI (with comments), COEFREL (with comments), and the
|
||||
function WIGMOD.
|
||||
|
||||
The subroutine OSBRAC calls for the ARR, COE, FLPHI, and COEFREL subroutines. FLPHI calls
|
||||
for the function WIGMOD. The meaning of the parameters of the OSBRAC subroutine is explained
|
||||
in comments at its beginning.
|
||||
-------------------------------------------------------------------------------------------------------------------------------------------
|
||||
While it is clear from the comments in the beginnings of ALLOSBRAC and OSBRAC how to
|
||||
implement these subroutines, this is also obvious from examples listed in the above mentioned
|
||||
TESTALLOSBRAC and TESTOSBRAC programs.
|
||||
|
||||
In general, to compile codes that include the present programs one may say e.g.,
|
||||
f95 main.f90 ... allosbrac.f90 ... or, alternatively, f95 main.f90 ... osbrac.f90 ....
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,451 @@
|
|||
SUBROUTINE ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC)
|
||||
|
||||
! COMMENTS ON USE:
|
||||
!---------------------------------------------------------------------------------------------------------------------------------------
|
||||
! THE OSCILLATOR BRACKETS <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\VARPHI ARE
|
||||
! CALCULATED HERE. THIS IS DONE WITH THE HELP OF EQS. (18), (20)-TYPE, AND (22) IN
|
||||
! THE ACCOMPANYING CPC PAPER. THE QUANTITIES L1, L2, L1P, AND L2P IN THE ABOVE
|
||||
! DEFINITION OF THE BRACKETS ARE THE PARTIAL ANGULAR MOMENTA, L IS THE TOTAL
|
||||
! ANGULAR MOMENTUM, AND N1, N2, N1P, N2P ARE THE RADIAL QUANTUM NUMBERS.
|
||||
! IN THE NOTATION LIKE L1P, ETC., "P" SYMBOLIZES "PRIMED" HERE AND BELOW.
|
||||
|
||||
! THE SUBROUTINE RETURNS THE ARRAY OF ALL THE BRACKETS SUCH THAT
|
||||
! L1+L2+2*(N1+N2).LE.NQMAX, LMIN \LE L \LE LMAX, AND THE BRACKETS PERTAIN TO STATES
|
||||
! OF THE SAME PARITY (-1)^(L1+L2) WHICH IS THE PARITY OF NQMAX.
|
||||
|
||||
! THE L1, L2, L1P, AND L2P ORBITAL MOMENTA ARE EXPRESSED IN THE SUBROUTINE IN
|
||||
! TERMS OF THE M, N, MP, AND NP VARIABLES DEFINED AS FOLLOWS
|
||||
! M=(L1+L2-L-NN)/2, N=(L1-L2+L-NN)/2, MP=(L1P+L2P-L-NN)/2, NP=(L1P-L2P+L-NN)/2
|
||||
! WHERE NN EQUALS 0 OR 1 WHEN, RESPECTIVELY, NQMAX-L IS EVEN OR ODD. ONE THEN HAS
|
||||
! L1 = M+N+NN, L2 = M-N+L, L1P = MP+NP+NN, L2P = MP-NP+L.
|
||||
! WHEN L1, L2, L1P AND L2P TAKE ALL THE VALUES ALLOWED AT GIVEN L AND GIVEN PARITY
|
||||
! THE N AND NP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO L-NN AND
|
||||
! THE M AND MP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO (NQMAX-L-NN)/2.
|
||||
! ONE ALSO HAS NQ=2*(M+N1+N2)+L+NN=2*(MP+N1P+N2P)+L+NN.
|
||||
|
||||
! ALL THE PARAMETERS OF THE SUBROUTINE BUT BRAC ARE INPUT ONES. THE MEANING
|
||||
! OF THE NQMAX, LMIN, AND LMAX PARAMETERS IS SEEN FROM THE ABOVE DESCRIPTION.
|
||||
|
||||
! THE BRAC PARAMETER DESIGNATES THE ARRAY OF OUTPUT BRACKETS. IT IS OF THE FORM
|
||||
! BRAC(NP,N1P,MP,N1,N2,N,M,L) . (AS SAID ABOVE, L1=L1(N,M), L2=L2(N,M), L1P=L1P(NP.MP),
|
||||
! AND L2P=L2P(NP,MP). THE QUANTITY N2P IS DETERMINED BY THE EQUALITY
|
||||
! MP+N1P+N2P=M+N1+N2.) THE ORDER OF THE ARGUMENTS OF BRAC CORRESPONDS TO
|
||||
! NESTING OF THE LOOPS AT ITS COMPUTATION.
|
||||
|
||||
! THE ROUTINE PARAMETERS CO AND SI ARE AS FOLLOWS, CO = COS(PHI) AND SI = SIN(PHI).
|
||||
! THESE QUANTITIES DEFINE THE PSEUDO(!!!)ORTHOGONAL TRANSFORMATION
|
||||
! XI1P = CO*XI1+SI*XI2, XI2P = SI*XI1-CO*XI2.
|
||||
|
||||
! (BRACKETS PERTAINING TO THE CASE OF THE ORTHOGONAL TRANSFORMATION
|
||||
! XI1P = CO*XI1+SI*XI2, XI2P = -SI*XI1+CO*XI2 ARE SIMPLY EXPRESSED IN TERMS OF
|
||||
! THOSE CALCULATED IN THE PRESENT PROGRAM, SEE THE ACCOMPANYING CPC PAPER.
|
||||
|
||||
! THE MEANING OF THE ARRAYS, OTHER THAN BRAC, ENTERING THE DIMENSION LIST IS AS
|
||||
! FOLLOWS. B IS A SUBSIDIARY ARRAY TO PERFORM THE RECURSION.
|
||||
|
||||
! FAC(I)=I!, DFAC(I)=(2I+1)!!, AND DEFAC(I)=(2I)!!. RFAC(I)=SQRT((2I)!)/I!
|
||||
! THESE ARRAYS ARE PRODUCED BY THE "ARR" ROUTINE WHICH IS CALLED FROM THE PRESENT
|
||||
! ROUTINE AND WHICH IS CONTAINED IN THE PRESENT FILE. THE SET UPPER BOUNDS OF THESE
|
||||
! ARRAYS ARE SUFFICIENT FOR THE COMPUTATION. THESE ARRAYS ARE USED IN THE PRESENT
|
||||
! ROUTINE AND IN THE "FLPHI" ROUTINE.
|
||||
|
||||
! THE A ARRAY IS A(N,L)=(-1)**N/SQRT(DEFAC(N)*DFAC(N+L)). IT APPEARS BOTH IN THE N1=N2=0
|
||||
! BRACKETS AND IN THE RELATION BETWEEN THE < | > AND [ | ] TYPE BRACKETS.
|
||||
|
||||
! THE BI ARRAY IS BI(M,N)=FAC(N)/[FAC(M)*FAC(N-M)]. IT ENTERS THE 3J SYMBOLS. IT IS USED
|
||||
! IN THE "FLPHI" ROUTINE AND IN THE FUNCTION WIGMOD.
|
||||
|
||||
! THE FL(NP,MP,N) ARRAY REPRESENTS THE QUANTITY
|
||||
! [(2*L1P+1)*(2*L2P+1)]^{1/2}*F_L^VARPHI WHERE F_L^VARPHI IS GIVEN BY EQ. (23) IN THE
|
||||
! ACCOMPANYING CPC PAPER. THIS ARRAY IS PRODUCED IN ADVANCE BY THE
|
||||
! "FLPHI" ROUTINE WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS
|
||||
! CONTAINED IN THE PRESENT FILE. THE VARIABLES MP, NP, AND N ARE DEFINED ABOVE ALONG
|
||||
! WITH THEIR UPPER BOUNDS.
|
||||
|
||||
! THE PSIP(P,Q) AND PSIM(P,Q) ARRAYS REPRESENT THE QUANTITIES (4) IN THE ACCOMPANYING
|
||||
! CPC PAPER. THEY ARE PRODUCED IN ADVANCE BY THE "COEFREC" ROUTINE
|
||||
! WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS CONTAINED IN THE PRESENT
|
||||
! FILE. THEIR ARGUMENTS P AND Q TAKE THE VALUES L1P, L1P+/-1 AND L2P, L2P+/-1,
|
||||
! RESPECTIVELY. THE BOUNDS OF THE PSIP AND PSIM ARRAYS ARE SUCH THAT ALL THE P AND Q
|
||||
! VALUES REQUIRED TO PERFORM THE RECURSION ARE PROVIDED.
|
||||
|
||||
! THE MENTIONED "FLPHI" ROUTINE USES THE FUNCTION WIGMOD WHICH IS ALSO CONTAINED
|
||||
! IN THE PRESENT FILE.
|
||||
!----------------------------------------------------------------------------------------------------------------------------------------
|
||||
DOUBLE PRECISION BRAC,FL,PSIP,PSIM,A,FAC,DFAC,DEFAC,RFAC,BI,CO,SI,T,CO2,SI2,SC,FA,&
|
||||
PLFACT,D2L12P,TN,PA,S,FA1,FA2
|
||||
DIMENSION&
|
||||
BRAC(0:LMAX,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,&
|
||||
0:(NQMAX-LMIN)/2,0:LMAX,0:(NQMAX-LMIN)/2,LMIN:LMAX),&
|
||||
! THE ARRAY IS OF THE FORM BRAC(NP,N1P,MP,N1,N2,N,M,L)
|
||||
FL(0:LMAX,0:(NQMAX-LMIN)/2,0:LMAX),&
|
||||
PSIP((NQMAX+LMAX)/2+1,(NQMAX+LMAX)/2+1),&
|
||||
PSIM((NQMAX+LMAX)/2+1,(NQMAX+LMAX)/2+1),A(0:NQMAX/2+1,0:NQMAX),&
|
||||
FAC(0:2*NQMAX+1),DFAC(0:NQMAX+1),DEFAC(0:NQMAX/2+1),RFAC(0:NQMAX),&
|
||||
BI(0:2*NQMAX+1,0:2*NQMAX+1)
|
||||
IF (NQMAX.GT.84) THEN
|
||||
WRITE(6,*)'IN THE R(*8) COMPUTATION NQMAX SHOULD NOT EXCEED 84'
|
||||
STOP
|
||||
ENDIF
|
||||
IF (CO.EQ.0.D0.OR.SI.EQ.0.D0) THEN
|
||||
WRITE(6,*)'THE PROGRAM IS OF NO USE AT ZERO COS(PHI) OR SIN(PHI) VALUES'
|
||||
WRITE(6,*)'COS(PHI)=',CO,' SIN(PHI)=',SI
|
||||
STOP
|
||||
ENDIF
|
||||
IF (NQMAX.LT.LMAX) THEN
|
||||
WRITE(6,*)'NQMAX=',NQMAX,' L=',LMAX
|
||||
WRITE(6,*)'L SHOULD NOT EXCEED NQMAX'
|
||||
STOP
|
||||
ENDIF
|
||||
T=SI/CO
|
||||
CO2=CO**2
|
||||
SI2=SI**2
|
||||
SC=SI*CO
|
||||
CALL ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
|
||||
CALL COE(2*NQMAX+1,FAC,BI)
|
||||
! THE A_{NL} ARRAY:
|
||||
DO NI=0,NQMAX/2+1
|
||||
NA=NI-2*(NI/2)
|
||||
K=1
|
||||
IF(NA.EQ.1)K=-1
|
||||
DO LI=0,NQMAX-NI
|
||||
A(NI,LI)=K/SQRT(DEFAC(NI)*DFAC(NI+LI))
|
||||
ENDDO
|
||||
ENDDO
|
||||
! BRACKETS [N_1'L_1'N_2'L_2'|0L_10L_2]_L^\PHI TO START THE RECURSION, EQ. (22) IN
|
||||
! THE ACCOMPANYING CPC PAPER.
|
||||
! THESE BRACKETS ARE REPRESENTED AS BRAC(NP,N1P,MP,0,0,N,M,L).
|
||||
! SUBSIDIARY QUANTITIES:
|
||||
DO L=LMIN,LMAX
|
||||
CALL FLPHI(NQMAX,LMAX,LMIN,L,T,BI,RFAC,FL)
|
||||
! THIS SUBROUTINE PRODUCES THE ARRAY FL USED BELOW.
|
||||
NQML=NQMAX-L
|
||||
NQMLD2=NQML/2
|
||||
NN=NQML-2*(NQMLD2)
|
||||
CALL COEFREC((NQMAX+LMAX)/2+1,L,NN,NQMAX,PSIP,PSIM)
|
||||
NMAX=L-NN
|
||||
MMAX=(NQML-NN)/2
|
||||
DO M=0,MMAX
|
||||
DO N=0,NMAX
|
||||
L1=M+N+NN
|
||||
L2=M-N+L
|
||||
L1L2=L1+L2
|
||||
L1L2ML=L1L2-L
|
||||
FA=SQRT((2.D0*L1+1)*(2*L2+1))*CO**L1L2
|
||||
PLFACT=SQRT(FAC(L1L2+L+1)*FAC(L1L2ML))
|
||||
! COMPUTATION OF THE N1=N2=0 BRACKETS. THEY ARE REPRESENTED AS
|
||||
! BRAC(NP,N1P,MP,0,0,N,M,L).
|
||||
MPMAX=(L1L2ML-NN)/2
|
||||
DO MP=0,MPMAX
|
||||
L12P=2*MP+L+NN
|
||||
! THIS IS BY DEFINITION, L12P=L1P+L2P
|
||||
D2L12P=1.D0/2.D0**L12P
|
||||
N12P=MPMAX-MP
|
||||
! INDEED, 2*N12P=2*(N1P+N2P)=NQ-L12P. FOR THE BRACKETS WE COMPUTE NOW
|
||||
! WE HAVE NQ=L1L2. THEN 2*N12P=L1L2-L12P. SINCE L1L2=2*MPMAX-L-NN AND
|
||||
! L12P=2*MP+L+NN, ONE GETS N12P=MPMAX-MP.
|
||||
DO N1P=0,N12P
|
||||
N1PA=N1P-2*(N1P/2)
|
||||
IF (N1PA.EQ.0) THEN
|
||||
K=1
|
||||
ELSE
|
||||
K=-1
|
||||
ENDIF
|
||||
! K=(-1)**N1P
|
||||
N2P=N12P-N1P
|
||||
TN=T**N12P*K
|
||||
DO NP=0,NMAX
|
||||
L1P=MP+NP+NN
|
||||
L2P=L12P-L1P
|
||||
! BY DEFINITION ONE ALSO HAS: L2P=L2P(MP,NP)=MP-NP+L
|
||||
PA=A(N1P,L1P)*A(N2P,L2P)
|
||||
BRAC(NP,N1P,MP,0,0,N,M,L)=FL(NP,MP,N)*FA*TN*D2L12P&
|
||||
*PLFACT*PA*PA
|
||||
ENDDO ! NP
|
||||
ENDDO ! N1P
|
||||
ENDDO ! MP
|
||||
! RECURSION TO OBTAIN THE GENERAL FORM BRACKETS
|
||||
! [N_1'L_1'N_2'L_2'|N_1L_1N_2L_2]_L^\PHI.
|
||||
N12MAX=(NQMAX-L1L2)/2
|
||||
MPMAX0=MPMAX
|
||||
! THE N1=0, N2-1-->N2 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER WITH THE
|
||||
! MODIFICATION POINTED OUT THERE:
|
||||
DO N2=1,N12MAX
|
||||
MPMAX=MPMAX+1
|
||||
! THIS MPMAX VALUE IS (NQ-L-NN)/2 WHERE NQ=L1+L2+2*N2.
|
||||
DO MP=0,MPMAX
|
||||
N12P=MPMAX-MP
|
||||
DO N1P=0,N12P
|
||||
! MP=(L1P+L2P-L-NN)/2. THEREFORE, N12P=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
|
||||
DO NP=0,NMAX
|
||||
L1P=MP+NP+NN
|
||||
L2P=MP-NP+L
|
||||
S=0.D0
|
||||
! BELOW THE RESTRICTIONS ARE IMPOSED ON THE BRAC ARRAY ENTERING THE RIGHT-HAND
|
||||
! SIDE OF THE RECURRENCE FORMULAE. THIS ARRAY IS OF THE FORM BRAC(K1,K2,K3,...) WHERE
|
||||
! K1=K1(NP), K2=N2(N1P), AND K3=K3(MP). NAMELY, K1=NP, OR NP+1, OR NP-1;
|
||||
! K2=N1P OR N1P-1; K3=MP, OR MP-1, OR MP+1. (THE COMBINATION K2=N1P AND K3=MP+1
|
||||
! DOES NOT ARISE IN THE RECURRENCE FORMULAE.) THE QUANTITIES K1, K2, AND K3
|
||||
! REPRESENT, RESPECTIELY, NP, N1P, AND MP VALUES AT THE PRECEDING STAGE OF THE
|
||||
! RECURSION. THE RESTRICTIONS ENSURE THAT K1, K2, AND K3 RANGE WITHIN THE LIMITS
|
||||
! PERTAINING TO THE BRAC ARRAY OBTAINED AT THAT PRECEDING STAGE OF THE RECURSION.
|
||||
|
||||
! THUS K1, K2, AND K3 SHOULD BE NON-NEGATIVE. THEREFORE, WHEN K1=NP-1 THE VALUE
|
||||
! NP=0 IS TO BE EXCLUDED, WHEN K2=N1P-1 THE VALUE N1P=0 IS TO BE EXCLUDED, AND
|
||||
! WHEN K3=MP-1 THE VALUE MP=0 IS TO BE EXCLUDED.
|
||||
|
||||
! FURTHERMORE, IT SHOULD BE K1 \LE NMAX AND AT THE SAME TIME ONE HAS NP \LE NMAX.
|
||||
! THEREFORE, WHEN K1=NP+1 THE VALUE NP=NPMAX IS TO BE EXCLUDED.
|
||||
|
||||
! IN ADDITION, IT SHOULD BE K2+K3 \LE MPMAX-1 WHILE ONE HAS N1P+MP \LE MPMAX. WHEN
|
||||
! K2=N1P AND K3=MP-1, OR K2=N1P-1 AND K3=MP-1, OR K2=N1P-1 AND K3=MP ONE
|
||||
! AUTOMATICALLY HAS N1P+MP \LE MPMAX. THEREFORE, IN THESE CASES THE CONDITION
|
||||
! K2+K3 \LE MPMAX-1 DOES NOT CREATE ANY RESTRICTIONS.
|
||||
! BUT WHEN K2=N1P AND K3=MP, OR K2=N1P-1 AND K3=MP+1 THE CASE N1P+MP=MPMAX
|
||||
! IS TO BE FORBIDDEN.
|
||||
|
||||
! IN THE RECURRENCE FORMULAE BELOW THE DESCRIBED RESTRICTIONS ARE
|
||||
! IMPOSED.
|
||||
IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,0,N2-1,N,M,L)*PSIP(L1P,L2P)
|
||||
IF (N1P.NE.0.AND.N1P.NE.N12P) S=S+&
|
||||
BRAC(NP,N1P-1,MP+1,0,N2-1,N,M,L)*PSIP(L1P+1,L2P+1)
|
||||
! RECALL THAT N12P=MPMAX-MP.
|
||||
IF (NP.NE.NMAX.AND.N1P.NE.0) S=S-&
|
||||
BRAC(NP+1,N1P-1,MP,0,N2-1,N,M,L)*PSIM(L1P+1,L2P)
|
||||
IF (NP.NE.0.AND.N1P.NE.N12P) S=S-&
|
||||
BRAC(NP-1,N1P,MP,0,N2-1,N,M,L)*PSIM(L1P,L2P+1)
|
||||
S=S*SC
|
||||
IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,0,N2-1,N,M,L)*SI2
|
||||
IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,0,N2-1,N,M,L)*CO2
|
||||
BRAC(NP,N1P,MP,0,N2,N,M,L)=S
|
||||
ENDDO ! NP
|
||||
ENDDO ! N1P
|
||||
ENDDO ! MP
|
||||
ENDDO ! N2
|
||||
! N1-1-->N1 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER:
|
||||
DO N2=0,N12MAX
|
||||
MPMAX=MPMAX0+N2
|
||||
! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2=MPMAX0+N2 SINCE NQ=L1+L2+2*N2.
|
||||
! THE RECURSION:
|
||||
DO N1=1,N12MAX-N2
|
||||
MPMAX=MPMAX+1
|
||||
! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2 AND NQ=L1+L2+2*N2+2*N1.
|
||||
DO MP=0,MPMAX
|
||||
N12P=MPMAX-MP
|
||||
DO N1P=0,N12P
|
||||
! MP=(L1P+L2P-L-NN)/2. THEREFORE, N1PMAX=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
|
||||
DO NP=0,NMAX
|
||||
L1P=MP+NP+NN
|
||||
L2P=MP-NP+L
|
||||
S=0.D0
|
||||
IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,N1-1,N2,N,M,L)*PSIP(L1P,L2P)
|
||||
IF (N1P.NE.0.AND.N1P.NE.N12P) S=S+&
|
||||
BRAC(NP,N1P-1,MP+1,N1-1,N2,N,M,L)*PSIP(L1P+1,L2P+1)
|
||||
IF (NP.NE.NMAX.AND.N1P.NE.0) S=S-&
|
||||
BRAC(NP+1,N1P-1,MP,N1-1,N2,N,M,L)*PSIM(L1P+1,L2P)
|
||||
IF (NP.NE.0.AND.N1P.NE.N12P) S=S-&
|
||||
BRAC(NP-1,N1P,MP,N1-1,N2,N,M,L)*PSIM(L1P,L2P+1)
|
||||
! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
|
||||
S=-S*SC
|
||||
IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,N1-1,N2,N,M,L)*CO2
|
||||
IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,N1-1,N2,N,M,L)*SI2
|
||||
! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
|
||||
BRAC(NP,N1P,MP,N1,N2,N,M,L)=S
|
||||
ENDDO ! NP
|
||||
ENDDO ! N1P
|
||||
ENDDO ! MP
|
||||
ENDDO ! N1
|
||||
ENDDO ! N2
|
||||
! RENORMALIZATION OF THE BRACKETS, EQ. (18) IN THE ACCOMPANYING CPC PAPER:
|
||||
DO N2=0,N12MAX
|
||||
MPMAX1=MPMAX0+N2
|
||||
DO N1=0,N12MAX-N2
|
||||
MPMAX=MPMAX1+N1
|
||||
FA1=A(N1,L1)*A(N2,L2)
|
||||
DO MP=0,MPMAX
|
||||
N12P=MPMAX-MP
|
||||
DO N1P=0,MPMAX-MP
|
||||
N2P=N12P-N1P
|
||||
DO NP=0,NMAX
|
||||
L1P=MP+NP+NN
|
||||
L2P=MP-NP+L
|
||||
FA2=A(N1P,L1P)*A(N2P,L2P)
|
||||
BRAC(NP,N1P,MP,N1,N2,N,M,L)=BRAC(NP,N1P,MP,N1,N2,N,M,L) &
|
||||
*FA1/FA2
|
||||
ENDDO ! NP
|
||||
ENDDO ! N1P
|
||||
ENDDO ! MP
|
||||
ENDDO ! N1
|
||||
ENDDO ! N2
|
||||
ENDDO ! N
|
||||
ENDDO ! M
|
||||
ENDDO ! L
|
||||
RETURN
|
||||
END
|
||||
|
||||
SUBROUTINE ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
|
||||
! FAC(I), DFAC(I),DEFAC(I), AND RFAC(I) ARE, RESPECTIVELY, THE QUANTITIES
|
||||
! I!, (2I+1)!!, (2I)!!, AND SQRT((2*I)!)/I!
|
||||
DOUBLE PRECISION FAC,DFAC,DEFAC,RFAC
|
||||
DIMENSION FAC(0:2*NQMAX+1),DFAC(0:NQMAX+1),DEFAC(0:NQMAX/2+1),RFAC(0:NQMAX)
|
||||
FAC(0)=1.D0
|
||||
DFAC(0)=1.D0
|
||||
DEFAC(0)=1.D0
|
||||
RFAC(0)=1.D0
|
||||
DO I=1,2*NQMAX+1
|
||||
FAC(I)=FAC(I-1)*I
|
||||
ENDDO
|
||||
DO I=1,NQMAX+1
|
||||
DFAC(I)=DFAC(I-1)*(2*I+1)
|
||||
ENDDO
|
||||
DO I=1,NQMAX/2+1
|
||||
DEFAC(I)=DEFAC(I-1)*2*I
|
||||
ENDDO
|
||||
DO I=1,NQMAX
|
||||
RFAC(I)=RFAC(I-1)*2*SQRT(1-0.5D0/I)
|
||||
ENDDO
|
||||
RETURN
|
||||
END
|
||||
|
||||
|
||||
SUBROUTINE FLPHI(NQMAX,LMAX,LMIN,L,T,BI,RFAC,FL)
|
||||
! PROVIDES THE QUANTITY SQRT((2L1P+1)*(2L2P+1))*F_L^\VARPHI, SEE EQ. (23) IN THE
|
||||
! ACCOMPANYING PAPER, IN THE FORM OF THE FL ARRAY.
|
||||
! USES THE FUNCTION WIGMOD.
|
||||
! THIS FLPHI DIFFERS FROM THAT IN THE OTHER FILE osbrac.f90.
|
||||
DOUBLE PRECISION BI,RFAC,FL,T,T2,SQP,F,WIGMOD,z
|
||||
DIMENSION&
|
||||
FL(0:LMAX,0:(NQMAX-LMIN)/2,0:LMAX),BI(0:2*NQMAX+1,0:2*NQMAX+1),RFAC(0:NQMAX)
|
||||
! THE OUTPUT IS FL(NP,MP,N) WHERE MP=(L1P+L2P-L-NN)/2, L1P+L2P=L1T+L2T,
|
||||
! NP=(L1P-L2P+L-NN)/2, AND N=(L1-L2+L-NN)/2, L1-L2=L1T-L2T.
|
||||
T2=T*T
|
||||
NQMAL=NQMAX-L
|
||||
MPMAX=NQMAL/2
|
||||
NN=NQMAL-2*MPMAX
|
||||
LMNN=L-NN
|
||||
LL3=2*L
|
||||
DO N=0,LMNN
|
||||
L1ML2=2*N-LMNN
|
||||
DO MP=0,MPMAX
|
||||
L1PL2P=2*MP+L+NN
|
||||
DO NP=0,LMNN
|
||||
L1PML2P=2*NP-LMNN
|
||||
L1P=(L1PL2P+L1PML2P)/2
|
||||
L2P=(L1PL2P-L1PML2P)/2
|
||||
M3=L2P-L1P
|
||||
SQP=SQRT((2*L1P+1.D0)*(2*L2P+1))
|
||||
L1T=(L1PL2P+L1ML2)/2
|
||||
L2T=(L1PL2P-L1ML2)/2
|
||||
N3=L1T+L2T-L
|
||||
IMIN=ABS(L1T-L1P)
|
||||
IMAX=MIN(L1P+L1T,L2P+L2T)
|
||||
NAL1=(L1T+L1P-IMIN)/2
|
||||
NAL2=(L1T-L1P+IMIN)/2
|
||||
NAL3=(L2T-L2P+IMIN)/2
|
||||
NAL4=(L2T+L2P-IMIN)/2
|
||||
NAL4P=NAL4-2*(NAL4/2)
|
||||
F=0.D0
|
||||
z=t**imin
|
||||
if(nal4p.ne.0)z=-z
|
||||
DO I=IMIN,IMAX,2
|
||||
F=F+RFAC(NAL1)*RFAC(NAL2)*RFAC(NAL3)*RFAC(NAL4)*&
|
||||
WIGMOD(L1T,L2T,L,L1P-I,I-L2P,BI,2*NQMAX+1)*z
|
||||
NAL1=NAL1-1
|
||||
NAL2=NAL2+1
|
||||
NAL3=NAL3+1
|
||||
NAL4=NAL4-1
|
||||
z=-z*t2
|
||||
ENDDO
|
||||
FL(NP,MP,N)=F*SQRT((2*L1P+1)*(2*L2P+1)&
|
||||
*BI(L1T+L-L2T,2*L1T)*BI(N3,2*L2T)/((LL3+1)*BI(N3,L1T+L2T+L+1)*BI(L+M3,LL3)))
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
RETURN
|
||||
END
|
||||
|
||||
FUNCTION WIGMOD(L1,L2,L3,M1,M2,BI,NMAX)
|
||||
! THE 3J SYMBOL IN TERMS OF BINOMIAL COEFFICIENTS WITHOUT THE FACTOR
|
||||
! SQRT(F) WHERE F IS AS FOLLOWS,
|
||||
! F=BI(L1+L3-L2,2*L1)*BI(N3,2*L2)/((LL3+1)*BI(N3,L1+L2+L3+1)*BI(L3+M3,LL3))
|
||||
! WITH N3=L1+L2-L3, LL3=2*L3, AND M3=-M1-M2.
|
||||
DOUBLE PRECISION WIGMOD,BI,S
|
||||
DIMENSION BI(0:NMAX,0:NMAX)
|
||||
M3=-M1-M2
|
||||
N3=L1+L2-L3
|
||||
LM1=L1-M1
|
||||
LP2=L2+M2
|
||||
KMIN=MAX(0,L1-L3+M2,L2-L3-M1)
|
||||
KMAX=MIN(N3,LM1,LP2)
|
||||
S=0.d0
|
||||
NPH=1
|
||||
DO K=KMIN,KMAX
|
||||
S=S+NPH*BI(K,N3)*BI(LM1-K,L1+L3-L2)*BI(LP2-K,L2+L3-L1)
|
||||
NPH=-NPH
|
||||
ENDDO
|
||||
WIGMOD=S/SQRT(BI(LM1,2*L1)*BI(LP2,2*L2))
|
||||
NY=KMIN+L1-L2-M3
|
||||
NYP=NY-2*(NY/2)
|
||||
IF(NYP.NE.0)WIGMOD=-WIGMOD
|
||||
RETURN
|
||||
END
|
||||
|
||||
SUBROUTINE COE(NMAX,FAC,BI)
|
||||
DOUBLE PRECISION FAC,BI
|
||||
DIMENSION FAC(0:NMAX),BI(0:NMAX,0:NMAX)
|
||||
DO N=0,NMAX
|
||||
DO M=0,N/2
|
||||
BI(M,N)=FAC(N)/(FAC(M)*FAC(N-M))
|
||||
BI(N-M,N)=BI(M,N)
|
||||
ENDDO
|
||||
ENDDO
|
||||
RETURN
|
||||
END
|
||||
|
||||
SUBROUTINE COEFREC(MAXCOE,L,NN,NQMAX,PSIP,PSIM)
|
||||
! CALCULATES THE COEFFICIENTS OF THE RECURSION FORMULA, EQ. (21) IN THE ACCOMPANYING
|
||||
! CPC PAPER
|
||||
! MAXCOE EQUALS INT((NQMAX+LMAX)/2)+1 AT CALLS OF THIS ROUTINE.
|
||||
DOUBLE PRECISION PSIP,PSIM
|
||||
DIMENSION PSIP(MAXCOE,MAXCOE),PSIM(MAXCOE,MAXCOE)
|
||||
LP=L+1
|
||||
LM=L-1
|
||||
DO M1=1,(NQMAX+L)/2+1
|
||||
DO M2=1,M1
|
||||
M1P2=M1+M2
|
||||
NA=M1P2+L
|
||||
NNA=NA-2*(NA/2)
|
||||
M1M2=M1-M2
|
||||
! CALCULATION OF PSIP. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P AND L2P OR
|
||||
! L1P+1 AND L2P+1.)
|
||||
IF (NNA.EQ.NN.AND.M1P2.GE.L.AND.M1P2.LE.NQMAX.AND.&
|
||||
ABS(M1M2).LE.L) THEN
|
||||
IF (M1P2.GT.LP) THEN
|
||||
M1P2L=M1P2+L
|
||||
M1P2ML=M1P2-L
|
||||
PSIP(M1,M2)=SQRT(M1P2L*(M1P2L+1.D0)*M1P2ML*(M1P2ML-1)/((4*M1*M1-1)&
|
||||
*(4*M2*M2-1)))
|
||||
ELSE
|
||||
PSIP(M1,M2)=0.D0
|
||||
ENDIF
|
||||
PSIP(M2,M1)=PSIP(M1,M2)
|
||||
ENDIF
|
||||
! CALCULATION OF PSIM. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P+1 AND L2P
|
||||
! OR L2P+1 AND L1P.)
|
||||
IF (NNA.NE.NN.AND.M1P2.GE.LP.AND.M1P2.LE.NQMAX&
|
||||
.AND.M1M2.LE.LP.AND.M1M2.GE.-LM) THEN
|
||||
IF (M1M2.LT.L) THEN
|
||||
M1M2L=M1M2+L
|
||||
M1M2ML=M1M2-L
|
||||
PSIM(M1,M2)=SQRT(M1M2L*(M1M2L+1.D0)*M1M2ML*(M1M2ML-1)/&
|
||||
((4*M1*M1-1)*(4*M2*M2-1)))
|
||||
ELSE
|
||||
PSIM(M1,M2)=0.D0
|
||||
ENDIF
|
||||
PSIM(M2,M1)=PSIM(M1,M2)
|
||||
ENDIF
|
||||
ENDDO
|
||||
ENDDO
|
||||
RETURN
|
||||
END
|
||||
|
||||
|
|
@ -0,0 +1,467 @@
|
|||
SUBROUTINE OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
|
||||
|
||||
! COMMENTS ON USE:
|
||||
!---------------------------------------------------------------------------------------------------------------------------------------
|
||||
! THE OSCILLATOR BRACKETS <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\VARPHI ARE
|
||||
! CALCULATED HERE. THIS IS DONE WITH THE HELP OF EQS. (18), (20)-TYPE, AND (22) IN
|
||||
! THE ACCOMPANYING CPC PAPER. THE QUANTITIES L1, L2, L1P, AND L2P IN THE ABOVE
|
||||
! DEFINITION OF THE BRACKETS ARE THE PARTIAL ANGULAR MOMENTA, L IS THE TOTAL
|
||||
! ANGULAR MOMENTUM, AND N1, N2, N1P, N2P ARE THE RADIAL QUANTUM NUMBERS.
|
||||
! IN THE NOTATION LIKE L1P, ETC., "P" SYMBOLIZES "PRIMED" HERE AND BELOW.
|
||||
|
||||
! N,M,L, AND NQMAX ARE INPUT PARAMETERS OF THE SUBROUTINE, AND BRAC IS THE
|
||||
! OUTPUT PARAMETER REPRESENTING THE ARRAY OF CALCULATED BRACKETS. THE PARAMETERS
|
||||
! N AND M REPRESENT THE L1 AND L2 QUANTUM NUMBERS, SEE BELOW, SO THAT L1=L1(N,M)
|
||||
! AND L2=L2(N,M).
|
||||
|
||||
! LET US USE THE NOTATION NQ=L1+L2+2*(N1+N2). THE ARRAY BRAC CONTAINS ALL BRACKETS
|
||||
! EXISTING AT GIVEN INPUT L, L1, AND L2 VALUES AND HAVING N1 AND N2 SUCH THAT
|
||||
! NQ.LE.NQMAX. IT IS CLEAR THAT NQ AND NQMAX PERTAINING TO ALL THESE BRACKETS ARE OF
|
||||
! THE SAME PARITY EQUAL TO (-1)^(L1+L2).
|
||||
|
||||
! IN OTHER WORDS, THE COMPUTED BRACKETS PERTAIN TO STATES WITH GIVEN N, M, AND L
|
||||
! VALUES SUCH THAT NQ.LE.NQMAX. AUTOMATICALLY, THESE STATES ARE OF THE SAME PARITY
|
||||
! AS THE PARITY OF NQMAX. (AND THIS PARITY EQUALS (-1)^(L1+L2).)
|
||||
|
||||
! THE L1, L2, L1P, AND L2P ORBITAL MOMENTA ARE EXPRESSED IN TERMS OF THE N, M, NP, AND
|
||||
! MP VARIABLES AS FOLLOWS
|
||||
! N=(L1-L2+L-NN)/2, M=(L1+L2-L-NN)/2, NP=(L1P-L2P+L-NN)/2, MP=(L1P+L2P-L-NN)/2
|
||||
! WHERE NN EQUALS 0 OR 1 WHEN, RESPECTIVELY, NQMAX-L IS EVEN OR ODD. ONE THEN HAS
|
||||
! L1 = M+N+NN, L2 = M-N+L, L1P = MP+NP+NN, L2P = MP-NP+L.
|
||||
! WHEN L1, L2, L1P AND L2P TAKE ALL THE VALUES ALLOWED AT GIVEN L AND GIVEN PARITY
|
||||
! THE N AND NP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO L-NN AND
|
||||
! THE M AND MP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO (NQMAX-L-NN)/2.
|
||||
! ONE ALSO HAS NQ=2*(M+N1+N2)+L+NN=2*(MP+N1P+N2P)+L+NN.
|
||||
|
||||
! THE SAID ABOVE DETERMINES THE MEANING OF THE N,M, AND NQMAX PARAMETERS OF
|
||||
! THE SUBROUTINE.
|
||||
|
||||
! THE BRAC ARRAY IS OF THE FORM BRAC(NP,N1P,MP,N1,N2). AS SAID ABOVE, NP AND MP
|
||||
! RERESENT L1P AND L2P. THE MISSING N2P VALUE IS DETERMNED BY THE EQUALITY
|
||||
! MP+N1P+N2P=M+N1+N2 THAT IS ANOTHER FORM OF THE RELATION
|
||||
! 2(N1+N2)+L1+L2=2(N1P+N2P)+L1P+L2P.
|
||||
|
||||
! THE ROUTINE PARAMETERS CO AND SI ARE AS FOLLOWS, CO = COS(PHI) AND SI = SIN(PHI).
|
||||
! THESE QUANTITIES DEFINE THE PSEUDO(!)ORTHOGONAL TRANSFORMATION
|
||||
! XI1P = CO*XI1+SI*XI2, XI2P = SI*XI1-CO*XI2.
|
||||
|
||||
! (BRACKETS PERTAINING TO THE CASE OF THE ORTHOGONAL TRANSFORMATION
|
||||
! XI1P = CO*XI1+SI*XI2, XI2P = -SI*XI1+CO*XI2 ARE SIMPLY EXPRESSED IN TERMS OF
|
||||
! THOSE CALCULATED IN THE PRESENT PROGRAM, SEE THE ACCOMPANYING CPC PAPER.
|
||||
|
||||
! THE MEANING OF THE FIRSTCALL PARAMETER IS EXPLAINED BELOW.
|
||||
|
||||
! THE MEANING OF THE ARRAYS, OTHER THAN BRAC, ENTERING THE DIMENSION LIST IS AS
|
||||
! FOLLOWS. B IS A SUBSIDIARY ARRAY TO PERFORM THE RECURSION.
|
||||
|
||||
! FAC(I)=I!, DFAC(I)=(2I+1)!!, AND DEFAC(I)=(2I)!!, AND RFAC(I)=SQRT((2I)!)/I!
|
||||
! IS A SUBSIDIARY ARRAY TO CALCULATE THE FL ARRAY, SEE BELOW. THESE ARRAYS ARE
|
||||
! PRODUCED BY THE "ARR" ROUTINE WHICH IS CALLED FROM THE PRESENT ROUTINE AND
|
||||
! WHICH IS CONTAINED IN THE PRESENT FILE. THE SET UPPER BOUNDS OF THESE
|
||||
! ARRAYS ARE SUFFICIENT FOR THE COMPUTATION. THESE ARRAYS ARE USED IN THE PRESENT
|
||||
! ROUTINE AND IN THE "FLPHI" ROUTINE.
|
||||
|
||||
! THE A ARRAY IS A(N,L)=(-1)**N/SQRT(DEFAC(N)*DFAC(N+L)). IT APPEARS BOTH
|
||||
! IN THE N1=N2=0 BRACKETS AND IN THE RELATION BETWEEN THE < | > AND [ | ] TYPE
|
||||
! BRACKETS.
|
||||
|
||||
! THE BI ARRAY IS BI(M,N)=FAC(N)/[FAC(M)*FAC(N-M)]. IT ENTERS THE 3J SYMBOLS. IT IS USED
|
||||
! IN THE "FLPHI" ROUTINE AND IN THE FUNCTION WIGMOD.
|
||||
|
||||
! THE FL(NP,MP,N) ARRAY REPRESENTS THE QUANTITY
|
||||
! [(2*L1P+1)*(2*L2P+1)]^{1/2}*F_L^VARPHI WHERE F_L^VARPHI IS GIVEN BY EQ. (23) IN THE
|
||||
! ACCOMPANYING CPC PAPER. THIS ARRAY IS PRODUCED IN ADVANCE BY THE
|
||||
! "FLPHI" ROUTINE WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS
|
||||
! CONTAINED IN THE PRESENT FILE. THE VARIABLES MP, NP, AND N ARE DEFINED ABOVE ALONG
|
||||
! WITH THEIR UPPER BOUNDS.
|
||||
|
||||
! THE PSIP(P,Q) AND PSIM(P,Q) ARRAYS REPRESENT THE QUANTITIES (4) IN THE ACCOMPANYING
|
||||
! CPC PAPER. THEY ARE PRODUCED IN ADVANCE BY THE "COEFREC" ROUTINE
|
||||
! WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS CONTAINED IN THE PRESENT
|
||||
! FILE. THEIR ARGUMENTS P AND Q TAKE THE VALUES L1P, L1P+/-1 AND L2P, L2P+/-1,
|
||||
! RESPECTIVELY. THE BOUNDS OF THE PSIP AND PSIM ARRAYS ARE SUCH THAT ALL THE P AND Q
|
||||
! VALUES REQUIRED TO PERFORM THE RECURSION ARE PROVIDED.
|
||||
|
||||
! THE MENTIONED "FLPHI" ROUTINE USES THE FUNCTION WIGMOD WHICH IS ALSO CONTAINED
|
||||
! IN THE PRESENT FILE.
|
||||
|
||||
! THE PARAMETER FIRSTCALL IS A LOGICAL VARIABLE WHICH ENSURES THAT FLPHI
|
||||
! IS CALLED ONLY ONCE AT A GIVEN L VALUE PROVIDING THAT ALL THE COMPUTATIONS ARE
|
||||
! PERFORMED WITH THE SAME NQMAX. UNNECESSARY REPEATED CALLS OF FLPHI ARE
|
||||
! SUPPRESSED VIA SAYING FIRSTCALL=.FALSE. IN A PROPER PLACE. IF THE CALLING PROGRAM
|
||||
! INCLUDES NESTED LOOPS OVER N-TYPE VARIABLES FROM N=0 UP TO N=NMAX AND THESE
|
||||
! LOOPS CONTAIN CALLS OF OSBRAC THEN FIRSTCALL=.FALSE. IS TO BE SAID AFTER THE MOST
|
||||
! INNER OF THESE LOOPS WITH THIS PURPOSE. AND AFTER A CHANGE OF L ONE SHOULD
|
||||
! CALL FLPHI ANEW VIA SAYING FIRSTCALL=.TRUE. (SEE ALSO THE APPENDED PROGRAM
|
||||
! TESTOSBRAC AS AN EXAMPLE.)
|
||||
!----------------------------------------------------------------------------------------------------------------------------------------
|
||||
DOUBLE PRECISION BRAC,FL,PSIP,PSIM,A,FAC,DFAC,DEFAC,RFAC,BI,CO,SI,T,CO2,SI2,SC,FA,&
|
||||
PLFACT,D2L12P,TN,PA,S,FA1,FA2
|
||||
LOGICAL FIRSTCALL
|
||||
DIMENSION BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2),&
|
||||
! THIS ARRAY IS OF THE FORM: BRAC(NP,N1P,MP,N1,N2)
|
||||
A(0:42,0:84),PSIP(85,85),PSIM(85,85),&
|
||||
FL(0:84,0:42,0:84),&
|
||||
! THIS ARRAY IS OF THE FORM: FL(NP,MP,N)
|
||||
FAC(0:169),DFAC(0:84),DEFAC(0:42),RFAC(0:84),BI(0:169,0:169)
|
||||
SAVE PSIP,PSIM,FL,FAC,DFAC,RFAC,A,BI
|
||||
IF (NQMAX.GT.84) THEN
|
||||
WRITE(6,*)'IN THE R(*8) COMPUTATION NQMAX SHOULD NOT EXCEED 84'
|
||||
STOP
|
||||
ENDIF
|
||||
IF (CO.EQ.0.D0.OR.SI.EQ.0.D0) THEN
|
||||
WRITE(6,*)'THE PROGRAM IS OF NO USE AT ZERO COS(PHI) OR SIN(PHI) VALUES'
|
||||
WRITE(6,*)'COS(PHI)=',CO,' SIN(PHI)=',SI
|
||||
STOP
|
||||
ENDIF
|
||||
IF (NQMAX.LT.L) THEN
|
||||
WRITE(6,*)'NQMAX=',NQMAX,' L=',L
|
||||
WRITE(6,*)'L SHOULD NOT EXCEED NQMAX'
|
||||
STOP
|
||||
ENDIF
|
||||
! BRACKETS [N_1'L_1'N_2'L_2'|0L_10L_2]_L^\PHI TO START THE RECURSION, EQ. (22) IN
|
||||
! THE ACCOMPANYING CPC PAPER.
|
||||
! SUBSIDIARY QUANTITIES:
|
||||
NQML=NQMAX-L
|
||||
NQMLD2=NQML/2
|
||||
NN=NQML-2*(NQMLD2)
|
||||
NPMAX=L-NN
|
||||
IF (M.GT.NQMLD2.OR.N.GT.NPMAX) THEN
|
||||
WRITE(6,*) 'M SHOULD NOT EXCEED (NQ-L)2 AND N SHOULD NOT EXCEED L-NN. BUT'
|
||||
WRITE(6,*) 'M=',M,' MMAX=',MPMAX,' N=',N,' NMAX=',NPMAX
|
||||
STOP
|
||||
ENDIF
|
||||
L1=M+N+NN
|
||||
L2=M-N+L
|
||||
L1L2=L1+L2
|
||||
L1L2ML=L1L2-L
|
||||
MPMAX=(L1L2ML-NN)/2
|
||||
FA=SQRT((2.D0*L1+1)*(2*L2+1))*CO**L1L2
|
||||
T=SI/CO
|
||||
IF (FIRSTCALL) THEN
|
||||
CALL ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
|
||||
CALL COE(2*NQMAX+1,FAC,BI)
|
||||
! THE A_{NL} ARRAY:
|
||||
DO NI=0,NQMAX/2
|
||||
NA=NI-2*(NI/2)
|
||||
K=1
|
||||
IF(NA.EQ.1)K=-1
|
||||
DO LI=0,NQMAX-2*NI
|
||||
A(NI,LI)=K/SQRT(DEFAC(NI)*DFAC(NI+LI))
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDIF
|
||||
PLFACT=SQRT(FAC(L1L2+L+1)*FAC(L1L2ML))
|
||||
IF (FIRSTCALL) CALL FLPHI(NQMAX,L,N,T,BI,RFAC,FL)
|
||||
! THIS SUBROUTINE PRODUCES THE ARRAY FL USED BELOW.
|
||||
! COMPUTATION OF THE N1=N2=0 BRACKETS. THEY ARE THE QUANTITIES BRAC(NP,N1P,MP,0,0):
|
||||
DO MP=0,MPMAX
|
||||
L12P=2*MP+L+NN
|
||||
! L12P=L1P+L2P
|
||||
D2L12P=1.D0/2.D0**L12P
|
||||
N12P=MPMAX-MP
|
||||
! N12P=N1P+N2P=(L1L2-L12P)/2
|
||||
DO N1P=0,N12P
|
||||
N1PA=N1P-2*(N1P/2)
|
||||
IF (N1PA.EQ.0) THEN
|
||||
K=1
|
||||
ELSE
|
||||
K=-1
|
||||
ENDIF
|
||||
! K=(-1)**N1P
|
||||
TN=T**N12P*K
|
||||
N2P=N12P-N1P
|
||||
DO NP=0,NPMAX
|
||||
L1P=MP+NP+NN
|
||||
L2P=L12P-L1P
|
||||
! BY DEFINITION ONE ALSO HAS: L2P=L2P(MP,NP)=MP-NP+L
|
||||
PA=A(N1P,L1P)*A(N2P,L2P)
|
||||
BRAC(NP,N1P,MP,0,0)=FL(NP,MP,N)*FA*TN*D2L12P*PLFACT*PA*PA
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
! RECURSION TO OBTAIN THE GENERAL FORM BRACKETS
|
||||
! [N_1'L_1'N_2'L_2'|N_1L_1N_2L_2]_L^\PHI.
|
||||
IF (FIRSTCALL) CALL COEFREC((NQMAX+L)/2+1,L,NN,NQMAX,PSIP,PSIM)
|
||||
N12MAX=(NQMAX-L1L2)/2
|
||||
MPMAX0=MPMAX
|
||||
CO2=CO**2
|
||||
SI2=SI**2
|
||||
SC=SI*CO
|
||||
! THE N1=0, N2-1-->N2 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER WITH THE
|
||||
! MODIFICATION POINTED OUT THERE:
|
||||
DO N2=1,N12MAX
|
||||
MPMAX=MPMAX+1
|
||||
! THIS MPMAX VALUE IS (NQ-L-NN)/2 WHERE NQ=L1+L2+2*N2.
|
||||
DO MP=0,MPMAX
|
||||
N12P=MPMAX-MP
|
||||
DO N1P=0,N12P
|
||||
! MP=(L1P+L2P-L-NN)/2. THEREFORE, N12P=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
|
||||
DO NP=0,NPMAX
|
||||
L1P=MP+NP+NN
|
||||
L2P=MP-NP+L
|
||||
S=0.D0
|
||||
! BELOW THE RESTRICTIONS ARE IMPOSED ON THE BRAC ARRAY ENTERING THE RIGHT-HAND
|
||||
! SIDE OF THE RECURRENCE FORMULAE. THIS ARRAY IS OF THE FORM BRAC(K1,K2,K3,...) WHERE
|
||||
! K1=K1(NP), K2=N2(N1P), AND K3=K3(MP). NAMELY, K1=NP, OR NP+1, OR NP-1;
|
||||
! K2=N1P OR N1P-1; K3=MP, OR MP-1, OR MP+1. (THE COMBINATION K2=N1P AND K3=MP+1
|
||||
! DOES NOT ARISE IN THE RECURRENCE FORMULAE.) THE QUANTITIES K1, K2, AND K3
|
||||
! REPRESENT, RESPECTIELY, NP, N1P, AND MP VALUES AT THE PRECEDING STAGE OF THE
|
||||
! RECURSION. THE RESTRICTIONS ENSURE THAT K1, K2, AND K3 RANGE WITHIN THE LIMITS
|
||||
! PERTAINING TO THE BRAC ARRAY OBTAINED AT THAT PRECEDING STAGE OF THE RECURSION.
|
||||
|
||||
! THUS K1, K2, AND K3 SHOULD BE NON-NEGATIVE. THEREFORE, WHEN K1=NP-1 THE VALUE
|
||||
! NP=0 IS TO BE EXCLUDED, WHEN K2=N1P-1 THE VALUE N1P=0 IS TO BE EXCLUDED, AND
|
||||
! WHEN K3=MP-1 THE VALUE MP=0 IS TO BE EXCLUDED.
|
||||
|
||||
! FURTHERMORE, IT SHOULD BE K1 \LE NMAX AND AT THE SAME TIME ONE HAS NP \LE NMAX.
|
||||
! THEREFORE, WHEN K1=NP+1 THE VALUE NP=NPMAX IS TO BE EXCLUDED.
|
||||
|
||||
! IN ADDITION, IT SHOULD BE K2+K3 \LE MPMAX-1 WHILE ONE HAS N1P+MP \LE MPMAX. WHEN
|
||||
! K2=N1P AND K3=MP-1, OR K2=N1P-1 AND K3=MP-1, OR K2=N1P-1 AND K3=MP ONE
|
||||
! AUTOMATICALLY HAS N1P+MP \LE MPMAX. THEREFORE, IN THESE CASES THE CONDITION
|
||||
! K2+K3 \LE MPMAX-1 DOES NOT CREATE ANY RESTRICTIONS.
|
||||
! BUT WHEN K2=N1P AND K3=MP, OR K2=N1P-1 AND K3=MP+1 THE CASE N1P+MP=MPMAX
|
||||
! IS TO BE FORBIDDEN.
|
||||
|
||||
! IN THE RECURRENCE FORMULAE BELOW THE DESCRIBED RESTRICTIONS ARE
|
||||
! IMPOSED.
|
||||
IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,0,N2-1)*PSIP(L1P,L2P)
|
||||
IF (N1P.NE.0.AND.N1P.NE.N12P) &
|
||||
S=S+BRAC(NP,N1P-1,MP+1,0,N2-1)*PSIP(L1P+1,L2P+1)
|
||||
! RECALL THAT N12P=MPMAX-MP.
|
||||
IF (NP.NE.NPMAX.AND.N1P.NE.0) &
|
||||
S=S-BRAC(NP+1,N1P-1,MP,0,N2-1)*PSIM(L1P+1,L2P)
|
||||
IF (NP.NE.0.AND.N1P.NE.N12P) &
|
||||
S=S-BRAC(NP-1,N1P,MP,0,N2-1)*PSIM(L1P,L2P+1)
|
||||
S=S*SC
|
||||
IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,0,N2-1)*SI2
|
||||
IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,0,N2-1)*CO2
|
||||
BRAC(NP,N1P,MP,0,N2)=S
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO ! N2
|
||||
! N1-1-->N1 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER:
|
||||
DO N2=0,N12MAX
|
||||
MPMAX=MPMAX0+N2
|
||||
! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2=MPMAX0+N2 SINCE NQ=L1+L2+2*N2.
|
||||
! THE RECURSION:
|
||||
DO N1=1,N12MAX-N2
|
||||
MPMAX=MPMAX+1
|
||||
! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2 AND NQ=L1+L2+2*N2+2*N1.
|
||||
DO MP=0,MPMAX
|
||||
N12P=MPMAX-MP
|
||||
DO N1P=0,N12P
|
||||
! MP=(L1P+L2P-L-NN)/2. THEREFORE, N1PMAX=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
|
||||
DO NP=0,NPMAX
|
||||
L1P=MP+NP+NN
|
||||
L2P=MP-NP+L
|
||||
S=0.D0
|
||||
IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,N1-1,N2)*PSIP(L1P,L2P)
|
||||
IF (N1P.NE.0.AND.N1P.NE.N12P) &
|
||||
S=S+BRAC(NP,N1P-1,MP+1,N1-1,N2)*PSIP(L1P+1,L2P+1)
|
||||
IF (NP.NE.NPMAX.AND.N1P.NE.0) &
|
||||
S=S-BRAC(NP+1,N1P-1,MP,N1-1,N2)*PSIM(L1P+1,L2P)
|
||||
IF (NP.NE.0.AND.N1P.NE.N12P) &
|
||||
S=S-BRAC(NP-1,N1P,MP,N1-1,N2)*PSIM(L1P,L2P+1)
|
||||
! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
|
||||
S=-S*SC
|
||||
IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,N1-1,N2)*CO2
|
||||
IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,N1-1,N2)*SI2
|
||||
! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
|
||||
BRAC(NP,N1P,MP,N1,N2)=S
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO ! N1
|
||||
ENDDO ! N2
|
||||
! RENORMALIZATION OF THE BRACKETS, EQ. (18) IN THE ACCOMPANYING CPC PAPER:
|
||||
DO N2=0,N12MAX
|
||||
MPMAX1=MPMAX0+N2
|
||||
DO N1=0,N12MAX-N2
|
||||
MPMAX=MPMAX1+N1
|
||||
FA1=A(N1,L1)*A(N2,L2)
|
||||
DO MP=0,MPMAX
|
||||
N12P=MPMAX-MP
|
||||
DO N1P=0,N12P
|
||||
N2P=N12P-N1P
|
||||
DO NP=0,NPMAX
|
||||
L1P=MP+NP+NN
|
||||
L2P=MP-NP+L
|
||||
FA2=A(N1P,L1P)*A(N2P,L2P)
|
||||
BRAC(NP,N1P,MP,N1,N2)=BRAC(NP,N1P,MP,N1,N2) &
|
||||
*FA1/FA2
|
||||
ENDDO ! NP
|
||||
ENDDO ! N1P
|
||||
ENDDO ! MP
|
||||
ENDDO ! N1
|
||||
ENDDO ! N2
|
||||
RETURN
|
||||
END
|
||||
|
||||
SUBROUTINE ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
|
||||
! FAC(I), DFAC(I),DEFAC(I), AND RFAC(I) ARE, RESPECTIVELY, THE QUANTITIES
|
||||
! I!, (2I+1)!!, (2I)!!, AND SQRT((2*I)!)/I!
|
||||
DOUBLE PRECISION FAC,DFAC,DEFAC,RFAC
|
||||
DIMENSION FAC(0:169),DFAC(0:84),DEFAC(0:42),RFAC(0:84)
|
||||
FAC(0)=1.D0
|
||||
DFAC(0)=1.D0
|
||||
DEFAC(0)=1.D0
|
||||
RFAC(0)=1.D0
|
||||
DO I=1,2*NQMAX+1
|
||||
FAC(I)=FAC(I-1)*I
|
||||
ENDDO
|
||||
DO I=1,NQMAX+1
|
||||
DFAC(I)=DFAC(I-1)*(2*I+1)
|
||||
ENDDO
|
||||
DO I=1,NQMAX/2+1
|
||||
DEFAC(I)=DEFAC(I-1)*2*I
|
||||
ENDDO
|
||||
DO I=1,NQMAX
|
||||
RFAC(I)=RFAC(I-1)*2*SQRT(1-0.5D0/I)
|
||||
ENDDO
|
||||
RETURN
|
||||
END
|
||||
|
||||
SUBROUTINE FLPHI(NQMAX,L,N,T,BI,RFAC,FL)
|
||||
! PROVIDES THE QUANTITY SQRT((2L1P+1)*(2L2P+1))*F_L^\VARPHI, SEE EQ. (23) IN THE CPC
|
||||
! PAPER, IN THE FORM OF THE FL ARRAY.
|
||||
! USES THE FUNCTION WIGMOD.
|
||||
! THIS FLPHI ROUTINE IS NOT IDENTICAL TO THAT IN THE OTHER FILE allosbrac.f90.
|
||||
! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
|
||||
DOUBLE PRECISION BI,RFAC,FL,T,T2,SQP,F,WIGMOD,z
|
||||
DIMENSION BI(0:169,0:169),RFAC(0:84),FL(0:84,0:42,0:84)
|
||||
! THE OUTPUT IS FL(NP,MP,N) WHERE NP=(L1P-L2P+L-NN)/2,
|
||||
! MP=(L1P+L2P-(L+NN))/2, L1P+L2P=L1T+L2T, AND N=(L1-L2+L-NN)/2, L1-L2=L1T-L2T.
|
||||
T2=T*T
|
||||
NQMAL=NQMAX-L
|
||||
MPMAX=NQMAL/2
|
||||
NN=NQMAL-2*MPMAX
|
||||
LMNN=L-NN
|
||||
L1ML2=2*N-LMNN
|
||||
LL3=2*L
|
||||
DO MP=0,MPMAX
|
||||
L1PL2P=2*MP+L+NN
|
||||
DO NP=0,LMNN
|
||||
L1PML2P=2*NP-LMNN
|
||||
L1P=(L1PL2P+L1PML2P)/2
|
||||
L2P=(L1PL2P-L1PML2P)/2
|
||||
M3=L2P-L1P
|
||||
SQP=SQRT((2*L1P+1.D0)*(2*L2P+1))
|
||||
L1T=(L1PL2P+L1ML2)/2
|
||||
L2T=(L1PL2P-L1ML2)/2
|
||||
N3=L1T+L2T-L
|
||||
IMIN=ABS(L1T-L1P)
|
||||
IMAX=MIN(L1P+L1T,L2P+L2T)
|
||||
NAL1=(L1T+L1P-IMIN)/2
|
||||
NAL2=(L1T-L1P+IMIN)/2
|
||||
NAL3=(L2T-L2P+IMIN)/2
|
||||
NAL4=(L2T+L2P-IMIN)/2
|
||||
NAL4P=NAL4-2*(NAL4/2)
|
||||
F=0.D0
|
||||
z=t**imin
|
||||
if(nal4p.ne.0)z=-z
|
||||
DO I=IMIN,IMAX,2
|
||||
F=F+RFAC(NAL1)*RFAC(NAL2)*RFAC(NAL3)*RFAC(NAL4)*&
|
||||
WIGMOD(L1T,L2T,L,L1P-I,I-L2P,BI)*z
|
||||
NAL1=NAL1-1
|
||||
NAL2=NAL2+1
|
||||
NAL3=NAL3+1
|
||||
NAL4=NAL4-1
|
||||
z=-z*t2
|
||||
ENDDO
|
||||
FL(NP,MP,N)=F*SQRT((2*L1P+1)*(2*L2P+1)&
|
||||
*BI(L1T+L-L2T,2*L1T)*BI(N3,2*L2T)/((LL3+1)*BI(N3,L1T+L2T+L+1)*BI(L+M3,LL3)))
|
||||
ENDDO
|
||||
ENDDO
|
||||
RETURN
|
||||
END
|
||||
|
||||
FUNCTION WIGMOD(L1,L2,L3,M1,M2,BI)
|
||||
! THE 3J SYMBOL IN TERMS OF BINOMIAL COEFFICIENTS WITHOUT THE FACTOR
|
||||
! SQRT(F) WHERE F IS AS FOLLOWS,
|
||||
! F=BI(L1+L3-L2,2*L1)*BI(N3,2*L2)/((LL3+1)*BI(N3,L1+L2+L3+1)*BI(L3+M3,LL3))
|
||||
! WITH N3=L1+L2-L3, LL3=2*L3, AND M3=-M1-M2.
|
||||
! THIS WIGMOD ROUTINE IS NOT IDENTICAL TO THAT IN THE OTHER FILE allosbrac.f90.
|
||||
DOUBLE PRECISION WIGMOD,BI,S
|
||||
DIMENSION BI(0:169,0:169)
|
||||
M3=-M1-M2
|
||||
N3=L1+L2-L3
|
||||
LM1=L1-M1
|
||||
LP2=L2+M2
|
||||
KMIN=MAX(0,L1-L3+M2,L2-L3-M1)
|
||||
KMAX=MIN(N3,LM1,LP2)
|
||||
S=0.d0
|
||||
NPH=1
|
||||
DO K=KMIN,KMAX
|
||||
S=S+NPH*BI(K,N3)*BI(LM1-K,L1+L3-L2)*BI(LP2-K,L2+L3-L1)
|
||||
NPH=-NPH
|
||||
ENDDO
|
||||
WIGMOD=S/SQRT(BI(LM1,2*L1)*BI(LP2,2*L2))
|
||||
NY=KMIN+L1-L2-M3
|
||||
NYP=NY-2*(NY/2)
|
||||
IF(NYP.NE.0)WIGMOD=-WIGMOD
|
||||
RETURN
|
||||
END
|
||||
|
||||
SUBROUTINE COE(NMAX,FAC,BI)
|
||||
! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
|
||||
DOUBLE PRECISION FAC,BI
|
||||
DIMENSION FAC(0:169),BI(0:169,0:169)
|
||||
DO N=0,NMAX
|
||||
DO M=0,N/2
|
||||
BI(M,N)=FAC(N)/(FAC(M)*FAC(N-M))
|
||||
BI(N-M,N)=BI(M,N)
|
||||
ENDDO
|
||||
ENDDO
|
||||
RETURN
|
||||
END
|
||||
|
||||
SUBROUTINE COEFREC(MAXCOE,L,NN,NQMAX,PSIP,PSIM)
|
||||
! CALCULATES THE COEFFICIENTS OF THE RECURSION FORMULA, EQ. (21) IN THE ACCOMPANYING
|
||||
! CPC PAPER
|
||||
DOUBLE PRECISION PSIP,PSIM
|
||||
DIMENSION PSIP(85,85),PSIM(85,85)
|
||||
! MAXCOE EQUALS INT((NQMAX+LMAX)/2)+1 AT CALLS OF THIS ROUTINE.
|
||||
LP=L+1
|
||||
LM=L-1
|
||||
DO M1=1,MAXCOE
|
||||
DO M2=1,M1
|
||||
M1P2=M1+M2
|
||||
NA=M1P2+L
|
||||
NNA=NA-2*(NA/2)
|
||||
M1M2=M1-M2
|
||||
! CALCULATION OF PSIP. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P AND L2P OR
|
||||
! L1P+1 AND L2P+1.)
|
||||
IF (NNA.EQ.NN.AND.M1P2.GE.L.AND.M1P2.LE.NQMAX.AND.&
|
||||
ABS(M1M2).LE.L) THEN
|
||||
IF (M1P2.GT.LP) THEN
|
||||
M1P2L=M1P2+L
|
||||
M1P2ML=M1P2-L
|
||||
PSIP(M1,M2)=SQRT(M1P2L*(M1P2L+1.D0)*M1P2ML*(M1P2ML-1)/((4*M1*M1-1)&
|
||||
*(4*M2*M2-1)))
|
||||
ELSE
|
||||
PSIP(M1,M2)=0.D0
|
||||
ENDIF
|
||||
PSIP(M2,M1)=PSIP(M1,M2)
|
||||
ENDIF
|
||||
! CALCULATION OF PSIM. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P+1 AND L2P
|
||||
! OR L2P+1 AND L1P.)
|
||||
IF (NNA.NE.NN.AND.M1P2.GE.LP.AND.M1P2.LE.NQMAX&
|
||||
.AND.M1M2.LE.LP.AND.M1M2.GE.-LM) THEN
|
||||
IF (M1M2.LT.L) THEN
|
||||
M1M2L=M1M2+L
|
||||
M1M2ML=M1M2-L
|
||||
PSIM(M1,M2)=SQRT(M1M2L*(M1M2L+1.D0)*M1M2ML*(M1M2ML-1)/&
|
||||
((4*M1*M1-1)*(4*M2*M2-1)))
|
||||
ELSE
|
||||
PSIM(M1,M2)=0.D0
|
||||
ENDIF
|
||||
PSIM(M2,M1)=PSIM(M1,M2)
|
||||
ENDIF
|
||||
ENDDO
|
||||
ENDDO
|
||||
RETURN
|
||||
END
|
||||
|
||||
|
|
@ -0,0 +1,260 @@
|
|||
program testallosbrac
|
||||
! Tests of the program that calculates the <N_1'L_1'N_2'L_2'|N_1L_1N_2L_2>_L^\varphi brackets.
|
||||
implicit double precision(a-h,o-z)
|
||||
dimension fac(0:201),dfac(0:201),defac(0:201),rfac(0:100)
|
||||
! These arrays enter the routine ARR called below. The first three arrays are used merely to
|
||||
! compute the exact expression for the bracket below. The array rfac is of no use.
|
||||
allocatable :: brac(:,:,:,:,:,:,:,:),bi(:,:)
|
||||
! Below the allocation/deallocation actions are done independently for each of the tests
|
||||
! performed as it would be required in the case when only a single test is being performed.
|
||||
open(13,file='allosoutput')
|
||||
! The input angular parameters are CO=COS(VARPHI) and SI=SIN(VARPHI). Set, for example,
|
||||
si=1/sqrt(2.d0)
|
||||
co=-si
|
||||
! This corresponds to Eq. (13) in the accompanying CPC paper in the equal mass case. The
|
||||
! corresponding transformation (7) in that paper arises at calculating standard two-body
|
||||
! shell-model matrix elements.
|
||||
! One might also test the program for COS(VARPHI) and SIN(VARPHI) values pertaining to
|
||||
! nonequal mass cases. For example, one might set
|
||||
! si=1/sqrt(3.d0)
|
||||
! co=-sqrt(2/3.d0)
|
||||
! 1. Comparison of the value of the <n1p,L00|n1,l1,n2,l2>_L^\varphi bracket with the result
|
||||
! of its analytic calculation. Here l1+l2 should be of the same parity as L.
|
||||
l=4
|
||||
l1=2
|
||||
l2=4
|
||||
n1=0
|
||||
n2=0
|
||||
nq=2*n1+2*n2+l1+l2
|
||||
! also nq=2*n1p+L
|
||||
NN=0
|
||||
! this is because nq-l=2*n1p is even
|
||||
mp=0
|
||||
! here l1p=l, mp=(l1p+l2p-l-nn)/2=0
|
||||
np=l
|
||||
! np=(l1p-l2p+l-nn)/2=l
|
||||
n1p=(nq-l)/2
|
||||
! Computation of the exact bracket.
|
||||
nfac=2*nq+1
|
||||
CALL ARR(FAC,DFAC,DEFAC,RFAC,NQ)
|
||||
! The array of binomial coefficients to calculate the Wigner coefficients:
|
||||
allocate (BI(0:NFAC,0:NFAC),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate BI failed"
|
||||
call coe(NFAC,FAC,BI)
|
||||
nm=n1p-n1-n2
|
||||
nfa=nm-2*(nm/2)
|
||||
LL3=2*L
|
||||
N3=L1+L2-L
|
||||
! exact bracket:
|
||||
b=co**(2*n1+l1)*si**(2*n2+l2)*(-1)**l*sqrt((2.d0*l1+1)*(2*l2+1))*(-1)**nfa*&
|
||||
WIGMOD(l1,l2,l,0,0,bi,nfac)*SQRT(BI(L1+L-L2,2*L1)*BI(N3,2*L2)/&
|
||||
((LL3+1)*BI(N3,L1+L2+L+1)*BI(L,LL3)))&
|
||||
*sqrt(defac(n1p)*dfac(n1p+l)/(defac(n1)*defac(n2)*dfac(n1+l1)*dfac(n2+l2)))
|
||||
! bracket from the code:
|
||||
ndel=4
|
||||
nqmax=nq+ndel
|
||||
! The result should not depend on ndel, provided that ndel is even, i.e., nq and nqmax are of the
|
||||
! same parity.
|
||||
NQML=NQMAX-L
|
||||
m=(l1+l2-l-nn)/2
|
||||
n=(l1-l2+l-nn)/2
|
||||
allocate (BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,&
|
||||
0:(NQMAX-L)/2,0:L,0:(NQMAX-L)/2,L:L),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
|
||||
CALL ALLOSBRAC(NQMAX,L,L,CO,SI,BRAC)
|
||||
a=brac(np,n1p,mp,n1,n2,n,m,l)
|
||||
write(6,*)'1. bracket <n1p L 0 0|n1 l1 n2 l2>_L^varphi'
|
||||
write(6,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l
|
||||
write(6,*)'exact bracket',b,' bracket from the code',a
|
||||
write(13,*)'1. bracket <n1p L 0 0|n1 l1 n2 l2>_L^varphi'
|
||||
write(13,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l
|
||||
write(13,*)'exact bracket',b,' bracket from the code',a
|
||||
! write(6,*)co,si
|
||||
deallocate (BRAC,STAT=istatus)
|
||||
! 2. Test of the symmetry of the brackets. The relation
|
||||
! <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\varphi=<N1,L1,N2,L2|N1P,L1P,N2P,L2P>_L^\varphi
|
||||
! should hold true. This symmetry relation may be rewritten as
|
||||
! BRAC(NP,N1P,MP,N1,N2,N,M,L)=BRAC1(N,N1,M,N1P,N2P,NP,MP,L) where N and
|
||||
! M represent the angular momenta L1 and L2 and NP and MP represent L1P and L2P.
|
||||
! An example:
|
||||
nq=30
|
||||
l=9
|
||||
! nn=1 in the present example.
|
||||
! calculate: nmax=l-nn.
|
||||
! nmax=8 in the present example. Any n and np that do not exceed nmax may be employed.
|
||||
n=4
|
||||
np=3
|
||||
! calculate: nqmld2=(nq-l-nn)/2.
|
||||
! nqmld2=10 in the present example. Any m, n1, n2, mp, n1p, and n2p may be employed
|
||||
! provided that m+n1+n2=mp+n1p+n2p=nqmld2.
|
||||
m=5
|
||||
n1=3
|
||||
n2=2
|
||||
mp=4
|
||||
n1p=3
|
||||
n2p=3
|
||||
nqmax=32
|
||||
! Recall that nq=30. Any nqmax of the same parity as nq and such that nqmax.ge.nq
|
||||
! may be employed. Then its choice does not influence the results.
|
||||
allocate (BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,&
|
||||
0:(NQMAX-L)/2,0:L,0:(NQMAX-L)/2,L:L),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
|
||||
CALL ALLOSBRAC(NQMAX,L,L,CO,SI,BRAC)
|
||||
write(6,*)'2. Symmetry test, L=',l
|
||||
write(6,*)'m =',m,' n =',n,' n1 =',n1,' n2 =',n2
|
||||
write(6,*)'mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
|
||||
write(6,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =',BRAC(np,n1p,mp,n1,n2,n,m,l)
|
||||
write(6,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =',BRAC(n,n1,m,n1p,n2p,np,mp,l)
|
||||
write(13,*)'2. Symmetry test, L=',l
|
||||
write(13,*)'m =',m,' n =',n,' n1 =',n1,' n2 =',n2
|
||||
write(13,*)'mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
|
||||
write(13,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =',BRAC(np,n1p,mp,n1,n2,n,m,l)
|
||||
write(13,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =',BRAC(n,n1,m,n1p,n2p,np,mp,l)
|
||||
deallocate (BRAC,STAT=istatus)
|
||||
! 3. Test of the relation
|
||||
! \sum_{i,i',L}\sum_j<j|i><j|i'>=\sum_L N_0(L).
|
||||
! Here i=(n_1,l_1,n_2,l_2), i'=(n_3,l_3,n_4,l_4), and j=(n_1p,l_1p,n_2p,l_2p). The l.h.s. sum
|
||||
! runs over all i and i' values such that nq(i)=nq(i')=nq.
|
||||
! In the r.h.s. sum, N_0(L) is the number of the (n_1,l_1,n_2,l_2) states pertaining to
|
||||
! given nq and L values. The sum runs over all the L values compatible with a given nq value
|
||||
! and thus represents the total number of states with this nq.
|
||||
nq=14
|
||||
ndel=4
|
||||
nqmax=nq+ndel
|
||||
! The result should not depend on ndel, i.e., on nqmax, provided that ndel is even.
|
||||
numb=0
|
||||
s=0.d0
|
||||
LMIN=0
|
||||
LMAX=NQ
|
||||
allocate (BRAC(0:LMAX,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,&
|
||||
0:(NQMAX-LMIN)/2,0:LMAX,0:(NQMAX-LMIN)/2,LMIN:LMAX),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
|
||||
CALL ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC)
|
||||
do L=0,nq
|
||||
nqml=nq-l
|
||||
nqmld2=nqml/2
|
||||
nn=nqml-2*nqmld2
|
||||
mmax=(nq-(l+nn))/2
|
||||
nmax=l-nn
|
||||
nst=(nmax+1)*(mmax+1)*(mmax+2)/2
|
||||
do na=0,nmax
|
||||
do ma=0,mmax
|
||||
! IN THE ABOVE NOTATION NA=NA(L1,L2), MA=MA(L1,L2).
|
||||
n12=mmax-ma
|
||||
do nb=0,nmax
|
||||
do mb=0,mmax
|
||||
! IN THE ABOVE NOTATION NB=NB(L3,L4), MB=MB(L4,L4).
|
||||
n34=mmax-mb
|
||||
ds=0.d0
|
||||
do n2=0,n12
|
||||
n1=n12-n2
|
||||
do n4=0,n34
|
||||
n3=n34-n4
|
||||
DO MP=0,mmax
|
||||
do n1p=0,mmax-mp
|
||||
DO NP=0,nmax
|
||||
DS=DS+BRAC(np,n1p,mp,n1,n2,na,ma,l)&
|
||||
*BRAC(np,n1p,mp,n3,n4,nb,mb,l)
|
||||
enddo ! NP
|
||||
ENDDO ! n1p
|
||||
ENDDO ! MP
|
||||
enddo ! n4
|
||||
enddo ! n2
|
||||
s=s+ds
|
||||
enddo ! MB
|
||||
enddo ! NB
|
||||
enddo ! MA
|
||||
enddo ! NA
|
||||
numb=numb+nst
|
||||
ENDDO !l
|
||||
deallocate (BRAC,STAT=istatus)
|
||||
write(6,*)'3. Test of the relation sum_{i,ip,L}sum_j<j|i><j|ip> = the number of states.'
|
||||
write(6,*)'See the comments in the text of the program and in more detail in the'
|
||||
write(6,*)'accompanying CPC paper for further explanations.'
|
||||
write(6,*)'nqmax =',nqmax,' nq =',nq
|
||||
write(6,*)'sum_{i1,i2,L}sum_j<j|i1><j|i2> =',s,' exact value (number of states) =', numb
|
||||
write(13,*)'3. Test of the relation sum_{i,ip,L}sum_j<j|i><j|ip> = the number of states.'
|
||||
write(13,*)'See the comments in the text of the program and in more detail in the'
|
||||
write(13,*)'accompanying CPC paper for further explanations.'
|
||||
write(13,*)'nqmax =',nqmax,' nq =',nq
|
||||
write(13,*)'sum_{i1,i2,L}sum_j<j|i1><j|i2> =',s,' exact value (number of states)=',numb
|
||||
! 4. Test of the relation \sum_{i,i',L}|\sum_j<j|i><j|i'>-\delta_{i,i'}|=0.
|
||||
! Here i=(n_1,l_1,n_2,l_2), i'=(n_1',l_1',n_2',l_2'), and j=(n_1'',l_1'',n_2'',l_2''). The outer
|
||||
! sum runs over all i and i' values such that nq(i)=nq(i') \le nqmax.
|
||||
! Also one has nq(i)=nq(i')=nq(j).
|
||||
nqmax=12
|
||||
lmin=0
|
||||
lmax=nqmax
|
||||
nqml=nqmax-lmin
|
||||
nqmld2=nqml/2
|
||||
allocate (BRAC(0:LMAX,0:nqmld2,0:nqmld2,0:nqmld2,0:nqmld2,0:LMAX,&
|
||||
0:nqmld2,LMIN:LMAX),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
|
||||
CALL ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC)
|
||||
! When nqmax is sufficiently high, this calculation is to be performed in somewhat another way
|
||||
! because of the memory restrictions. In such cases one may use ALLOSBRAC at lmin=lmax=l.
|
||||
! The BRAC array is to be allocated and deallocated inside the loop over l below, and
|
||||
! ALLOSBRAC is to be called inside this loop.
|
||||
s=0.d0
|
||||
! The quantity s represents the expression \sum_{i,i',L}|\sum_j<j|i><j|i'>-\delta_{i,i'}| we are
|
||||
! calculating.
|
||||
do l=0,nqmax
|
||||
nqml=nqmax-l
|
||||
nqmld2=nqml/2
|
||||
nn=nqml-2*nqmld2
|
||||
nmax=l-nn
|
||||
lpnn=l+nn
|
||||
do nq=lpnn,nqmax,2
|
||||
mmax=(nq-lpnn)/2
|
||||
do ma=0,mmax
|
||||
! The equality mmax=(nq-lpnn)/2 follows from the definition of m-type variables.
|
||||
n12=mmax-ma
|
||||
! The quantity n12 designates n1+n2 and here it is taken into account that
|
||||
! n1+n2+ma=mmax.
|
||||
do na=0,nmax
|
||||
! In the case of odd nqmax and l equal to zero, one has nmax=-1 and no contribution arises as it
|
||||
! should be.
|
||||
do n1=0,n12
|
||||
n2=n12-n1
|
||||
do mb=0,mmax
|
||||
n34=mmax-mb
|
||||
do nb=0,nmax
|
||||
do n3=0,n34
|
||||
n4=n34-n3
|
||||
ds=0.d0
|
||||
! The quantity ds represents the contribution of \sum_j|<j|i>_L<j|i'>_L-\delta_{i,i'}| to the net
|
||||
! result.
|
||||
DO MP=0,MMAX
|
||||
do n1p=0,mmax-mp
|
||||
DO NP=0,NMAX
|
||||
DS=DS+BRAC(np,n1p,mp,n1,n2,na,ma,l)&
|
||||
*BRAC(np,n1p,mp,n3,n4,nb,mb,l)
|
||||
enddo ! NP
|
||||
ENDDO ! n1p
|
||||
ENDDO ! MP
|
||||
if(ma.eq.mb.and.na.eq.nb.and.n1.eq.n3.and.n2.eq.n4)then
|
||||
ds=abs(ds-1.d0)
|
||||
else
|
||||
ds=abs(ds)
|
||||
endif
|
||||
s=s+ds
|
||||
enddo ! n3
|
||||
enddo ! nb
|
||||
enddo ! mb
|
||||
enddo ! n1
|
||||
enddo ! na
|
||||
enddo ! ma
|
||||
enddo ! nq
|
||||
ENDDO ! l
|
||||
deallocate (BRAC,STAT=istatus)
|
||||
write(6,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
|
||||
write(6,*)'For further explanations see the comments in the text of the program'
|
||||
write(6,*)'and, in more detail, in the accompanying CPC paper.'
|
||||
write(6,*)'nqmax =',nqmax
|
||||
write(6,*)'sum_{i1,i2,L}|sum_j<j|i1><j|i2>-delta_{i1,i2}| =',s
|
||||
write(13,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
|
||||
write(13,*)'For further explanations see the comments in the text of the program'
|
||||
write(13,*)'and, in more detail, in the accompanying CPC paper.'
|
||||
write(13,*)'nqmax =',nqmax
|
||||
write(13,*)'sum_{i1,i2,L}|sum_j<j|i1><j|i2>-delta_{i1,i2}| =',s
|
||||
end
|
||||
|
|
@ -0,0 +1,266 @@
|
|||
program testosbrac
|
||||
! Tests of the program that calculates the <N_1'L_1'N_2'L_2'|N_1L_1N_2L_2>_L^\varphi brackets.
|
||||
implicit double precision(a-h,o-z)
|
||||
logical firstcall
|
||||
dimension fac(0:169),dfac(0:84),defac(0:42),rfac(0:84),bi(0:169,0:169)
|
||||
! These arrays enter the routine ARR called below. The first three arrays are used merely to
|
||||
! compute the exact expression for the bracket below. The array rfac is of no use.
|
||||
allocatable :: brac(:,:,:,:,:),brac1(:,:,:,:,:)
|
||||
open(13,file='osoutput')
|
||||
! The input angular parameters are CO=COS(VARPHI) and SI=SIN(VARPHI). Set, for example,
|
||||
si=1/sqrt(2.d0)
|
||||
co=-si
|
||||
! This corresponds to Eq. (13) in the accompanying CPC paper in the equal mass case. The
|
||||
! corresponding transformation (7) in that paper arises at calculating standard two-body
|
||||
! shell-model matrix elements.
|
||||
! One might also test the program for COS(VARPHI) and SIN(VARPHI) values pertaining to
|
||||
! nonequal mass cases. For example, one might set
|
||||
! si=1/sqrt(3.d0)
|
||||
! co=-sqrt(2/3.d0)
|
||||
! 1. Comparison of the value of the <n1p,L00|n1,l1,n2,l2>_L^\varphi bracket with the result
|
||||
! of its analytic calculation. Here l1+l2 should be of the same parity as L.
|
||||
l=4
|
||||
l1=2
|
||||
l2=4
|
||||
n1=3
|
||||
n2=2
|
||||
mp=0
|
||||
np=l
|
||||
nq=2*n1+2*n2+l1+l2
|
||||
n1p=(nq-l)/2
|
||||
! Computation of the exact bracket.
|
||||
nfac=2*nq+1
|
||||
CALL ARR(FAC,DFAC,DEFAC,RFAC,NQ)
|
||||
! The array of binomial coefficients to calculate the Wigner coefficients:
|
||||
call coe(NFAC,FAC,BI)
|
||||
nm=n1p-n1-n2
|
||||
nfa=nm-2*(nm/2)
|
||||
LL3=2*L
|
||||
N3=L1+L2-L
|
||||
! exact bracket:
|
||||
b=co**(2*n1+l1)*si**(2*n2+l2)*(-1)**l*sqrt((2.d0*l1+1)*(2*l2+1))*(-1)**nfa*&
|
||||
WIGMOD(l1,l2,l,0,0,bi)*SQRT(BI(L1+L-L2,2*L1)*BI(N3,2*L2)/&
|
||||
((LL3+1)*BI(N3,L1+L2+L+1)*BI(L,LL3)))&
|
||||
*sqrt(defac(n1p)*dfac(n1p+l)/(defac(n1)*defac(n2)*dfac(n1+l1)*dfac(n2+l2)))
|
||||
! bracket from the code:
|
||||
ndel=4
|
||||
nqmax=nq+ndel
|
||||
! The result should not depend on ndel, provided that ndel is even, i.e., nq and nqmax are of the
|
||||
! same parity.
|
||||
NQML=NQMAX-L
|
||||
NN=0
|
||||
m=(l1+l2-l-nn)/2
|
||||
n=(l1-l2+l-nn)/2
|
||||
allocate (BRAC(0:L,0:NQML/2,0:NQML/2,0:NQML/2,0:NQML/2),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
|
||||
FIRSTCALL=.TRUE.
|
||||
CALL OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
|
||||
a=brac(np,n1p,mp,n1,n2)
|
||||
write(6,*)'1. bracket <n1p,L00|n1l1n2l2>_L^phi'
|
||||
write(6,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l
|
||||
write(13,*)'1. bracket <n1p,L00|n1l1n2l2>_L^phi'
|
||||
write(13,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l
|
||||
write(6,*)'exact bracket',b,' bracket from the code',a
|
||||
write(13,*)'exact bracket',b,' bracket from the code',a
|
||||
deallocate (BRAC,STAT=istatus)
|
||||
! 2. Test of the symmetry of the brackets. The relation
|
||||
! <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\varphi=<N1,L1,N2,L2|N1P,L1P,N2P,L2P>_L^\varphi
|
||||
! should hold true. Let the bracket in the left-hand side of this relation is contained in the BRAC
|
||||
! array and that in its right-hand side in the BRAC1 array. Then the above symmetry relation may
|
||||
! be rewritten as BRAC(NP,N1P,MP,N1,N2)=BRAC1(N,N1,M,N1P,N2P) where N and M represent
|
||||
! the angular momenta L1 and L2 and NP and MP represent L1P and L2P.
|
||||
! An example:
|
||||
nq=30
|
||||
l=9
|
||||
! nn=1 at this choice.
|
||||
! calculate: nmax=l-nn.
|
||||
! nmax=8 in the present example. Any n and np values
|
||||
! that do not exceed nmax may be employed.
|
||||
n=4
|
||||
np=3
|
||||
! calculate: nqmld2=(nq-l-nn)/2.
|
||||
! nqmld2=10 in the present example. Any m, n1, n2, mp, n1p, and n2p may be employed
|
||||
! provided that m+n1+n2=mp+n1p+n2p=nqmld2.
|
||||
m=5
|
||||
n1=3
|
||||
n2=2
|
||||
mp=4
|
||||
n1p=3
|
||||
n2p=3
|
||||
nqmax=36
|
||||
! Recall that nq=30. Any nqmax of the same parity as nq and such that nqmax.ge.nq
|
||||
! may be employed. Then its choice does not influence the results.
|
||||
nqml=nqmax-l
|
||||
allocate (BRAC(0:L,0:NQML/2,0:NQML/2,0:NQML/2,0:NQML/2),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
|
||||
CALL OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
|
||||
allocate (BRAC1(0:L,0:NQML/2,0:NQML/2,0:NQML/2,0:NQML/2),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC1) failed"
|
||||
CALL OSBRAC(NP,MP,L,NQMAX,CO,SI,FIRSTCALL,BRAC1)
|
||||
write(6,*)'2. Symmetry test, L =',l
|
||||
write(6,*)'m =',m,' n =',n,' n1=',n1,' n2=',n2,' mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
|
||||
write(6,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =', BRAC(np,n1p,mp,n1,n2)
|
||||
write(6,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =', BRAC1(n,n1,m,n1p,n2p)
|
||||
write(13,*)'2. Symmetry test, L=',l
|
||||
write(13,*)'m =',m,' n =',n,' n1=',n1,' n2=',n2,' mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
|
||||
write(13,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =', BRAC(np,n1p,mp,n1,n2)
|
||||
write(13,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =', BRAC1(n,n1,m,n1p,n2p)
|
||||
deallocate (BRAC1,STAT=istatus)
|
||||
deallocate (BRAC,STAT=istatus)
|
||||
! 3. Test of the relation
|
||||
! \sum_{i,i',L}\sum_j<j|i>_L<j|i'>_L=\sum_L N_0(L).
|
||||
! Here i=(n_1,l_1,n_2,l_2), i'=(n_3,l_3,n_4,l_4), and j=(n_1p,l_1p,n_2p,l_2p). The l.h.s. sum
|
||||
! runs over all i and i' values such that nq(i)=nq(i')=nq.
|
||||
! In the r.h.s. sum, N_0(L) is the number of the (n_1,l_1,n_2,l_2) states pertaining to
|
||||
! given nq and L values. The sum runs over all the L values compatible with a given nq value
|
||||
! and thus represents the total number of states with this nq.
|
||||
nq=14
|
||||
ndel=6
|
||||
nqmax=nq+ndel
|
||||
! The result should not depend on ndel, i.e., on nqmax, provided that ndel is even.
|
||||
numb=0
|
||||
s=0.d0
|
||||
do l=0,nq
|
||||
FIRSTCALL=.true.
|
||||
nqml=nqmax-l
|
||||
nqmld2=nqml/2
|
||||
nn=nqml-2*nqmld2
|
||||
mmax=(nq-(l+nn))/2
|
||||
nmax=l-nn
|
||||
nst=(nmax+1)*(mmax+1)*(mmax+2)/2
|
||||
allocate (BRAC(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
|
||||
allocate (BRAC1(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC1) failed"
|
||||
! The allocation/deallocation of BRAC and BRAC1, at nqmld2-->nqmax/2, might also be done
|
||||
! outside the loop.
|
||||
do ma=0,mmax
|
||||
n12=mmax-ma
|
||||
do na=0,nmax
|
||||
! IN THE ABOVE NOTATION NA=NA(L1,L2) AND MA=MA(L1,L2).
|
||||
CALL OSBRAC(NA,MA,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
|
||||
! The order of summations below aims to avoid unnecessary calls of osbrac.
|
||||
do mb=0,mmax
|
||||
n34=mmax-mb
|
||||
do nb=0,nmax
|
||||
! IN THE ABOVE NOTATION NB=NB(L3,L4) AND MB=MB(L3,L4).
|
||||
CALL OSBRAC(NB,MB,L,NQMAX,CO,SI,FIRSTCALL,BRAC1)
|
||||
ds=0.d0
|
||||
do n2=0,n12
|
||||
n1=n12-n2
|
||||
do n4=0,n34
|
||||
n3=n34-n4
|
||||
DO MP=0,mmax
|
||||
do n1p=0,mmax-mp
|
||||
DO NP=0,nmax
|
||||
DS=DS+BRAC(NP,n1p,MP,n1,n2)*BRAC1(NP,n1p,MP,n3,n4)
|
||||
enddo ! NP
|
||||
ENDDO ! n1p
|
||||
ENDDO ! MP
|
||||
enddo ! n4
|
||||
enddo ! n2
|
||||
s=s+ds
|
||||
enddo ! nb
|
||||
FIRSTCALL=.false.
|
||||
! If the nb loop is nested inside the mb loop the FIRSTCALL=.false. command is more efficient.
|
||||
! In practice, usually the gain in the running time is relatively not large.
|
||||
enddo ! mb
|
||||
enddo ! ma
|
||||
enddo ! na
|
||||
numb=numb+nst
|
||||
deallocate (BRAC,STAT=istatus)
|
||||
deallocate (BRAC1,STAT=istatus)
|
||||
ENDDO ! l
|
||||
write(6,*)'3. Test of the relation sum_{i,ip,L,j}<j|i>_L<j|ip>_L = the number of states.'
|
||||
write(6,*)'See the comments in the text of the program and in more detail in the'
|
||||
write(6,*)'accompanying CPC paper for further explanations.'
|
||||
write(6,*)'nqmax =',nqmax,' nq =',nq
|
||||
write(6,*)'sum_{i,ip,L,j}<j|i>_L<j|ip>_L =',s,' exact value (number of states) =',numb
|
||||
write(13,*)'3. Test of the relation sum_{i,ip,L}sum_j<j|i><j|ip> = the number of states.'
|
||||
write(13,*)'See the comments in the text of the program and in more detail in the'
|
||||
write(13,*)'accompanying CPC paper for further explanations.'
|
||||
write(13,*)'nqmax =',nqmax,' nq =',nq
|
||||
write(13,*)'sum_{i,ip,L,j}<j|i>_L<j|ip>_L =',s,' exact value (number of states) =',numb
|
||||
! 4. Test of the relation \sum_{i,i',L}\sum_j|<j|i>_L<j|i'>_L-\delta_{i,i'}|=0.
|
||||
! Here i=(n_1,l_1,n_2,l_2), i'=(n_1',l_1',n_2',l_2'), and j=(n_1'',l_1'',n_2'',l_2''). The outer
|
||||
! sum runs over all i and i' values such that nq(i)=nq(i') \le nqmax.
|
||||
! Also one has nq(i)=nq(i')=nq(j).
|
||||
nqmax=12
|
||||
s=0.d0
|
||||
! The quantity s represents the expression \sum_{i,i',L}|\sum_j<j|i><j|i'>-\delta_{i,i'}| we are
|
||||
! calculating.
|
||||
do l=0,nqmax
|
||||
FIRSTCALL=.true.
|
||||
nqml=nqmax-l
|
||||
nqmld2=nqml/2
|
||||
nn=nqml-2*nqmld2
|
||||
nmax=l-nn
|
||||
lpnn=l+nn
|
||||
allocate (BRAC(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
|
||||
allocate (BRAC1(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
|
||||
IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC1) failed"
|
||||
! The allocation/deallocation of BRAC and BRAC1, at nqmld2-->nqmax/2, might also be done
|
||||
! outside the loop.
|
||||
do ma=0,nqmld2
|
||||
! The above upper limit of summation corresponds to the definition of m-type
|
||||
! variables.
|
||||
do na=0,nmax
|
||||
! In the case of odd nqmax and l equal to zero, one has nmax=-1 and no contribution arises as it
|
||||
! should be.
|
||||
CALL OSBRAC(NA,MA,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
|
||||
! The order of summations below aims to avoid unnecessary calls of osbrac.
|
||||
do mb=0,nqmld2
|
||||
do nb=0,nmax
|
||||
CALL OSBRAC(NB,MB,L,NQMAX,CO,SI,FIRSTCALL,BRAC1)
|
||||
do nq=l+nn,nqmax,2
|
||||
nqcmld2=(nq-l-nn)/2
|
||||
n12=nqcmld2-ma
|
||||
n34=nqcmld2-mb
|
||||
! We consider nq to be the number of quanta pertaining to states entering the brackets in our
|
||||
! sum. Then one has: ma+n1+n2=mb+n3+n4=mp+n1p+n2p=(nq-l-nn)/2 that is
|
||||
! nqcmld2.
|
||||
do n1=0,n12
|
||||
n2=n12-n1
|
||||
do n3=0,n34
|
||||
n4=n34-n3
|
||||
ds=0.d0
|
||||
! The quantity ds represents the contribution of \sum_j|<j|i>_L<j|i'>_L-\delta_{i,i'}| to the net
|
||||
! result.
|
||||
DO MP=0,nqcmld2
|
||||
do n1p=0,nqcmld2-mp
|
||||
DO NP=0,NMAX
|
||||
DS=DS+BRAC(NP,n1p,MP,n1,n2)*BRAC1(NP,n1p,MP,n3,n4)
|
||||
enddo ! NP
|
||||
ENDDO ! n1p
|
||||
ENDDO ! MP
|
||||
if(ma.eq.mb.and.na.eq.nb.and.n1.eq.n3.and.n2.eq.n4)then
|
||||
ds=abs(ds-1.d0)
|
||||
else
|
||||
ds=abs(ds)
|
||||
endif
|
||||
s=s+ds
|
||||
enddo ! n3
|
||||
enddo ! n1
|
||||
enddo ! nq
|
||||
enddo ! nb
|
||||
FIRSTCALL=.false.
|
||||
! If the nb loop is nested inside the mb loop the FIRSTCALL=.false. command is more efficient.
|
||||
! In practice, usually the gain in the running time is relatively not large.
|
||||
enddo ! mb
|
||||
enddo ! ma
|
||||
enddo ! na
|
||||
deallocate (BRAC,STAT=istatus)
|
||||
deallocate (BRAC1,STAT=istatus)
|
||||
ENDDO ! l
|
||||
write(6,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
|
||||
write(6,*)'For further explanations see the comments in the text of the program'
|
||||
write(6,*)'and, in more detail, in the accompanying CPC paper.'
|
||||
write(6,*)'nqmax =',nqmax
|
||||
write(6,*)'sum_{i1,i2,L}|sum_j<j|i1>_L<j|i2>_L-delta_{i1,i2}| =',s
|
||||
write(13,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
|
||||
write(13,*)'For further explanations see the comments in the text of the program'
|
||||
write(13,*)'and, in more detail, in the accompanying CPC paper.'
|
||||
write(13,*)'nqmax =',nqmax
|
||||
write(13,*)'sum_{i1,i2,L}|sum_j<j|i1>_L<j|i2>_L-delta_{i1,i2}| =',s
|
||||
end
|
||||
|
|
@ -0,0 +1,43 @@
|
|||
# compile using
|
||||
# gfortran -shared -fPIC osbrac.f90 -o osbrac.so
|
||||
# gfortran -shared -fPIC allosbrac.f90 -o allosbrac.so
|
||||
|
||||
println("OSBRACKETS test against Brody et al.")
|
||||
|
||||
Emax = 3
|
||||
Λ = 1
|
||||
l1 = 1
|
||||
l2 = 2
|
||||
|
||||
ϵ = (Emax - Λ) % 2
|
||||
N = (l1 - l2 + Λ - ϵ) ÷ 2
|
||||
M = (l1 + l2 - Λ - ϵ) ÷ 2
|
||||
L = Λ
|
||||
NQMAX = Emax
|
||||
CO = 1/sqrt(2)
|
||||
SI = 1/sqrt(2)
|
||||
FIRSTCALL = true
|
||||
# from source: BRAC(NP,N1P,MP,N1,N2) with dimensions BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2)
|
||||
BRAC = zeros(Float64, L + 1, (NQMAX - L) ÷ 2 + 1,(NQMAX - L) ÷ 2 + 1,(NQMAX - L) ÷ 2 + 1,(NQMAX - L) ÷ 2 + 1)
|
||||
|
||||
@ccall "./OSBRACKETS/osbrac.so".osbrac_(N::Ref{Int32},M::Ref{Int32},L::Ref{Int32},NQMAX::Ref{Int32},CO::Ref{Float64},SI::Ref{Float64},FIRSTCALL::Ref{UInt8},BRAC::Ptr{Array{Float64}})::Cvoid
|
||||
|
||||
function get_bracket(n1′, l1′, n2′, l2′, n1, n2)
|
||||
Nq = l1 + l2 + 2 * (n1 + n2)
|
||||
Nq′ = l1′ + l2′ + 2 * (n1′ + n2′)
|
||||
if Nq ≠ Nq′; return 0; end
|
||||
N′ = (l1′ - l2′ + Λ - ϵ) ÷ 2
|
||||
N1′ = n1′
|
||||
M′ = (l1′ + l2′ - Λ - ϵ) ÷ 2
|
||||
N1 = n1
|
||||
N2 = n2
|
||||
return BRAC[N′ + 1, N1′ + 1, M′ + 1, N1 + 1, N2 + 1]
|
||||
end
|
||||
|
||||
test_n = [0 0 0 0 1 1]
|
||||
test_l = [0 1 1 2 0 1]
|
||||
test_N = [1 0 1 0 0 0]
|
||||
test_L = [1 2 0 1 1 0]
|
||||
|
||||
bracs = get_bracket.(test_n, test_l, test_N, test_L, 0, 0)
|
||||
display(bracs)
|
||||
Loading…
Reference in New Issue