** Nebular Luminosity Spectrum **
November 22, 2008: E.Komatsu
Ref: Fernandez & Komatsu, ApJ, 646, 703 (2006)

In "nebular_spectrum.f90" we provide double-precision functions that 
return the luminosity spectra of various radiation processes in the 
nebular around a star for a given nu and QH, where:

- nu: frequency, in units of Hz
- QH: the number of hydrogen-ionizing photons per second, in units of s^-1

The available functions are:

- Lnu_ff(nu,QH): Free-free emission [erg s^-1 Hz^-1] 
- Lnu_fb(nu,QH): Free-bound emission [erg s^-1 Hz^-1] 
- Lnu_2g(nu,QH): Two-photon emission [erg s^-1 Hz^-1] 
- L_lya(nu,QH): Lyman-alpha line emission [erg s^-1] 

Note that the Lyman-alpha routine returns the luminosity, while
the other routines return the luminosity per frequency.

- To compile and use the sample program (given below), edit Makefile
and simply "./make"
- It will generate an executable called "compute_nebular_spectrum"

=========================================================================
A simple program, "compute_nebular_spectrum," for generating nu*Lnu as 
a function of nu is

PROGRAM Compute_Nebular_Spectrum
  ! A sample program for computing nu*L_nu (in units of erg s^-1), where nu
  ! is in units of Hz, of nebular emission around a star, including:
  ! - Free-free
  ! - Free-bound
  ! - Two-photon
  ! - Lyman-alpha line
  ! The output file contains:
  ! 1st column: h*nu [eV]
  ! 2nd column: nu*L_nu for free-free  [erg s^-1]
  ! 3rd column: nu*L_nu for free-bound [erg s^-1]
  ! 4th column: nu*L_nu for two-photon [erg s^-1]
  ! 5th column: L for Lyman-alpha line [erg s^-1]
  ! REF: Fernandez & Komatsu, ApJ, 646, 703 (2006)
  ! November 22, 2008: E.Komatsu
  IMPLICIT none
  double precision :: Lnu_ff,Lnu_fb,Lnu_2g,L_lya
  double precision :: hnu,nu,QH
  double precision :: h=4.1356668d-15! Planck's constant [eV s]
  QH=1d49 ! the number of hydrogen-ionizing photons per second
  open(1,file='hnu_ff_fb_2g_lya.txt',status='unknown')
  do hnu=1,15,0.01
     nu=hnu/h
     write(1,'(5E13.5)'),hnu,nu*Lnu_ff(nu,QH),nu*Lnu_fb(nu,QH),&
          nu*Lnu_2g(nu,QH),L_lya(nu,QH)
  enddo
  close(1)
END PROGRAM Compute_Nebular_Spectrum