c Rozk³ad cz±steczek o okre¶lonym stopniu jonizacji i poziomie
c wzbudzenia zgodnie z prawami Sahy i Boltzmanna
c dla cz±steczki Schadeenium (E) - nieistniej±cego pierwiastka
c o nastêpuj±cych w³a¶ciwo¶ciach:
c 1) Energie jonizacji: chi_1=7 eV dla neutralnego E
c                       chi_2=16eV dla E^{+}
c                       chi_3=31eV dla E^{2+}
c                       chi_4=51eV dla E^{3+}
c 2) Energie wzbudzenia rosn± o 1eV na ka¿dy stopieñ: chi_{r,s}=s-1 eV
c 3) Wagi statystyczne wszystkich poziomów s± takie same: g_{r,s}=1

       implicit none

       real*8 k,kerg,h,me,Pi,T
       integer n(5)

       real*8 Pe,Ur,Boltz,Saha
       integer r,s,i       

       common k,kerg,h,me,Pi                       
       common /xx/ n
       data n/7,16,31,51,71/ 

c Sta³e fizyczne
       k=8.61734d-5      
c [eV/K]  - Sta³a Boltzmanna
       kerg=1.38065d-16  
c [erg/K] - Sta³a Boltzmanna
       h=6.62607d-27     
c [erg*s] - Sta³a Plancka
       me=9.10939d-28    
c [g]     - masa elektronu
       Pi=3.1415926      
c []      - liczba Pi



c Liczba poziomów dla kolejnych stopni jonizacji, ostatni jest
c ekstrapolowany z pozosta³ych czterech za pomoc± wielomianu st. 2

       Pe=131 
c [dyn/(cm*cm)] taka jednostka cisnienia spowoduje, ¿e gêsto¶æ elektronów
c wewn±trz funkcji Saha - bêdzie w elektronach/cm^3

       Pe=1d+3

       write(*,'(/a)') ' Part I - partition function - Ur'
       write(*,'(/4a9/)') 'r','5000K','10000K','20000K'
       do r=1,4
         write(*,'(i9,3f9.2)') r, Ur(r,5000d+0), Ur(r,10000d+0), 
     & Ur(r,20000d+0)
       end do

       write(*,'(/a)') ' Part II - Boltzmann law - n_(r,s)/N_r'
       write(*,'(/a9,3a15/)') 's','5000K','10000K','20000K'
       r=1
       do s=1,15
         write(*,'(i9,3d15.2)') s, Boltz(r,s,5000d+0), 
     &   Boltz(r,s,10000d+0), Boltz(r,s,20000d+0)
       end do


       write(*,'(/a,f7.1,a)') ' Part III - Saha law - N_r/N 
     &(for pressure Pe=',Pe,' [dyne/cm^2])'
       write(*,'(/a9,3a15/)') 'r','5000K','10000K','20000K'
       do r=1,5
         write(*,'(i9,3d15.2)') r , Saha(r,5000d+0,Pe), 
     & Saha(r,10000d+0,Pe), Saha(r,20000d+0,Pe)
       end do

       write(*,'(/a)') ' Part IV - Payne Curves for schadeenium '
       write(*,'(/a)') ' calculate ground-state populations n_(r,1)/N '
       write(*,'(a)') ' for pressure Pe=131 [dyne/cm^2]' 
       write(*,'(a)') ' file: payne_E.dat' 
       write(*,'(a)') ' one column for temperature i'
       write(*,'(a)') ' and four columns for populations 
     &in stages r=1,2,3,4 '
       write(*,'(a)') ' 300 raws for temperature T=100,...,30000K 
     &in steps of 100K'

      Pe=131
      open(1,file='payne_E.dat',status='unknown')
      do i=1,300
        T=100.0*i
      write(1,'(f10.0,4f20.15)') T, (Saha(r,T,Pe)*Boltz(r,1,T),r=1,4)  
      end do
      close(1)
      end

c--------- statistical weight
      real*8 function g(r,s)
      integer r,s
      g=1d0
      end 

c--------- excitation energy
      real*8 function chi(r,s)
      integer r,s
      chi=s-1d0
      end 

c--------- partition function (funkcja podzia³u)
      real*8 function Ur(r,T)
      integer r,s
      real*8 g,chi
       real*8 k,kerg,h,me,Pi,T
       integer n(1:5)
      common k,kerg,h,me,Pi
       common /xx/ n
      Ur=0.0
      do s=1,n(r)
        Ur=Ur+g(r,s)*exp(-chi(r,s)/(k*T))
      end do
      end 

c--------- Boltzmann function
c Oblicza zawarto¶æ jonów na s-tym poziomie wzbudzenia w¶ród jonów
c o r-tym stopniu jonizacji w temperaturze T
      real*8 function Boltz(r,s,T)
      integer r,s
      real*8 g,chi,Ur
       real*8 k,kerg,h,me,Pi,T
       integer n(1:5)
      common k,kerg,h,me,Pi
       common /xx/ n
      Boltz=g(r,s)*exp(-chi(r,s)/(k*T))/Ur(r,T)
      end 

c--------- Saha function
c Oblicza zawarto¶æ jonów o r-tym stopniu jonizacji w¶ród wszystkich jonów
c w temperaturze T, przy ci¶nieniu elektronowym Pe
      real*8 function Saha(r,T,Pe)
      integer r,rr
      real*8 Pe,Ne,Sc,Ur
      real*8 Saha_r(4)
       real*8 k,kerg,h,me,Pi,T,sum
       integer n(1:5)
      common k,kerg,h,me,Pi
       common /xx/ n
      Saha_r(1)=0d0
      Saha_r(2)=0d0
      Saha_r(3)=0d0
      Saha_r(4)=0d0
      Ne=Pe/(kerg*T)
      Sc=(2/Ne)*(((2*Pi*me*kerg*T)/(h*h))**1.5)
      Saha_r(1)=1.
      do rr=1,3
         Saha_r(rr+1)=Saha_r(rr)*Ur(rr+1,T)*Sc*exp(-n(rr)/
     &   (k*T))/(Ur(rr,T))
      end do
      sum=0d0
      do i=1,4
       sum=sum+Saha_r(i)
      end do
      Saha=Saha_r(r)/sum
      end 


