program eliovar                        !inizio programma principale: nome
      parameter (nmax=1000)               !serve nelle dichiarazioni successive
      dimension r(nmax),chi(nmax),q(nmax),aux(nmax),vh(nmax)  !dimensioni array
      open(unit=70,status='unknown',file='f-vs-Zs.txt')  !70 = file f-vs-Zs.txt
      open(unit=80,status='unknown',file='E-vs-Zs.txt')  !80 = file E-vs-Zs.txt
      Z=2                                             !numero atomico dell'elio
      r1=0.01/Z                                !parametri della griglia radiale
      grigl=1.05
      al=alog(grigl)
      nx=alog(Z*75./0.005)/al
      r0=r1/sqrt(grigl)           !per il primo pezzetto dell'integrale radiale
      if(nx.gt.nmax) then                   !controllo dimensioni: ho sfondato?
         write(70,100) nx,nmax                                   !scrivo avviso
 100     format(1x,'stop: nx = ',i5,3x,'nmax = ',i5 )           !formato avviso
         stop                                              !stop se ho sfondato
      endif                                                      !fine ciclo if
      call griglia(r1,grigl,nx,nmax,r)         !costruisco griglia radiale r(i)
         do izz=1,21                             !ciclo su diversi valori di Z*
         iz=izz-16
         Zstar=Z+(float(iz)/16.)
            do i=1,nx                         !definisco chi(i), di simmetria s
            chi(i)=2.*(Zstar**1.5)*r(i)*exp(-Zstar*r(i))
            end do 
         coef=2.*(Zstar**1.5)      !coef = derivata prima radiale di chi in r=0
         an=(coef**2)*(r0**3)/3.     !pezzetto da r=0 al primo punto di griglia
            do i=1,nx                                    !ciclo: carica radiale
            an=an+(chi(i)*chi(i)*al*r(i))
            q(i)=an
            end do                                   !fine ciclo carica radiale
         an=(coef**2)*(r0**2)/2.     !pezzetto da r=0 al primo punto di griglia
            do i=1,nx              !ciclo: aux = integrale in dr di (chi*chi/r)
            an=an+(chi(i)*chi(i)*al)
            aux(i)=an 
            end do                                              !fine ciclo aux
         auxx=aux(nx)
            do i=1,nx          !ciclo: sottraggo auxx (costante d'integrazione)
            aux(i)=aux(i)-auxx !per avere campo elettrico -q(r)/r**2 a grandi r
            end do                                              !fine ciclo aux
            do i=1,nx                          !ciclo: potenziale di hartree vh
            vh(i)=(q(i)/r(i))-aux(i)
            end do                               !fine ciclo potenziale hartree
         hartree=0.0             !qui per il primo pezzetto va bene hartree = 0
            do i=1,nx                       !ciclo: energia di hartree e totale
            hartree=hartree+vh(i)*chi(i)*chi(i)*al*r(i)      !basate su fz. chi
            end do                                       !fine hartree e totale
      etotal=-(Zstar**2)+hartree+2.*(Zstar-Z)*auxx    !energia totale (variaz.)
      write(70,1000) 
      write(70,1001)(Zstar,i,r(i),chi(i)/r(i),q(i),vh(i),i=1,nx)!scrivo file 70
      write(80,1002) izz,Zstar,hartree, etotal                  !scrivo file 80
 1000 format(//)                                          !formato di scrittura
 1001 format(1x,f6.4,i5,5x,4e15.4)                        !formato di scrittura
 1002 format(1x,i5,5x,1p4e15.4)                           !formato di scrittura
      end do                               !fine ciclo sui diversi valori di Z*
      stop                                           !stop quando si arriva qui
      end                                                !fine modulo programma
C      
      subroutine griglia(r1,grigl,nx,nmax,r) !inizio subroutine: nome argomenti
      dimension r(nmax)                                     !dimensione array r
      r(1)=r1                                        !primo punto della griglia
         do i=2,nx                                  !ciclo: costruzione griglia
         r(i)=grigl*r(i-1)
         end do                                 !fine ciclo costruzione griglia
      return                                  !ritorna nel programma principale
      end                                               !fine modulo subroutine