C-----CAMATREL----------------MAIN----------------------C.E THEODOSIOU--
C 18-APR-1982 *** adapted for VAX 11.apr.1989
C +----------+ +-----------------+ +----------------------------------+
C | CAMATREL | | C.E. THEODOSIOU | | VERSION OF 14.AP.1983 / AT 21:00 |
C +----------+ +-----------------+ +----------------------------------+
C +-------------------------------------------------------------------+
C | |
C | THIS PROGRAM CALCULATES RADIAL MATRIX ELEMENTS FOR MULTIPOLE |
C | TRANSITIONS BETWEEN TWO RYDBERG-STATE ORBITALS: |
C | |
C | R = ( N1*(N1),L1 |R**IP| N2*(N2),L2 ) |
C | |
C +-------------------------------------------------------------------+
C | REQUIRED SUBPROGRAMS : SUBROUTINES: NONE REQUIRED |
C | FUNCTIONS : RADLME, GMMMA |
C +-------------------------------------------------------------------+
C
C EXPLANATION OF PARAMETERS AND SYMBOLS
C -------------------------------------
C
C ID : LABELING OF THE RUN. PROVIDES NONCRUCIAL INFORMATION
C CHARGE : EFFECTIVE CHARGE OF THE ATOM OR ION
C XN1 : N1*, ONE OF THE TWO EFFECTINE N'S
C XN2 : N2*, THE OTHER OF THE TWO N*'S
C IP : ORDER OF THE TRANSITION. I.E.:
C IP=1 --> DIPOLE,
C IP=2 --> QUADRUPOLE, ETC.
C L1 : ANGULAR MOMENTUM VALUE OF N1* STATE
C L2 : ANGULAR MOMENTUM VALUE OF N2* STATE
C
C-----------------------------------------------------------------------
C
C THE SUBPROGRAMS WERE DEVELOPED BY THE COLUMBIA U. GROUP.
C THEY ARE TAKEN FROM THE PROGRAM LISTING SENT BY W. HAPPER TO CET.
C
C-----------------------------------------------------------------------
C
IMPLICIT REAL*8 (A-H,O-X,Z)
CHARACTER*4 ENDATA,RDOLD,IDATA,IEND,ID
C
DIMENSION JL(9),ID(20)
C
DATA JL/1HS,1HP,1HD,1HF,1HG,1HI,1HJ,1HK,1HL/
DATA ENDATA/'&&&&'/,RDOLD/'- '/,IDATA/'DATA'/,IEND/'END '/
CHARGE=1.D0
ISKIP=0
c
write(6,*) ' EXPLANATION OF PARAMETERS AND SYMBOLS'
write(6,*) '----------------------------------------------------'
write(6,*) ' ID : LABELING OF THE RUN. (NONCRUCIAL)'
write(6,*) ' CHARGE : EFFECTIVE CHARGE OF THE ATOM OR ION'
write(6,*) ' XN1 : N1*, ONE OF THE TWO EFFECTINE N* s'
write(6,*) ' XN2 : N2*, THE OTHER OF THE TWO N* s'
write(6,*) ' IP : ORDER OF THE TRANSITION. I.E.:'
write(6,*) ' IP=1 : DIPOLE, IP=2 : QUADRUPOLE, ETC.'
write(6,*) ' L1 : ANGULAR MOMENTUM VALUE OF N1* STATE'
write(6,*) ' L2 : ANGULAR MOMENTUM VALUE OF N2* STATE'
write(6,*) '----------------------------------------------------'
write(6,*)
C write(6,*) ' what is the output unit?'
C read(5,*) iw
iw=1
c
C open (unit=iw,file='camatrel.out',status='new')
c
C WRITE(iw,5000)

5000 FORMAT (1H1)
C
1 write(6,*) ' what is the new ID ?'
READ (5,1000,END=999) ID
1000 FORMAT (20A4)
ISKIP=0
IF (ID(1).EQ.ENDATA) STOP
C WRITE(iw,2000) ID

2000 FORMAT (1H1,2(/),5X,20A4,/1X,74(1H-))
2 write(6,*) ' what is the CHARGE of the ion ?'
READ (5,*) CHARGE
c 1010 FORMAT (2F8.5,6I5)

