	SUBROUTINE SUMRUL
	IMPLICIT REAL*8(A-H,O-Z)
C
C    THIS SUBROUTINE CALCULATES THE ISOSCALAR SUM RULES FOR THE
C    MONOPOLE,DIPOLE AND HIGHER MULTIPOLE OPERATORS.
C
	COMMON AMASS,Z,A13,EPSILON,AR(3),R(3),SR(3),L_MULT
	DIMENSION R2(3),R4(3),RLN(3)
C   
C    GENERAL PARAMETERS
C
	PARAMETER (TWOPI = 6.283185D0)
	PARAMETER (HBARC2 = 38927.29D0)
	PARAMETER (MC2 = 938.256D0)
C
C
	Z2A=Z*Z/AMASS
C
C    GO HERE FOR MONOPOLE TRANSITIONS
C
	IF(L_MULT.EQ.0)THEN
	TEMP=HBARC2*Z2A/(2*MC2)
	DO 1 I=1,3
   	CALL RFERMI(R(I),AR(I),2,R2(I))
   1	SR(I)=TEMP*R2(I)
C
C    GO HERE FOR DIPOLE TRANSITIONS
C
	ELSE IF(L_MULT.EQ.1)THEN
	EGQR=65./A13
	EGMR=80./A13
	EPSILON=(4./EGQR+5./EGMR)*HBARC2/(3.*MC2*AMASS)
	TEMP=3.*HBARC2*Z2A/(16.*MC2*TWOPI)
	DO 2 I=1,3
   	CALL RFERMI(R(I),AR(I),2,R2(I))
   	CALL RFERMI(R(I),AR(I),4,R4(I))
   2	SR(I)=TEMP*(11.*R4(I)-(25./3.)*R2(I)*R2(I)-10.*EPSILON*R2(I))
C
C    GO HERE FOR L_MULT.GE.2
C
	ELSE IF(L_MULT.GE.2)THEN
	XL=FLOAT(L_MULT)
	LN=2*L_MULT-2
	TEMP=XL*(2.*XL+1.)*(2.*XL+1.)*HBARC2*Z2A/(4.*TWOPI*MC2)
	DO 3 I=1,3
   	CALL RFERMI(R(I),AR(I),LN,RLN(I))
   3	SR(I)=TEMP*RLN(I)
	ENDIF
	RETURN
	END
  
	SUBROUTINE RFERMI(RR,AR,N,RN)
	IMPLICIT REAL*8(A-H,O-Z)
C
C    THIS SUBROUTINE CALCULATES THE MEAN TO THE Lth POWER OF THE RADIUS 
C    FOR A FERMI OR A UNIFORM MASS DISTRIBUTION {i.e.<R**L>}.
C
	CHARACTER DEF*1
C	Changed by JE 14.04.2010
C	COMMON /DEFSR/ DEF
	COMMON /INT/ DR,KR
	DIMENSION U(500)
	DATA DEF /'N'/
	IF (AR.EQ.0.) THEN
  
C   UNIFORM DISTRIBUTION
  
	  RN=(3./(3.+FLOAT(N)))*RR**N
	ELSE
  
C  FERMI DISTRIBUTION
  
	  R=0.D0
	  SUM1=0.D0
	  SUM2=0.D0
	  IF(DEF.EQ.'Y') THEN
	    CALL POTENT(U,'GR',0,1,RR,AR)
	  ELSE
	    F2=EXP(-DR/AR)
	    F1=EXP(RR/AR)
	  ENDIF
	  DO 1 M=1,KR
	  R=R+DR	
	  RPWR=R**(N+2)
	  R2=R*R
	  X=1.D0
	  IF(M.EQ.1.OR.M.EQ.KR)X=0.5D0
	  IF(DEF.EQ.'Y') THEN
	    TEMP=X*U(M)
	  ELSE
	    F1=F1*F2
	    TEMP=X*F1/(1.0+F1)
	  ENDIF
	  SUM1=SUM1+TEMP*R2
	  SUM2=SUM2+TEMP*RPWR
  1  	  CONTINUE
	  RN=SUM2/SUM1
	ENDIF
	RETURN
	END
  
	SUBROUTINE INTEG(BAND,L_MULT,NX,RR,AR,MEL,MELCOR)
	IMPLICIT REAL*8(A-H,O-Z)
