program atmosfera

! Rozk³ad cz±steczek o okre¶lonym stopniu jonizacji i poziomie
! wzbudzenia zgodnie z prawami Sahy i Boltzmanna
! dla cz±steczki wodoru
! 1) Energia jonizacji: chi_1=13.598 eV 
! 2) Energie wzbudzenia dla neutralnego wodoru: chi_{1,s}=13.598*(1-1/(s*s)) eV
! 3) Wagi statystyczne dla neutralnego wodoru: g_{1,s}=2*s*s
! 4) Waga statystyczna i funkcja podzia³u dla jednokrotnie zjonizowanego wodoru (proton):
!     U_2=g_{2,1}=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
real*8, parameter         :: chi1=13.598       ![eV]    - energia jonizacji atomu H

integer, parameter        :: n=100
                        ! Liczba poziomów wzbudzenia (s_max)  
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+2

!---------
write(*,'(/a)') ' Part I - partition function - Ur'
write(*,'(/4a9/)') 'r','5000K','10000K','20000K'
do r=1,2
  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,100,10
  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,2
  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_H.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_H.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,2)  
end do
close(1)

!--------- Koniec programu g³ównego!

contains

!--------- statistical weight
real*8 function g(r,s)
integer :: r,s
  if (r.eq.1) then
    g=2.0*s*s
  else
    if ((r.eq.2).and.(s.eq.1)) then
      g=1.0
    else
      g=0.0
    endif
  endif
end function

!--------- excitation energy
real*8 function chi(r,s)
integer :: r,s
  if (r.eq.1) then 
    chi=13.598*(1.0-1.0/(s*s))
  else
    chi=0.0
  endif
end function

!--------- partition function (funkcja podzia³u)
real*8 function Ur(r,T)
integer :: r,s
real*8  :: T
  if (r.eq.1) then 
    Ur=0.0
    do s=1,n
      Ur=Ur+g(r,s)*exp(-chi(r,s)/(k*T))
    end do
  else
    Ur=1.0
  endif
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.
  Saha_r(2)=Saha_r(1)*Ur(2,T)*Sc*exp(-chi1/(k*T))/(Ur(1,T))
  Saha=Saha_r(r)/sum(Saha_r)
end function

end program