3 write(6,*) ' what are "xn1,xn2,ip,l1,l2,n1,n2" ?'
READ (5,*,END=999) XN1,XN2,IP,L1,L2,N1,N2
IF (XN1.LT.0) GO TO 1
ISKIP=1
C
C..........................CALCULATION OF THE MATRIX ELEMENT :
C..........................(N1*,L1 ; N2*,L2),WHERE N1* WAS CALCULATED IN
JL1ID=JL(L1+1)
JL2ID=JL(L2+1)
F1= RADLME (IP,L1,L2,XN1,XN2,CHARGE)
WRITE(6,2010) N1,JL1ID,IP,N2,JL2ID,F1,XN1,XN2
C WRITE(iw,2010) N1,JL1ID,IP,N2,JL2ID,F1,XN1,XN2

2010 FORMAT (/3X,'(',I2,A1,'|R**',I1,'|',I2,A1,') =',F7.4,5X
* ,'N1*=',F6.3,' N2*=',F6.3)
GO TO 3
C-----------------------------------------------------------------------
999 CONTINUE
STOP
END
c.cp

C
C-----BLOCK DATA----------SUBROUTINE-------------------C.E.THEODOSIOU---
C 24.DEC.1985
BLOCK DATA
C
C COMPILATION OF ALL DATA NECESSARY FOR MAIN PROGRAM AND SUBS
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /NUMBRS/ D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,HALF,QUART,D12
COMMON /CRYDBR/ RYDBRG,RYDINF,RYDRAT
C
DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,HALF,QUART,D12
X /0.D0,1.D0,2.D0,3.D0,4.D0,5.D0,6.D0,7.D0,8.D0,9.D0,10.D0
2 ,0.5D0,0.25D0,12.D0/
DATA RYDINF/109737.31529D+00/
C
C END
C
C-----ERRSET/DUMMY-------SUBROUTINE-------------------------------------
C
C SUBROUTINE ERRSET(I,J,K,L)
C
C RETURN
END
c.cp

C
C-----GMMMA--------------FUNCTION-------------------------MAR. 19,1982--
C
SUBROUTINE GMMMA (XX,GX,IER)
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 XX,GX
C
C01=+1.0D0
C02=+5.772156649015329D-01
C03=-6.558780715202538D-01
C04=-4.20026350340952D-02
C05=+1.665386113822915D-01
C06=-4.21977345555443D-02
C07=-9.621971527877D-03
C08=+7.218943246663D-03
C09=-1.1651675918591D-03
C10=-2.152416741149D-04
C11=+1.280502823882D-04
C12=-2.01348547807D-05
C13=-1.2504934821D-06
C14=+1.133027232D-06
C15=-2.056338417D-07
C16=+6.116095D-09
C17=+5.0020075D-09
C18=-1.1812746D-09
C19=+1.043427D-10
C20=+7.7823D-12
C21=-3.6968D-12
C22=+5.1D-13
C23=-2.06D-14
C24=-5.4D-15
C25=+1.4D-15
C26=+1.0D-16
IF(XX-57.D0)6,6,4
4 IER=2
GX=1.D36
C GINV=1.D0/GX
RETURN
6 X=XX
ERR=1.0D-6
IER=0
GX=1.0D0
IF(X-2.0D0)50,50,15
10 IF(X-2.0D0)110,110,15
15 X=X-1.0D0
GX=GX*X
GO TO 10
50 IF(X-1.0D0)60,120,110
C
C--------------------------- SEE IF X IS NEAR NEGATIVE INTEGER OR ZERO
C
60 IF(X-ERR)62,62,80
62 Y=DFLOAT(IDINT(X))-X
IF(DABS(Y)-ERR)130,130,64
64 IF(1.0D0-Y-ERR)130,130,70
C
C------------------------------- X NOT NEAR A NEGATIVE INTEGER OR ZERO
C
70 IF(X-1.0D0)80,80,110
80 GX=GX/X
X=X+1.0D0
GO TO 70
110 Y=X-1.0D0
31 TMP=C19+Y*(C20+Y*(C21+Y*(C22+Y*(C23+Y*(C24+Y*(C25+Y*C26))))))
TMP=C12+Y*(C13+Y*(C14+Y*(C15+Y*(C16+Y*(C17+Y*(C18+Y*TMP))))))
TMP=C05+Y*(C06+Y*(C07+Y*(C08+Y*(C09+Y*(C10+Y*(C11+Y*TMP))))))
GY=1.D0/(C01+Y*(C02+Y*(C03+Y*(C04+Y*TMP))))
GX=GX*GY
120 RETURN
130 IER=1
RETURN
END
c.cp