C
C   SUBROUTINE TO CALCULATE TRANSITION MATRIX ELEMENTS FROM THE GROUND 
C   STATE TO STATES OF THE GROUND STATE, GAMMA, BETA AND OCTUPOLE BANDS
C
	COMMON /DEF/ VBETA,GAMMA,BETACOR,XK,BETA0_DEF,BETA3_DEF 
	1           ,BETA(3),LDFRM(3)
	COMMON /INT/ DR,KR
	DIMENSION U(500),V(500)
	CHARACTER BAND*2
	REAL*8 MEL,MELCOR
	PARAMETER (FOURPI=12.56637062D0)
C
C   GROUND STATE IS ASSUMED TO HAVE SPIN AND PARITY  0+
C
	XI = 0.D0
	XJ = L_MULT
	XKI = 0.D0
	XKJ=XK
	XLAM = XJ
  
C   THE SIGN SHOULD BE TAKEN POSITIVE TO CONFIRM WITH THE OTHER CONVENTIONS
  
	TEMP=+SQRT((2.*XJ+1.)/(2.*XI+1.))
	IF (BAND.EQ.'GA') TEMP=TEMP*GAMMA*BETA(1)
	IF (BAND.EQ.'BE'.OR.BAND.EQ.'OC') TEMP=TEMP*VBETA
  
	IF(BAND.EQ.'GR')THEN
C
C   GROUND STATE ROTATIONAL TRANSITION DENSITY
C
	  IF(XKJ.NE.XKI) TEMP=0.D0
	  TEMP=TEMP*CG(XJ,XKJ,XLAM,0.D0,XI,XKI)*PAR(NINT(XJ-XI))
C
C   GAMMA, BETA OR OCTUPOLE BAND TRANSITION DENSITY
C
	ELSE 
	  TEMP=TEMP*CG(XJ,-XKJ,XLAM,XKJ,XI,0.D0)*PAR(L_MULT)
	ENDIF
  
	CALL POTENT(U,BAND,L_MULT,NX,RR,AR)
	INX=NX
	IF(BAND.NE.'GR') INX=NX-1
	CALL POTENT(V,'GR',0,INX,RR,AR)
	L = L_MULT + 2
	IF(L_MULT.EQ.0) L=4
	IF(L_MULT.EQ.1) L=5
	R =0.D0
	MEL = 0.D0
	MELCOR = 0.D0
	SUM=0.D0
	DO 2 I=1,KR
	R=R+DR
	RPWR= R**L
	R2=R*R
	X=1.D0
	IF(I.EQ.1.OR.I.EQ.KR) X=0.5D0
	MEL=MEL+X*RPWR*U(I)
	IF(L_MULT.GT.1) GO TO 3
	MELCOR=MELCOR+X*RPWR*U(I)/R2
  3	SUM=SUM+X*R2*V(I)
  2	CONTINUE
	MEL=TEMP*MEL/(FOURPI*SUM)
	MELCOR=TEMP*MELCOR/(FOURPI*SUM)
	IF(L_MULT.EQ.0) THEN
  
C   THE MONOPOLE TRANSITION IS CALCULATED FOR THE OPERATOR (1/2)*R**2
  
	MEL=MEL*SQRT(FOURPI)/2.
	MELCOR=MELCOR*1.772454
	ELSE IF(L_MULT.EQ.1) THEN
  
C   THE DIPOLE TRANSITION IS CALCULATED FOR THE OPERATOR (1/2)*R**3*Y1
  
	MEL=MEL*0.5
	MELCOR=MELCOR*0.5
	ENDIF
	RETURN
	END
  
	SUBROUTINE POTENT(U,BAND,LAM,NX,RR,AR)
	IMPLICIT REAL*8(A-H,O-Z)
C
C    DEFORMED VOLUME,SURFACE OR SECOND DERIVATIVE POTENTIAL BY YLM
C    EXPANSION
C
	COMMON /DEF/ VBETA,GAMMA,BETACOR,XK,BETA0_DEF,BETA3_DEF
	1           ,BETA(3),LDFRM(3)
	COMMON /INT/ DR,KR
	DIMENSION U(500)
	DIMENSION CN(16),YLAM(16),CP(16),FACTR(3),X(8),W(8)
	CHARACTER BAND*2
	DATA X/0.0950125098D0,0.2816035507D0,0.4580167776D0,0.6178762444
	1D0,0.7554044083D0,0.8656312023D0,0.9445750230D0,0.9894009439D0/
	DATA W/0.1894506104D0,0.1826034150D0,0.1691565193D0,0.1495959888
	1D0,0.1246289712D0,0.0951585116D0,0.0622535239D0,0.0271524594D0/
	DATA SQRPI/1.772453851D0/
