program atmosfera

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

implicit none

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

integer, dimension (1:5), parameter  :: n=(/7,16,31,51,71/) 
                        ! Liczba poziomów dla kolejnych stopni jonizacji, ostatni jest
			! ekstrapolowany z pozosta³ych czterech za pomoc± wielomianu st. 2
integer :: r,s,i
real*8  :: T,Pe

Pe=131 ![dyn/(cm*cm)] taka jednostka cisnienia spowoduje, ¿e gêsto¶æ elektronów
       ! 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)

!--------- Koniec programu g³ównego - dalej same definicje funkcji!

contains

!--------- statistical weight
real*8 function g(r,s)
integer :: r,s
  g=1.0
end function

!--------- excitation energy
real*8 function chi(r,s)
integer :: r,s
  chi=s-1.0
end function

!--------- partition function (funkcja podzia³u)
real*8 function Ur(r,T)
integer :: r,s
real*8  :: T
  Ur=0.0
  do s=1,n(r)
    Ur=Ur+g(r,s)*exp(-chi(r,s)/(k*T))
  end do
end function

!--------- Boltzmann function
! Oblicza zawarto¶æ jonów na s-tym poziomie wzbudzenia w¶ród jonów
! o r-tym stopniu jonizacji w temperaturze T
real*8 function Boltz(r,s,T)
integer :: r,s
real*8  :: T
  Boltz=g(r,s)*exp(-chi(r,s)/(k*T))/Ur(r,T)
end function

!--------- Saha function
! Oblicza zawarto¶æ jonów o r-tym stopniu jonizacji w¶ród wszystkich jonów
! w temperaturze T, przy ci¶nieniu elektronowym Pe
real*8 function Saha(r,T,Pe)
integer :: r,rr
real*8  :: T,Pe,Ne,Sc
real*8, dimension(4) :: Saha_r
  Saha_r=0.
  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
  Saha=Saha_r(r)/sum(Saha_r)
end function

end program

