 program ThomasFermi
!Objective: calculate the Thomas-Fermi universal function (and a form factor)
 Implicit none
 integer i,N,jq
!Number of points on the grid
 parameter (N=50000)
 double precision xgrid(N),h,ymax,Phi(N),PhiDer(N),ygrid(N),mu,  &
                 rho_q(N),rgrid(N),a,a0,Z,pi,qgrid(N),alpha,me,qmax,hq,FF(N)


!maximum possible value of the y-variable (to be varied), dsqrt(x)
 ymax=dsqrt(1d1)
!the step in the y-variable
 h=ymax/dble(N-1)
!prepare the y-grid for integration
 do i=1,N
   ygrid(i)=h*dble(i-1)
!convert back to the x-variable
   xgrid(i)=ygrid(i)**2 
 enddo
!First part: solve the Thomas-Fermi equation
!initial conditions at the edge of the atom, the condition on the derivative 
!is meaningless, just to start running with something, to be rescaled later
 Phi(N)=1d-2  
 PhiDer(N)=-3d0*Phi(N)/ygrid(N)**2
!step backwards over the entire grid down to smallest x (almost zero)
!construct, by Euler's method, the solution to the system
 do i=N-1,1,-1
   Phi(i) = Phi(i+1)     -h *2d0*ygrid(i+1)*PhiDer(i+1)
   PhiDer(i)=PhiDer(i+1) -h *2d0*Phi(i+1)**1.5d0 
!      print* , xgrid(i), Phi(i)
 enddo
 mu=Phi(1)
! wrap things up. Here a file to save the function to disk
 open(unit=5,file='ThomasFermiPhi.dat')
 do i=1,N      
!Rescale the solution to satisfy the boundary condition 
   xgrid(i)=xgrid(i)*mu**0.33333333333d0
   Phi(i)=Phi(i)/mu
!Print to screen or save to disk
   print*, mu,xgrid(i),Phi(i)
!   write(5,*), xgrid(i),Phi(i)
 enddo
 close(5)
! Print the slope near the origin
! print*, (Phi(10)-Phi(1))/(xgrid(10)-xgrid(1))
!check that the solution does actually solve the Thomas-Fermi equation
!by comparing left and right hand side: the ratio should be 1
 do i=2,N-1,100
   print*, dsqrt(Phi(i)**3 / xgrid(i))/ ( (Phi(i+1)+Phi(i-1)-  &
   2d0*Phi(i))/((xgrid(i+1)-xgrid(i-1))/2d0)**2)
 enddo

end

















!Second part: compute and plot the form factor of the charge distribution
!First we return to the physical momentum grid by changing the variable
! pi=3.14159265d0
!choose the atom: Lead, Gold...
! Z=79
! To be completed by the student: computation of form factor as an integral 
! over the Phi(i) function
! do i=1,N
!   rgrid(i)= ???
!   rho_q(i)=FUNCION_DE(Phi)_A_OBTENER
! enddo
!now we set up the computation of the form factor in momentum space
!some physical parameters: coupling constant, electron mass in electronvolts,light
! alpha=1.d0/137d0
! me=511d3 
!length scale
! a0= 197d0/(alpha*me)
! a=a0*(9*pi**2/(128d0*Z))**0.3333333333d0
! print*, 'a0(nm) = ',a0, '  a_Z=  ',a
! do i=2,N
!   rgrid(i)= a * xgrid(i)
!charge density without the -eZ factor and with r**2 Jacobian 
!   rho_q(i)=(Phi(i)*2d0*me*197d0*alpha)**(3d0/2d0)*  &
!   dsqrt(Z*rgrid(i))*(8d0*pi/(3d0*(2d0*pi*197d0)**3))
!   print*, i,rgrid(i),rho_q(i)
! enddo
! open(unit=19,file='ChargeDensityThomasfermi.dat')
! do i=2,N
!   write(19,*) rgrid(i), rho_q(i)
! enddo
! close(19)
!
!momentum loop (in electron volts)
! qmax=1d3
! hq=qmax/dble(N+1)
! do jq=1,N 
! qgrid(jq)=dble(jq)*hq
!for each q, compute the form factor
! FF(jq)=0d0
! do i=2,N
!   FF(jq)=FF(jq)+rho_q(i)*dsin(qgrid(jq)*rgrid(i))  &
!   /(qgrid(jq)*rgrid(i))*((rgrid(i+1)-rgrid(i-1))/2.d0)
! enddo
! FF(jq)=FF(jq)*4d0*pi
! print*, qgrid(jq), FF(jq)
! enddo
!save the form factor to disk
! open(unit=7,file='Formfactor.dat')
! do jq=1,N
!   write(7,*) qgrid(jq),FF(jq)
! enddo
! close(7)


