C ARTURO QUIRANTES SIERRA C Department of Applied Physics, Faculty of Sciences C University of Granada, 18071 Granada (SPAIN) C http://www.ugr.es/local/aquiran/codigos.htm C aquiran@ugr.es C c Last update: 10 December 2004 C C Subroutine APU C to calculate sets of Clebsch-Gordan coefficients C C First of all, a warning: C This subroutine has been taylor-made for other programs C That's why it doesn't calculate one single C-G coefficient C For single C-G calculations, I recommend subroutine NED C (also available at the above URL) C C APU calculates a set of CG by using recurrende relations. C Specifically, we will calculate C-G(aj,bj,cj;am,bm,am+bm) C For given values of aj,bj,cj,bm C and am restricted to IAM1<=am<=IAM2 C Reference: Schulten y Gordon, J. Math. Phys. 16, 1961-1970 (1975) C You can also find interesting the paper of C Wielaard, Mishchenko, Macke, and Carlson: Appl. Opt. 36, 4305-4331 (1997) C C You need to add a "APU(IAJ,IBJ,ICJ,IAM1,IAM2,IBM,CGN)" in your main routine C Input: C IAJ=aj, IBJ=bj, ICJ=cj; IAM1,IAM2=range of am values; IBM=bm C Output: C CGN(n+1)=C-G(IAJ,IBJ,ICJ;n+IAM1,IBM,IAM+IBM) C n=0,1, ... ,IAM2-IAM1 C SUBROUTINE APU(IAJ,IBJ,ICJ,IAM1,IAM2,IBM,CGN) double precision cc1,cc2,dd1,dd2,cgn(500),CX(500),C1,SC,CMAX,CMAX2 INTEGER IAJ,IBJ,ICJ,IBM,IAM1,IAM2,IAM,IAMEDIO INTEGER DN1,KP,KP2,IAMA,IAMB,DNTOPE,IAMIN,IAMAX C First, let's check that icj and ibm falls within limits if(icj.lt.abs(iaj-ibj).or.icj.gt.(iaj+ibj)) go to 100 if(abs(ibm).gt.ibj) go to 200 C IAMIN, IAMAX are the lowest and highest possible values for am IAMIN=-MIN(IAJ,ICJ+IBM) IAMAX=MIN(IAJ,ICJ-IBM) C IAMIN, IAMAX might be such that CG=0 C If so happens, let's jump up till CG<>0 KP=MAX(0,IAMIN-IAM1) KP2=MAX(0,IAM2-IAMAX) DNTOPE=10 C Calculation of every possible CG value IAMA=IAMIN IAMB=IAMAX DN1=IAMAX-IAMIN CX(1)=1.0D0 CMAX=CX(1) SC=CX(1)*CX(1) if(dn1.eq.0) go to 100 dd2=dfloat(iaj*(iaj+1)-ibj*(ibj+1)+icj*(icj+1)) cc2=dfloat((iaj-iama+1)*(iaj+iama)) cc2=cc2*dfloat((icj-iama-ibm+1)*(icj+iama+ibm)) cc2=-sqrt(cc2) IF(DN1.GE.DNTOPE) IAMB=IAMIN+INT(DN1/2)-1 C If there are not many CG to calculate, upward recurrence alone will do iam=iama 5 iam=iam+1 icm=iam+ibm if(icm.gt.icj) go to 8 cc1=cc2 cc2=dfloat((iaj-iam+1)*(iaj+iam)) cc2=cc2*dfloat((icj-icm+1)*(icj+icm)) cc2=-sqrt(cc2) if(-icm.gt.icj) go to 8 dd1=dd2-2.0d0*dfloat((iam-1)*(icm-1)) CX(IAM-IAMIN+1)=-CX(IAM-IAMIN)*DD1/CC2 if(dn1.eq.1) then SC=SC+CX(IAM-IAMIN+1)*CX(IAM-IAMIN+1) go to 100 end if if(iam.eq.(iama+1)) go to 7 CX(IAM-IAMIN+1)=CX(IAM-IAMIN+1)-CX(IAM-IAMIN-1)*CC1/CC2 7 SC=SC+CX(IAM-IAMIN+1)*CX(IAM-IAMIN+1) CMAX2=ABS(CX(IAM-IAMIN+1)) IF(CMAX2.GT.CMAX) CMAX=CMAX2 C Now let's see if we need go on with recurrence if(iam.lt.iamb) go to 5 c CX(middle) is used for re-scaling. If equals zero, then disaster! c So let us just check and, if zero, go for the next one. c Due to round-off errors, we have to be a bit more tolerant: X<10^-14 c (X being the smallest/largest C-G ratio in our series) c The problem is, for high aj, bj... values, CG can have very small values c But then you would need extended precision anyway, so you might change c the tolerance criterion to, say, (c2max/cmax)<10¯30 if(dn1.ge.dntope.and.iam.ge.iamb.and.(cmax2/cmax).le.(1.0D-14)) go to 5 8 continue if(dn1.lt.dntope) go to 100 iamedio=iam C1=CX(IAMEDIO-IAMIN+1) SC=SC-C1*C1 C Downward recurrence IAMA=IAMEDIO iamb=iamax CX(IAMB-IAMIN+1)=1.0D0 cc1=dfloat((iaj-iamb)*(iaj+iamb+1)) cc1=cc1*dfloat((icj-iamb-ibm)*(icj+iamb+ibm+1)) cc1=-sqrt(cc1) do iam=iamb-1,iama,-1 icm=iam+ibm if(icm.gt.icj) go to 100 cc2=cc1 cc1=dfloat((iaj-iam)*(iaj+iam+1)) cc1=cc1*dfloat((icj-icm)*(icj+icm+1)) cc1=-dsqrt(cc1) if(-icm.gt.icj) go to 18 dd1=dd2-2.0d0*dfloat((iam+1)*(icm+1)) CX(IAM-IAMIN+1)=-CX(IAM-IAMIN+2)*DD1/CC1 if(iam.eq.(iamb-1)) go to 18 CX(IAM-IAMIN+1)=CX(IAM-IAMIN+1)-CX(IAM-IAMIN+3)*CC2/CC1 18 end do C1=C1/CX(IAMEDIO-IAMIN+1) C Now, rescaling do iam=IAMEDIO,IAMAX cx(iam-iamin+1)=cx(iam-iamin+1)*c1 sc=sc+cx(iam-iamin+1)*cx(iam-iamin+1) end do C Here's where the combined (upward+downward) recurrence ends 100 CONTINUE SC=(dfloat(2*ICJ+1)/dfloat(2*IBJ+1))/SC SC=SQRT(SC) C Now, let's set the sign for SC if(cx(iamax-iamin+1).ne.abs(cx(iamax-iamin+1))) SC=-SC if(mod((iaj+iamax),2).ne.0) SC=-SC if(kp.ne.0) then do iam=0,kp-1 cgn(iam+1)=0.0d0 end do end if if(kp2.ne.0) then do iam=iam2-kp2+1,iam2 cgn(iam-iam1+1)=0.0d0 end do end if do iam=iam1+kp,iam2-kp2 cgn(iam-iam1+1)=cx(iam-iamin+1)*SC end do 200 CONTINUE RETURN END