By Reinaldo Baretti Machin
e-mail : reibaretti2004@yahoo.com
URL :  www.oocities.org/reibaretti2004
       www.oocities.org/serienumerica

c coulomb potential in p-space
      g(p)=anorm/(p**2+alfa**2)**2
      Q(u)= -.5*log((u-1.)/(u+1.))
      pi=2.*asin(1.)
      z=1.
      alfa=z
      c=137.11
      edirac=c**2*(sqrt(1.0-(z/c)**2)-1.0)
      nstep1=2000
      pf=60.*z
      dp1=pf/float(nstep1)
      dp2=dp1
      nstep2=int(pf/dp2)
      call  consta(z,alfa,nstep1,dp1,pi,anorm)
      call akin(alfa,c,nstep1,dp1,pi,anorm,ecin)
      sum=0.
      do 10 j=1,nstep2
      p2=dp2*float(j)
      do 10 i=1,nstep1
      p1=dp1*float(i)
      if(p1.ne.p2)then
      u=(p1**2+p2**2)/(2.*p1*p2)
      sum=sum+(p1-dp1/2.)*(p2-dp2/2.)*q(u)*g(p1-dp1/2.)*
     $ g(p2-dp2/2.)
      sumres=(p1-dp1/2.)*(p2-dp2/2.)*q(u)*g(p1-dp1/2.)*
     $ g(p2-dp2/2.)
      endif
      if(p1.eq.p2)then
      sum=sum+sumres
      endif
10    continue
      sum=sum*(-4.*z)*dp1*dp2
      vc=sum
      et=ecin+vc
c      print*,'edirac=',edirac
      eh=-z**2/2.
      print*,'eh=',eh
      print 100, ecin, vc, et
100   format(2x,'Ecin, Vc, Et=',3(2x,e10.3))
      stop
      end
      
      subroutine consta(z,alfa,nstep1,dp1,pi,anorm)
      psi(p)=1.0/(p**2+alfa**2)**2
      sum=0.
      do 10 i=1,nstep1
      p=dp1*float(i)
      sum=sum + (Psi(p)*p)**2+(psi(p-dp1)*(p-dp1))**2
10    continue
      sum=sum*4.*pi*dp1/2.
      anorm=sqrt(1./sum)
      return
      end
      
      subroutine  akin(alfa,c,nstep1,dp1,pi,anorm,ecin)
      psi(p)=anorm/(p**2+alfa**2)**2
c      tk(p)=0.5*p**2 -(p**4)/(8.*c**2)
      tk(p)=0.5*p**2
      sum=0.
      do 10 i=1,nstep1
      p=dp1*float(i)
      sum=sum + (p*psi(p))**2*tk(p)+ ((p-dp1)*psi(p-dp1))**2*tk(p-dp1)
10    continue
      sum=sum*4*pi*dp1/2.
      ecin=sum
      return
      end
      
      Run
       eh= -0.5
      Ecin, Vc, Et=   0.500E+00  -0.993E+00  -0.493E+00