module my_type
INTEGER, PARAMETER :: IB = SELECTED_INT_KIND(9) !integer*4
INTEGER, PARAMETER :: SP = SELECTED_REAL_KIND(6,37) !real*4
INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15,307) !real*8
end module my_type
program limb
! shall calculate the angular dependence of photon's emitted from
! a plane-parallel, grey atmosphere of radial optical depth taumax
! Note that tau and taumax are RADIAL optical depths.
! Since photons have a direction of MU with respect to the radial
! direction, we have to account for projection effects:
! in order to calculate this quantity, i.e., to decide whether a
! photon is
! - inside the atmoshere (0 < tau < taumax), &
! - has left the atmosphere (tau < 0)
! - or has been scattered back into the stellar core (tau > taumax)
! the actual value of tau(radial) has to be updated according to
! tau(radial) = tau(radial) - tau_i * mu
! if tau_i is the actual path-length for the photon as drawn from the
! corresponding pdf and mu its direction, also drawn from the corresponding
! pdf.
! If a photon was found to have been scattered back into the core,
! just release a new photon at the inner boundary, however without
! updating the number of photons since we want to have NPHOT photons LEAVING
! the OUTER atmosphere!!!
use my_type
implicit none
! this is the array into which we will sort the emergent photons
integer(ib),allocatable,dimension(:) :: muarray
integer(ib) :: na,nphot,l,i,isum,norm
real(sp) :: dmu,taumax,mu,dist
print*,'give in number of channels' ! to check for resolution effects
read*,na
! allocate the angular array and initialize it accounting for
! the fact that no photons have been accumulated so far
dmu=??? !widths of channels
!provide input for photon number and total optical depth
do l=1,nphot !big loop for all photons
! follow the photon's path from emission to final escape under
! angle MU
! Finally (tau le 0), the photon should have left the photosphere, and
! we can accumulate the angular distribution
i=??? ! How do we find the index of the channel according to angle mu?
! Remember, we have NA channels, with widths DMU
! the first channel should be used for all photons with MU = 0...DMU,
! the last one for photons with MU= 1-DMU...1
! Note: sometimes, mu is exactly unity (or very close to it, so that
! numerically it IS unity), which would
! result in i=na+1. Thus: if i > na, set i=na.
! careful testing required.
muarray(i) = muarray(i) + 1 !one more photon into this channel
enddo
!test
isum=0
do i=1,na
isum=isum+muarray(i)
enddo
! Just a test for securacy. Have we missed any photon?
if(isum.ne.nphot) stop' not all photons counted'
!preparation of output
open(1,file='limb.out')
norm=muarray(na) !normalizing to the last channel
do i=1,na
mu=(i-0.5)*dmu !to be centered!
dist=muarray(i)/float(norm)
write(*, 99) mu,dist
write(1,100) mu,dist
enddo
close(1)
99 format(' angle ',f6.4,' photon-dist ',f6.4)
100 format(f6.4,2x,f6.4)
end