C
	IND=0
	DO 6 I=1,KR
  6	U(I) = 0.D0
	XLAM=LAM
  
C  FOR UNIFORM DENSITY AR MUST BE FINITE
	IF (AR.EQ.0.) AR=0.03*RR/1.2
  
	F2=EXP(-DR/AR)
	VR = 1.D0
	IF(BAND.NE.'GR') VR=RR/AR
	GO TO 11
  10	CONTINUE
	IND=1
	VR=VR*BETACOR
  11	CONTINUE
	IF(BAND.EQ.'GR'.OR.BAND.EQ.'GA') GO TO 12
	IF(IND.EQ.0.AND.(BETA0_DEF.EQ.1..OR.BETA3_DEF.EQ.1.)) GO TO 70
  12	CONTINUE
     	DO 1 I=1,3
	FACTR(I)=SQRT(DFLOAT(LDFRM(I)+LDFRM(I)+1))*BETA(I)/(SQRPI*2.)
  1	CONTINUE
	LMAX=MAX(LAM,LDFRM(3))+1
	LTEMP=2
 	IF(BAND.EQ.'OC') LTEMP=3
	IF(IND.EQ.1.AND.BAND.EQ.'BE') LTEMP=0
	IF(IND.EQ.1.AND.BAND.EQ.'OC') LTEMP=1
	XTEMP=LTEMP
	IF(BAND.NE.'GR') LMAX=MAX(LAM+LTEMP,LDFRM(3))+1
	DO 3 I=1,8
	YLAM(I  )=0.D0
  3	YLAM(I+8)=0.D0
	T2=(-1.)**LAM
	T1=(-1.)**LDFRM(1)
	DO 9 I=1,8
	CN(I  )=0.D0
	CN(I+8)=0.D0
	P2=0.D0
	P3=1.D0
	Z=X(I)
	DO 8 M=1,LMAX
	L=M-1
	XL = L
	FL=XL-1.
	IF(L.EQ.0)GOTO 2
	P3=((2.*FL+1.)*Z*P2-FL*P1)/(FL+1.)
  2	CONTINUE
	DO 4 J=1,3
	IF(L.NE.LDFRM(J)) GO TO 4
	CN(I  )=CN(I  )+P3*FACTR(J)
	CN(I+8)=CN(I+8)+P3*FACTR(J)*T1
  4	CONTINUE
	IF(BAND.EQ.'GR') GO TO 5
	IF(IABS(L-LAM).GT.LTEMP.OR.L+LAM.LT.LTEMP) GOTO 7
	T3=SQRT((2.*XTEMP+1.)*(2.*XL+1.)/(2.*XLAM+1.))*0.5/SQRPI
	1*CG(XL,0.D0,XTEMP,0.D0,XLAM,0.D0)*CG(XL,0.D0,XTEMP,XK,XLAM,XK)
	2 *P3*W(I)*SQRT(2.*XL+1.)*SQRPI
	YLAM(I  )=YLAM(I  )+T3
	YLAM(I+8)=YLAM(I+8)+T3
	GOTO 7
  5	IF(L.NE.LAM)GOTO 7
	YLAM(I  )=P3*W(I)*SQRT(2.*XL+1.)*SQRPI
	YLAM(I+8)=YLAM(I  )*T2
  7	CONTINUE
	P1=P2
	P2=P3
  8	CONTINUE
  9	CONTINUE
	F1=1.D0
	DO 35 J=1,16
  35	CP(J)=EXP((1.+CN(J))*RR/AR)
	IF(LAM.EQ.0.AND.BAND.EQ.'GR') VR=VR/(2.*SQRPI)
	DO 60 M=1,KR
	VTEMP=0.D0
	F1=F1*F2
	DO 55 J=1,16
	F3=F1*CP(J)
	GO TO (40,45,50),NX
  40	VTEMP=VTEMP+YLAM(J)*VR*F3/(1.+F3)
	GO TO 55
  45	VTEMP=VTEMP-YLAM(J)*VR*(F3/(1.+F3))/(1.+F3)
	GO TO 55
  50	VTEMP=VTEMP+YLAM(J)*VR*(F3/(1.+F3))*((1.-F3)/(1.+F3))/(1.+F3)
  55	CONTINUE
	U(M)=U(M)+VTEMP
  60	CONTINUE
  70	CONTINUE
	IF(BAND.EQ.'GR'.OR.BAND.EQ.'GA') RETURN
	IF(IND.EQ.0) GO TO 10
	RETURN
	END
