diff --git a/OSBRACKETS/README b/OSBRACKETS/README deleted file mode 100755 index c5e55f0..0000000 --- a/OSBRACKETS/README +++ /dev/null @@ -1,71 +0,0 @@ -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 .... - - - - - - - - - - diff --git a/OSBRACKETS/allosbrac.f90 b/OSBRACKETS/allosbrac.f90 deleted file mode 100755 index e256fdd..0000000 --- a/OSBRACKETS/allosbrac.f90 +++ /dev/null @@ -1,451 +0,0 @@ - SUBROUTINE ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC) - -! COMMENTS ON USE: -!--------------------------------------------------------------------------------------------------------------------------------------- -! THE OSCILLATOR BRACKETS _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 - diff --git a/OSBRACKETS/osbrac.f90 b/OSBRACKETS/osbrac.f90 deleted file mode 100755 index 6bc1c41..0000000 --- a/OSBRACKETS/osbrac.f90 +++ /dev/null @@ -1,467 +0,0 @@ - SUBROUTINE OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC) - -! COMMENTS ON USE: -!--------------------------------------------------------------------------------------------------------------------------------------- -! THE OSCILLATOR BRACKETS _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 - diff --git a/OSBRACKETS/testallosbrac.f90 b/OSBRACKETS/testallosbrac.f90 deleted file mode 100755 index 02c3761..0000000 --- a/OSBRACKETS/testallosbrac.f90 +++ /dev/null @@ -1,260 +0,0 @@ - program testallosbrac -! Tests of the program that calculates the _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 _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 _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 _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 -! _L^\varphi=_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=\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 = 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 =',s,' exact value (number of states) =', numb - write(13,*)'3. Test of the relation sum_{i,ip,L}sum_j = 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 =',s,' exact value (number of states)=',numb -! 4. Test of the relation \sum_{i,i',L}|\sum_j-\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-\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|_L_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|_L_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-delta_{i1,i2}| =',s - write(13,*)'4. Test of the relation sum_{i,ip,L}sum_j|_L_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-delta_{i1,i2}| =',s - end diff --git a/OSBRACKETS/testosbrac.f90 b/OSBRACKETS/testosbrac.f90 deleted file mode 100755 index c791cd4..0000000 --- a/OSBRACKETS/testosbrac.f90 +++ /dev/null @@ -1,266 +0,0 @@ - program testosbrac -! Tests of the program that calculates the _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 _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 _L^phi' - write(6,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l - write(13,*)'1. bracket _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 -! _L^\varphi=_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_L_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}_L_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}_L_L =',s,' exact value (number of states) =',numb - write(13,*)'3. Test of the relation sum_{i,ip,L}sum_j = 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}_L_L =',s,' exact value (number of states) =',numb -! 4. Test of the relation \sum_{i,i',L}\sum_j|_L_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-\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|_L_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|_L_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_L_L-delta_{i1,i2}| =',s - write(13,*)'4. Test of the relation sum_{i,ip,L}sum_j|_L_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_L_L-delta_{i1,i2}| =',s - end