c coulsubs.f, version 2.0, November 1, 1995. subroutine coulset(kk,mom,eta) implicit none include 'cparams.f' double precision mom,eta integer i0,kk double complex acouly(nka),acfact1(nka),acfact2(nka), + achype(40,nka),achype1(5,nka),achype2(5,nka) double complex couly,cfact1,cfact2,chype(40), + chype1(5),chype2(5) common /couljunk/acouly,acfact1,acfact2,achype,achype1,achype2 call coulpar(couly,cfact1,cfact2, + chype,chype1,chype2,eta) acouly(kk)=couly acfact1(kk)=cfact1 acfact2(kk)=cfact2 do 993 i0=1,40 achype(i0,kk)=chype(i0) 993 continue do 994 i0=1,5 achype1(i0,kk)=chype1(i0) achype2(i0,kk)=chype2(i0) 994 continue return end c ___________________________________________ double complex function cgamma(c) c This calc.s gamma functions which are in the form gamma(n+i*y) c where n is an integer and y is real. implicit none double complex c,ci,cg,cphase parameter (ci=(0.d0,1.d0)) integer mm,j double precision x,y,phase,euler,delp,cgmag,pi parameter (euler=.5772156649d0,pi=3.14159265358979323844d0) x=c y=-ci*c if(dabs(y).lt.1E-5) then cg=(1.d0,0.d0) go to 888 endif phase=-euler*y do 123 j=1,5000 delp=(y/dfloat(j))-atan(y/dfloat(j)) phase=phase+delp if(dabs(delp).lt.1E-5) go to 223 123 continue write(66,*) 'oops did not make it in cgamma' 223 phase=phase-2.d0*pi*int(phase/(2.d0*pi)) cphase=cdexp(ci*phase) cgmag=sqrt(pi*y/sinh(pi*y)) cg=cgmag*cphase 888 mm=nint(x) if(mm.lt.1) then do 444 j=1,-mm+1 cg=cg/(1.d0+dfloat(-j)+ci*y) 444 continue endif if(mm.gt.1) then do 333 j=1,mm-1 cg=cg*(dfloat(j)+ci*y) 333 continue endif cgamma=cg return end ****************************************** double complex function cgamma2(c) double complex cy,cf,c double precision rc,sign,drb,rb(8) integer j,mm data (rb(j), j=1,8)/1.d0,1.64493d0,1.20205d0,1.08232d0, + 1.03692d0,1.01734d0,1.00834d0,1.00407d0/ rc=c mm=int(.01d0+rc) cy=c cy=cy-mm cf=cy*(1.d0-.5772156649d0)-cdlog(1.d0+cy) sign=-1.d0 do 99 j=2,8 sign=-sign drb=rb(j) cf=cf + sign*(drb-1.d0)*(cy**j)/dfloat(j) 99 continue cf=cdexp(cf) if(mm.gt.0) then do 101 j=1,mm-1 cy=cy+dfloat(j) cf=cf*cy 101 continue else cf=cf/cy endif cgamma2=cf return end ***************************** subroutine coulpar(couly,cfact1,cfact2,chype,chype1, + chype2,eta) implicit none double complex a,b,a0,b0,ci,chype(40),chype1(5),chype2(5) double complex c1top1,c2top2,cgamma,cfact1,cfact2,bot double complex c1top2,c2top1,couly double precision eta integer j ci=(0.d0,1.d0) c See appendix of Messiah. "cgamma" is the gamma function. "chyper" is c appropriate hyperbolic function. This is explained in Messiah's c section on coul wave func.s. c Notation should be explanatory when reading the book. couly=cgamma(1.d0+ci*eta) couly=couly*dexp(-.5d0*eta*3.14159265358979323844d0) a=-ci*eta b=(1.d0,0.d0) a0=a b0=b do 175 j=1,40 chype(j)=a/(b*dfloat(j)) a=a+1.d0 b=b+1.d0 175 continue c1top1=1.d0 c2top1=1.d0 c1top2=1.d0 c2top2=1.d0 bot=1.d0 do 75 j=1,4 c1top1=c1top1*(dfloat(j)+a0-1.d0) c2top1=c2top1*(dfloat(j)-a0) c1top2=c1top2*(dfloat(j)-b0+a0) c2top2=c2top2*(dfloat(j)+b0-a0-1.d0) bot=bot*dfloat(j) chype1(j)=(c1top1*c1top2)/(bot) chype2(j)=(c2top1*c2top2)/(bot) 75 continue cfact1=(cgamma(b0))/(cgamma(b0-a0)) cfact2=(cgamma(b0))/(cgamma(a0)) return end ****************** double complex function coul(mom,r,zk,eta,kk) implicit none include 'cparams.f' integer kk double complex ci,bb,coulguy,chyper c double complex acouly(nka),acfact1(nka),acfact2(nka), + achype(40,nka),achype1(5,nka),achype2(5,nka) common /couljunk/acouly,acfact1,acfact2,achype,achype1,achype2 c external cgamma,chyper double precision mom,eta,arg,zk,r ci=(0.d0,1.d0) c See appendix of Messiah. "cgamma" is the gamma function. c "chyper" is c appropriate hyperbolic function. This is explained in Messiah's c section on coul wave func.s. c Notation should be explanatory when reading the book. bb=(1.d0,0.d0) coulguy=acouly(kk)*chyper(-ci*eta,bb,ci*mom*(r-zk)/197.323d0,kk) arg=zk*mom/197.323d0 arg=arg-2.d0*3.14159265358979323844d0 + *int(arg/(2.d0*3.14159265358979323844d0)) coulguy=coulguy*cdexp(ci*arg) coul=coulguy return end ********************************************************************** double complex function chyper(a,b,cz,kk) implicit none include 'cparams.f' double complex acouly(nka),acfact1(nka),acfact2(nka), + achype(40,nka),achype1(5,nka),achype2(5,nka) common /couljunk/acouly,acfact1,acfact2,achype,achype1,achype2 double complex cz,a,b,cw1,cw2,cf1,cf2,czarg,czstarj,cfact,ci double precision dmag,rcrit integer kk,j parameter(ci=(0.d0,1.d0),rcrit=10.d0) dmag=cdabs(cz) if(dmag.lt.rcrit) then cf1=(1.d0,0.d0) czstarj=(1.d0,0.d0) do 175 j=1,40 czstarj=czstarj*cz*achype(j,kk) cf1=cf1+czstarj if(cdabs(czstarj).lt.0.0001d0) go to 275 175 continue write(66,*) 'Holy cow batman I made it here.' 275 chyper=cf1 else cf1=1.d0 cf2=1.d0 do 75 j=1,5 cf1=cf1+achype1(j,kk)/((-cz)**j) cf2=cf2+achype2(j,kk)/((cz)**j) 75 continue cfact=acfact1(kk)*((-cz)**(-a)) cw1=cfact*cf1 czarg=cz-ci*2.d0*3.14159265358979323844d0 + *int(-ci*cz/(2.d0*3.14159265358979323844d0)) cfact=acfact2(kk)*((cz)**(a-b))*cdexp(czarg) cw2=cfact*cf2 chyper=cw1+cw2 endif return end