GRAVITATIONAL POTENTIAL

c   gravitational potential calculations
c   observation point  P   (rp, thetap, phip)
      rho(r, theta, phi)=1.
      pi=2.*asin(1.)
      rp=2.
      thetap=pi/2.
      phip=0.
      r0=0.
      rlimit=1.
      theta0=0.
      theta2=pi
      phi0=0.
      phi2=2.*pi
      nstep=20
      dr=(rlimit-r0)/float(nstep)
      dtheta=(theta2-theta0)/float(nstep)
      dphi=(phi2-phi0)/float(nstep)
      sum=0.
      do 10 ir=1,nstep
      r=dr*float(ir)
      do 10 itheta=1,nstep
      theta=dtheta*float(itheta)
      do 10 iphi=1,nstep
      phi=dphi*float(iphi)
      sum=sum+rho(r-dr/2.,theta-dtheta/2.,phi-dphi/2.)*(r-dr/2.)**2*
     $sin(theta-dtheta/2.)
10    continue
      amass=sum*dr*dtheta*dphi
c    amtest massa for constant density rho=1.
      amtest=4.*pi/3.
      print*,'masa (num) , amtest=',amass,  amtest
c computes the potential at observation point (P)
      sum=0.
      do 20 ir=1,nstep
      r=dr*float(ir)
      do 20 itheta=1,nstep
      theta=dtheta*float(itheta)
      do 20 iphi=1,nstep
      phi=dphi*float(iphi)
      cosgam=sin(theta-dtheta/2.)*cos(phi-dphi/2.)*sin(thetap)*cos(phip)
     $+sin(theta-dtheta/2.)*sin(phi-dphi/2.)*sin(thetap)*sin(phip)
     $+cos(theta-dtheta/2.)*cos(thetap)
      denom=sqrt(rp**2+r**2-2.*rp*r*cosgam)
      sum=sum+rho(r-dr/2.,theta-dtheta/2.,phi-dphi/2.)*(r-dr/2.)**2*
     $sin(theta-dtheta/2.)/denom
20    continue
      sum=sum*dr*dtheta*dphi
      print*,'potnum, sphericalpot=', sum , amass/rp
      stop
      end