Calculates disk spectrum using Frank, King, Raine 2002

disk_spek.f90

module Kind_Types
implicit none
save
  integer, parameter, public::    &
      single = selected_real_kind(p=6,r=37),   &
      double = selected_real_kind(p=13,r=200)
	!  p = precision
	!  r = range
	! return smallest kind of real value with a minimum
	!  of p decimal digits and maximum range >10^r
end module kind_types            
!*********************************************
module Natural_Constants
use Kind_Types
implicit none
save
real(kind=double), parameter, public::     &
PI = 3.141592654d0,    &
G = 6.67259d-8,        &   ! gravitational constant	(dyne.cm^2/s)
M_sun = 1.989d33,      &   ! mass of the Sun   ( g )
R_sun = 6.96d10,       &   ! Radius of the Sun
c = 2.99792458d10,     &   ! speed of light	  (cm/s)
k_B = 1.3806513d-16,     &   ! Boltzmann's constant (erg/K)
planck = 6.6252d-27,    &   ! Plancks constant
m_e = 9.11d-28,        &   ! electrons mass
m_p = 1.6725231d-24,   &   ! Mass of proton	  (g)
N_A = 6.02d23,         &   ! Avagadro's number
SB = 5.670399d-5,      &   ! Stefan-Boltzmann constant	(erg/cm^2.K^4.s)
thompson = 6.7d-25,    &   ! Thompson Cross Section
kappa = thompson/m_p,  &   ! opacity
day = 8.64d4,          &   ! 1 day is 86400 seconds
week = 7.d0*day,       &   !
month = 30.d0*day,     &   ! convert time units to seconds
year = 365.d0*day,     &   !
ly = c*year,           &
eV = 1.602177d-12,     &
pc = 3.26d0*ly,        &
kpc = pc*1.d3,         &
ms = 1.d-3,            &   ! millisecond in seconds
km = 1.d5,             &  ! km in cm's
v_min = 1.0d13,        &  ! minimum frequency
kappa_0 = 4.d25,      & ! *Z*(1+X) for Z metal and X hydrogen fraction http://arxiv.org/abs/astro-ph/0205212v1
alpha_SS = 0.1,       &  ! turbulent viscosity alpha
alpha = 19./16.,       &  ! alpha for electron scattering
lambda_max = c / v_min, &  ! minimum wavelength
kT_min = planck * v_min  !
!
end module Natural_Constants
!***************************************************
module Main_Parameters
use Natural_Constants
implicit none
save

real(kind=double), public::        &
                     !
Mass,     &	      ! mass of the star in "g"
!
Radius,            &
!
period,   & ! 0142
!
Jdisk = 2.80d50,       & ! assumed initial disk angular momentum
!
inclination = pi/3.d0,  &
!
L_X = 1.d35,   &
!
fract = 1.d0,   &
!
w_c = 0.9,      &
!
n = 1.d0,    & ! formerly jdot
!                      !
Mdot,  &
!
distance = 1.d0*kpc,   &
!
eta_d = 0.5d0,   & ! WCK
!
mmw = 2.0,    & ! mean molecular weight for metallic disk, used to be 0.65
!
Mdisk,	&
!
j_0, t0,nu_0,sigma_0,Mdot_0,C_0,r_0

integer, parameter, public::      &
order = 8,                        & ! then v_max = v_min*10^order
iMAX = 200

real(kind=double), public::	 v, Rin, Rout, t

end module Main_Parameters

!********************************
double precision function planck_dist(x)
!use Natural_Constants
use Main_Parameters
implicit none
real(kind=double), intent(in) :: x
real(kind=double):: hv_kT, T_0, T_1, T_diss, T_irr
real(kind=double):: planck_dist1, planck_dist2


T_0 = (3.d0*G*Mass*Mdot/(8.d0*PI*Rin**3*SB))**0.25

T_diss = T_0*((x**(-3))*(1.d0-n/x**0.5))**0.25

hv_kT = planck*v/(k_B*T_diss)

    if ( hv_kT <1.d-3) then

	    planck_dist1 =  x / hv_kT

    else if ( hv_kT > 12.d0 ) then

       planck_dist1 = x * dexp(-hv_kT)

	else

	     planck_dist1 = x / ( dexp(hv_kT) - 1.d0)

	end if




T_1 = (sqrt(k_B/(mmw*m_p))*L_X/(14.d0*pi*SB*sqrt(G*Mass*Rin**3)))**(2.d0/7.d0)

