function bplanck(temp,nu) implicit none double precision :: temp double precision :: nu double precision :: bplanck if(temp.eq.0.d0) then bplanck = 0.d0 return endif bplanck = 1.47455d-47 * nu * nu * nu / & (exp(4.7989d-11 * nu / temp)-1.d0) + 1.d-290 return end function bplanck program problem_setup implicit none ! ! Declarations of functions ! double precision :: bplanck ! ! Constants ! doubleprecision :: GG,mp,kk,cc,ss parameter(GG = 6.673d-8) ! Gravitational constant parameter(mp = 1.6726d-24) ! Mass of proton [g] parameter(kk = 1.3807d-16) ! Bolzmann's constant [erg/K] parameter(cc = 2.997924580000000d10) ! Light speed [cm/s] parameter(ss = 5.6703d-5) ! Stefan-Boltzmann const [erg/cm^2/K^4/s] doubleprecision :: muh2 parameter(muh2= 2.3000d0) ! Mean molec weight H2+He doubleprecision :: LS,RS,MS,AU,pc parameter(LS = 3.8525d33) ! Solar luminosity [erg/s] parameter(RS = 6.96d10) ! Solar radius [cm] parameter(MS = 1.98892d33) ! Solar mass [g] parameter(AU = 1.49598d13) ! Astronomical Unit [cm] parameter(pc = 3.08572d18) ! Parsec [cm] doubleprecision pi parameter(pi = 3.1415926535897932385d0) ! ! Parameters of the model ! double precision, parameter :: rstar=5.3*RS ! Radius of the star double precision, parameter :: mstar=MS ! Mass of the star double precision, parameter :: tstar=5000. ! Temperature of the star integer, parameter :: nr=200 ! Nr of radial grid points double precision, parameter :: rin=200*AU ! Inner radius of the dust shell double precision, parameter :: rout=1.d4*AU ! Outer radius of the dust shell double precision, parameter :: r0=1d3*AU ! Reference radius double precision, parameter :: plrho=-1.8 ! Powerlaw index for density double precision, parameter :: rhodust0=1d-18 ! Dust density at r0 integer, parameter :: nl=100 ! Nr of wavelengths double precision, parameter :: lambda0=1d-1 ! Shortest wavelength in micron double precision, parameter :: lambda1=1d+4 ! Longest wavelength in micron integer, parameter :: nphot=100000 ! Nr of photon packages ! ! The arrays ! double precision :: rc(nr) ! Cell center radii double precision :: ri(nr+1) ! Cell wall radii double precision :: rhodust(nr) ! Dust density array double precision :: lambda(nl) ! Wavelength grid double precision :: starspectrum(nl) ! F_nu of star at d=1*parsec ! ! Some other variables ! double precision :: nu integer :: ir,il ! ! ================== MAKE THE MODEL ================ ! ! Make the radial grid ! do ir=1,nr+1 ri(ir) = rin * (rout/rin)**((ir-1.d0)/(1.d0*nr)) enddo do ir=1,nr rc(ir) = sqrt(ri(ir)*ri(ir+1)) enddo ! ! Make the density ! do ir=1,nr rhodust(ir) = rhodust0 * (rc(ir)/r0)**plrho enddo ! ! Make the wavelength grid ! do il=1,nl lambda(il) = lambda0 * (lambda1/lambda0)**((il-1.d0)/(nl-1.d0)) enddo ! ! Make the stellar spectrum (a simple blackbody spectrum, flux seen ! at a distance of 1 parsec) ! do il=1,nl nu = 1d4*cc/lambda(il) starspectrum(il) = pi * rstar**2 * bplanck(tstar,nu) / pc**2 enddo ! ! ============ WRITE THE MODEL INPUT FILES FOR RADMC-3D ========== ! ! Write the grid file ! open(unit=1,file='amr_grid.inp') write(1,*) 1 ! Format identifier write(1,*) 0 ! 0 = regular grid write(1,*) 100 ! 100 = spherical coordinates write(1,*) 0 ! (always 0) write(1,*) 1,0,0 ! Include only r-coordinate, not theta nor phi write(1,*) nr,1,1 ! Grid size do ir=1,nr+1 write(1,*) ri(ir) ! The cell interface radii enddo write(1,*) 0.d0 ! Theta_min write(1,*) pi ! Theta_max write(1,*) 0.d0 ! Phi_min write(1,*) 2*pi ! Phi_max close(1) ! ! Write the density file ! open(unit=1,file='dust_density.inp') write(1,*) 1 ! Format identifier write(1,*) nr ! Number of cells write(1,*) 1 ! Nr of dust species do ir=1,nr write(1,*) rhodust(ir) enddo close(1) ! ! Write the wavelength file ! open(unit=1,file='wavelength_micron.inp') write(1,*) nl do il=1,nl write(1,*) lambda(il) enddo close(1) ! ! Write the stars.inp file ! open(unit=1,file='stars.inp') write(1,*) 2 ! Format identifier write(1,*) 1,nl ! One star, nl wavelengths write(1,*) ' ' write(1,*) rstar,mstar,0.d0,0.d0,0.d0 ! Radius, mass and position of the star write(1,*) ' ' do il=1,nl write(1,*) lambda(il) enddo write(1,*) ' ' do il=1,nl write(1,*) starspectrum(il) enddo close(1) ! ! Dust opacity control file ! open(unit=1,file='dustopac.inp') write(1,*) '2 Format number of this file' write(1,*) '1 Nr of dust species' write(1,*) '============================================================================' write(1,*) '1 Way in which this dust species is read' write(1,*) '0 0=Thermal grain' write(1,*) 'silicate Extension of name of dustkappa_***.inp file' write(1,*) '----------------------------------------------------------------------------' close(1) ! ! The RADMC-3D control file ! open(unit=1,file='radmc3d.inp') write(1,*) 'nphot = ',nphot close(1) ! end program problem_setup