Sph_infall is a program to model the intensity of an optically-thin, spherically-symmetric envelope with infall motion. The envelope, with an inner radius R_inn and an outer radius R_out, has an infall velocity and a temperature that are power laws of the radius with indices −α and −β. As a consequence of the continuity equation, the density is also a power law with index α−2.
V = V_0 (R/R_0)^(−α),
T = T_0 (R/R_0)^(−β),
ρ = ρ_0 (R/R_0)^(α−2).
The usual case is α = β = 1/2. R_0 is an arbitrary reference radius. In general, a line-of sight at a projected radius p will intersect the equal-line-of-sight-velocity surface V_z at two points. The contribution to the observed intensity I (p, V_z) of any of these points, at a distance R_i from the center, is evaluated as
I(p, V_z) = T(R_i) ρ(R_i) |dz/dV_z|(R_i)
That is, the intensity coming from each point at radial velocity V_z is proportional to the temperature, T(R), times the column density, n(R)∆z, and we approximate ∆z =|dz/dV_z |∆V_z , where ∆V_z is constant and can be omitted.
The model of the infalling spherical envelope is described by the following parameters.
Npix, ∆x: Map size and pixel size (arcsec).
Nchan, ∆Vchan : Number of velocity channels and channel width (km s^(−1) ).
HPBW: Beamwidth (arcsec). The beam is assumed to be circular because the 2D-convolution is calculated along one direction only.
∆V : Linewidth (without infall broadening) (km s^(−1) ).
R_inn , R_out: Envelope inner and outer radii (arcsec).
α, β: Velocity and temperature power-law indices, V ∝ R^ (−α), T ∝ R^( −β ).
V_0: Infall velocity (km s^(−1) ) at the reference radius (1" ).
Distances are in units of R_0, r = R/R_0, velocities in units of V_0, v = V/V_0, and temperatures in units of T_0, t = T/T_0. A radial cut of the envelope is calculated first, up to the outer radius r_out . We know that for any projected radius p, the minimum and maximum temperatures along the line of sight are the temperatures at the outer and inner radius of the envelope, t_min = r_out^(−β) and t_max = r_inn^(−β).
For each value of the projected radius p, all the velocity channels are initially set to zero. The maximum absolute value of the radial velocity is v_m = α^(α/2) (1 + α)^[−(1+α)/2] p^(−α) (or v_m = (4/27)^(1/4) p^(−1/2) for α = β = 1/2). For each channel with a velocity |v_z | below v_m, the two temperatures t_low and t_high are calculated. If any of the two temperatures is t_min ≤ t ≤ t_max , the point is inside the envelope and contributes to the intensity. Once the intensity at this velocity is calculated, a Gaussian line with the calculated intensity, with linewidth ∆V , centered on the velocity of the channel, is added to the velocity channels corresponding to the projected radius p.
Once the radial cut is calculated, the 2D channel maps (npix × npix) are estimated by interpolation onto a regular grid in the plane of the sky. Finally, the channels maps are 2D-convolved with a Gaussian beam. The 2D-convolution is calculated along one direction only, because the model and beam are circularly symmetric. The convolution is made through direct calculation and not via Fourier transform because the intensity has abrupt changes at the edges of the envelope and the Fourier transform has strong lobes. Finally, the radial cut output from the convolution process is again interpolated onto a regular grid in the plane of the sky to obtain the final model data cube.
Figure 1a: Typical spectrum for an infalling shell of radius R = 20", ∆V = 0.1 km s^(−1), V_0 = 0.9 km s^(−1) , α = β = 0.5. The spectra show emission from the blueshifted back half-shell and redshifted front half-shell.
Figure 1b: Typical spectrum for an infalling envelope with R_inn = 0.1", R_out = 20", ∆V = 0.1 km s^(−1), V_0 = 0.9 km s^(−1) , α = β = 0.5. The spectra show emission from the blueshifted back half-envelope and redshifted front half-envelope.
Fig. 2: Position-velocity plots for an infalling shell of radius R = 20" (left) and an envelope with R_inn = 0.1" , R_out = 20" (right). The rest of parameters are ∆V = 0.1 km s^(−1) , V_0 = 0.9 km s^(−1) , α = β = 0.5. The beam size is 2" .
Linux or Mac OS X system.
Fortran compiler, e.g. gfortran.
PGplot Graphics Subroutine Library (www.astro.caltech.edu/~tjp/pgplot/) compiled with gfortran. For the installation of PGplot with gfortran, follow the installation procedure described in pgplot_gfortran. Mac users can find a guide for installing PGplot on Mac OS X in mingus.as.arizona.edu/~bjw/software/pgplot_fix.html. You can find a local copy of the same page in pgplot_macosx.
Once downloaded, untar the tar file and compile the program with the PGplot and X11 libraries
gfortran -o sph_infall -l X11 -l pgplot sph_infall.f90
You may need to indicate the location of the PGplot or X11 libraries, adding the options, for instance,
-L /usr/local/pgplot -L /usr/X11/lib