T_irr  = T_1*x**(-3.d0/7.d0)


	hv_kT = planck*v/(k_B*T_irr)

    if ( hv_kT <1.d-3) then

	    planck_dist2 =  x / hv_kT

    else if ( hv_kT > 12.d0 ) then

       planck_dist2 = x * dexp(-hv_kT)

	else

	     planck_dist2 = x / ( dexp(hv_kT) - 1.d0)

	end if


  !planck_dist = planck_dist1 + planck_dist2
  planck_dist = planck_dist1 ! shut off irradiation


end function planck_dist
!**************************************
program spektrum
use Main_Parameters
implicit none
double precision:: a, b, res
real(kind=double):: omega, R_co, R_L, R0_initial
real(kind=double):: const,  base, Fv, total
INTEGER ::  i, j, lines
character(len=8) :: fmt !format
character(len=12) :: filename
character(len=5) :: x1
external planck_dist

fmt = '(I5.5)'

!  reading no of lines
OPEN (unit=21, file="datain.dat",status="old")
READ (21, *) lines ! no of lines in input file
WRITE (*, *)  lines

do j=1, lines

READ (21, *) t, Mass, Radius, Mdot, Period, Mdisk

! setting up the filenames

write (x1,fmt) j
filename='wck'//x1//'.dat'

OPEN (unit=13, file=filename,status="replace")


omega = 2.d0*pi/period
!
R_co = (G*Mass / omega**2)**(1.d0/3.d0)

!initialize
!
j_0 = Jdisk/Mdisk
!
r_0 = (3.3558/1.7030)**2*j_0**2/(G*Mass)  !  replace
!
Sigma_0 = Mdisk / (4.d0*pi*r_0**2*3.3558d-5)
!
C_0 = alpha_SS**(8./7.)*(27.*kappa_0/SB)**(1./7.)    &
      *(k_B/(mmw*m_p))**(15./14.)*(G*Mass)**(-5./14.)
!
nu_0 = C_0*r_0**(15./14.)*Sigma_0**(3./7.)
!
t0 = r_0**2/(0.75d0*nu_0)
!
Mdot_0 = (alpha -1.0)*Mdisk/t0
!
!

if (j==1) then
  R0_initial = r_0 ! for the first step
endif 


Rin = Radius
!
Rout = r_0*(Mdot/Mdot_0)**(-2.*(1.-1./alpha))
!
!if (Rout < R0_initial .OR. Rout < Rin) then
!   Rout = r_0
!endif

!write(*,*) Mdot_0, Mdot, r_0/Rin, Rout/Rin
WRITE (*, *)  t, Rout/Rin
!
n = 1-(Rin/R_co)**(1.5)

const = 4.d0*pi*planck*cos(inclination)*(Rin/distance)**2/c**2


base = 10.d0**(dble(order)/dble(iMAX))


a = 1.d0
b = Rout/Rin

v = v_min

write(13,104) '# t = ',t, '/year       inclination = ', &
	inclination/pi , '*pi      distance = ', distance/kpc, '/kpc' 
write(13,'(A35,A53)') '#  lambda/cm      energy/eV   F_nue', &
		  '/(erg cm^-2 s^-1 Hz^-1)   lam*F_lam /(erg cm^-2 s^-1)'
	          
do i=1, imax

     
     call simpson(a,b,planck_dist,res)
	 

     Fv = res*const*v**3
 
     write (unit=13,fmt=103) c/v, planck*v/eV, Fv, v*Fv

     v = v_min*base**i

	 if (dlog10(v*Fv) < -40) then
	  exit
     end if

end do


100 FORMAT (1x, ES14.5, 2x, ES14.5, 2x, ES14.5)
103 FORMAT (1x, ES14.5, 2x, ES14.5, 2x, ES14.5, 2x, ES14.5)
104 FORMAT (1x,A6,ES12.4,1x,A26,F8.4,1x,A20,F6.2,1x,A4)

end do

end program spektrum

!********************************
!************************************
subroutine simpson(a,b,func,integral)
use Kind_Types
implicit none
integer, parameter:: n=400
integer:: i
real(kind=double), intent(in):: a,b
real(kind=double), intent(out)::integral
real(kind=double):: dx
real(kind=double), dimension(0:n):: c, x, y
real(kind=double), external :: func

! calculate the step, coefficients of formula

    dx = (b-a)/dble(n)
    c(0) = dx/3.d0
    do i = 1, n/2
        c(2*i-1) = (4.d0/3.d0)*dx
        c(2*i) = (2.d0/3.d0)*dx
    end do
    c(n) = c(0)
    integral = 0

! calculate the interpolate points and value on them

    do i = 0, n
        x(i) = a + i*dx
        y(i) = func(x(i))
    end do

! calculate the approximate integral

    do i = 0, n
        integral = integral + c(i)*y(i)
    end do

end subroutine simpson