C
C-----RADLME-----------------FUNCTION------------------C.E. THEODOSIOU--
C 18-APR-1982
FUNCTION RADLME(IP,LG,LE,NG,NE,Q)
C
C ADAPTED FROM W. HAPPER'S PROGRAM LISTING.
C NG,NE : EFFECTIVE QUANTUM NUMBERS OF GROUND AND EXCITED
C STATES
C LG,LE : ANGULAR MOMENTA OF GROUND AND EXCITED STATES
C AG,AE : A-COEFFICIENTS OF BATES AND DAMGAARD
C CG,CE : C-COEFFICIENTS OF BATES AND DAMGAARD
C Q : IONICITY (=Z-N+1)
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 RADLME,NG,NE,NG0,Q
C
DIMENSION AG(50),AE(50),CG(50),CE(50)
DATA NG0/0.D0/,LG0/9999/
C
C------------ CALCULATE THE A AND C COEFFICIENTS FOR THE GROUND STATE --
C
IF(NG.EQ.NG0.AND.LG.EQ.LG0) GO TO 20
NG0=NG
LG0=LG
AG(1)=1.D0
L1=1
CALL GMMMA(NG+LG+1.D0,G1,IER)
CALL GMMMA(NG-LG,G2,IER)
c G1=1.D0/GAMINV(NG+LG+1.D0)

c G2=1.D0/GAMINV(NG-LG)

CG(L1)=(2.D0/NG)**NG*AG(L1)/(NG*DSQRT(G1*G2))
C
C-------- MAXIMUM INDEX FOR A AND C COEFFICIENTS IN THE GROUND STATE --
C
I1=NG+1.D0
DO 10 L1=2,I1
AG(L1)=AG(L1-1)*(NG/(2.D0*(L1-1.D0)*Q))*(LG*(LG+1.D0)
+ -(NG-L1+1.D0)*(NG-L1+2.D0))
CG(L1)=(2.D0/NG)**NG*Q**(L1-1)*AG(L1)/(NG*DSQRT(G1*G2))
10 CONTINUE
C
C------------------ CHECK THE NORMALIZATION OF THE RADIAL WAVEFUNCTIONS
C
CT SN(K)=0.0D0
CT DO 85 L=1,I1
CT DO 86 J=1,I1
CTC CALL GMMMA(NG+NG-L-J+3.D0,G5,IER)
CT G5=GAMINV(NG+NG-L-J+3.D0)
CT SN(K)=SN(K)+CG(L)*CG(J)*(NG*NG/(NG+NG))**(NG+NG-L-J+3.D0)*G5/Q
CT 86 CONTINUE
CT 85 CONTINUE
C
C--------------- CALCULATION DONE FOR ELECTRIC MULTIPOLE TRANSITIONS
C
CT IF(DABS(DABS(LE-LG)-1.D0).LT.1.D-05.AND.DABS(JE-JG).LT.1.0005D0)
CT 1 GO TO 20
CT GO TO 21
C
C--- CALCULATE THE A AND C COEFFICIENTS FOR THE EXCITED (HIGHER) STATES
C
20 AE(1)=1.D0
L2=1
CALL GMMMA(NE+LE+1.D0,G3,IER)
CALL GMMMA(NE-LE,G4,IER)
c G3=1.D0/GAMINV(NE+LE+1.D0)

c G4=1.D0/GAMINV(NE-LE)

CE(L2)=(2.0D0/NE)**NE*AE(L2)/(NE*DSQRT(G3*G4))
C
C--------- MAXIMUM INDEX FOR A AND C COEFFICIENTS IN THE HIGHER STATE
C
I2=NE+1.D0
DO 11 L2=2,I2
AE(L2)=AE(L2-1)*(NE/(2.D0*(L2-1.D0)*Q))*(LE*(LE+1.D0)
+ -(NE-L2+1.D0)*(NE-L2+2.D0))
CE(L2)=(2.D0/NE)**NE*Q**(L2-1)*AE(L2)/(NE*DSQRT(G3*G4))
11 CONTINUE
C
C-------------------------------------- CALCULATE THE RADIAL INTEGRAL
C
SUM=0.0D0
DO 12 L1=1,I1
DO 13 L2=1,I2
CALL GMMMA(NG+NE-L1-L2+IP+3.D0,G5,IER)
SUM=SUM+CG(L1)*CE(L2)*(NG*NE/(NG+NE))**(NG+NE-L1-L2+IP+3.D0)*G5
+ /Q**IP
13 CONTINUE
12 CONTINUE
RADLME=SUM
RETURN
END