Embedded Files
Envelope_plummer. Version 2014
Envelope_plummer. Version 2014
Envelope_plummer is a program to fit simultaneously the intensity radial profile and the spectral energy distribution (SED) of an envelope. The envelope is assumed to be spherically symmetric, with a Plummer-like radial profile of density
rho = rho_c [1+(r/r_c)^2]^(-qd/2)
and a power-law radial profile of temperature,
T = T_0 (r/r_0)^(-qt).
The degeneracy between qd and qt, and between the density and temperature scales rho_c and T_0 in the observed intensity profiles is broken by fitting the SED, which gives unambiguously the temperature power-law index qt and the temperature scale T_0. No assumptions of optically thin emission are made. Additional free-free emission or dust emission from an unresolved central disk are also included in the model. For a description of the model see Palau et al. (2014).
The program code (Fortran) can be downloaded here. The publication you have to cite if you use this program is:
Palau, A., Estalella, R., Girart, J. M. et al. (2014), ApJ, 785, 42
If you want to use the program and need support, you can contact me at robert.estalella(at)ub.edu
Model parameters
Model parameters
beta, dust opacity power-law index, which gives the temperature power-law index, qt= 2/(4+beta).
T_0, temperature at r_0 = 1 kau = 1000 au
rho_0, nominal mass density at r_0 = 1 kau (extrapolated from the power law for r >> r_c, true density if r_c << r_0)
rc/r0, critical radius, in kau. For r << r_c, rho(r) = rho_c (constant)
qd, density power-law index for r >> r_c
Derived parameter
Derived parameter
rho_c = rho_0 (r_c/r_0)^(-qd), central density at r = 0
Outputs
Outputs
Radial intensity profile convolved with a Gaussian beam (main and error beams) including optional beam chopping
Spectral energy distribution
Fig, 1: Example of graphic output from envelope_plummer. Left: Intensity radial profiles at 850 and 350 um; observations in blue with error bars, model in black solid line, beam in red dotted line. Right: SED (top) and fit logarithmic residual (bottom); observations (top) or residual (bottom) in blue with error bars, model values as red dots, model for a constant aperture in black solid line.
Requisites
Requisites
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.
Installation
Installation
Once downloaded, untar the tar file and compile the program with the PGplot and X11 libraries
gfortran -o envelope_plummer -l X11 -l pgplot envelope_plummer.f
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
Versions history
Versions history
2001 Jul: 1st version.
2002 Oct: Adapted for L723 at 850 and 450 um, using PGplot.
2003 Feb: Added SED fit at other frequencies.
2003 Jul: Added an unresolved central disk..
2007 May: Units changed, rewriting most code.
2013 Feb: Adapted to any source. All internal units cgs. Fitting of model parameters, program documentation.
2013 Jun: New improved fit engine, with Halton sampling
2014 Oct: Beam chop optional. Chop distance in par file. Dist =< 0: no chop.
2014 Oct: Plummer-like density profile.
2014 Nov: Error taking into account actual number of parameters fitted. Added calculation of Herschel fluxes at 160, 250, 350, 500 um and average column density.
2014 Dec: Added plot update for every best fit found.
Example of input parameter file: test_plummer.par
Example of input parameter file: test_plummer.par
"test_source" !Source name
1980.0 !Distance (pc)
345.0 !870 um: Freq for radial profile (GHz)
1.00 !870 um: Main beam relative height (rhmb+rheb= 1)
22.0 !870 um: HP main BW for radial profile (arcsec)
0.00 !870 um: Error beam relative height (rhmb+rheb= 1)
0.0 !870 um: HP error BW for radial profile (arcsec)
0.10 !870 um: relative error in flux scale
0.00 !870 um chopthrow
"test_source-ellint-15oct-east.dat" !850 um: file with radial profile
856.55 !350 um: Freq for radial profile (GHz)
1.00 !350 um: Main beam relative height (rhmb+rheb= 1)
9.60 !350 um: HP main BW for radial profile (arcsec)
0.00 !350 um: Error beam relative height (rhmb+rheb= 1)
0.00 !350 um: HP error BW for radial profile (arcsec)
0.50 !350 um: relative error in flux scale
0.00 !350 um chopthrow
"test_source-ellint-15oct-east.dat" !350 um: file with radial profile
0.20 !Relative error in flux scale for SED data
"test_source_sed_Fnu_lam_scaled.dat" !file with SED
230.0 !Reference frequency for dust opacity (GHz)
0.008991 !Dust opacity at reference frequency (cm^2 g^-1)
1.5 !Dust opacity power-law index: freq^beta
300.0 !Temperature at r_0= 1000 AU (K)
2e-16 !Central density (for r << r_c) (g cm^-3)
2.0 !Radius r_c/r_0 (kAU)
1.8 !Density power-law index (for r >> r_c) r^{-qd}
12.0 !Envelope minimum temperature (K)
5.0e3 !Envelope minimum number density (cm^-3)
0.0 !Disk flux density at 7 mm (mJy)
43.31 !Reference frequency for disk flux (GHz)
1.0 !Disk dust opacity power-law index, spect_indx= 2+beta_disk
100 !Disk radius (AU)
8.3 !Reference frequency for free-free emission (GHz)
0.00 !Free-free emission at reference frequency (mJY)
1.9 !Free-free spectral index
Example of intensity radial profile
Example of intensity radial profile
#***** Integrate a Miriad image in elliptical annuli *****
# bunit: Jy/beam frequency: 0.000000 K/Jy: 1.0000
# Axes: RA---GLS x DEC--GLS Beam: 22.00 x 22.00 arcsecs
# Center: -69.6 -333.9 Ellipse pa: 0. Incline: 0. Not P.B. Corrected
# Radius(") Pixels Average rms Ann. Sum Cum. Sum
5.000 4.000 6.232 0.101 24.928 24.928
10.000 9.000 5.533 0.251 49.793 74.721
15.000 13.000 4.438 0.308 57.690 132.411
20.000 21.000 3.206 0.494 67.331 199.742
25.000 22.000 2.452 0.597 53.950 253.692
30.000 27.000 1.890 0.496 51.024 304.716
35.000 33.000 1.462 0.476 48.260 352.976
40.000 35.000 1.139 0.503 39.857 392.833
45.000 25.000 0.720 0.199 18.010 410.843
50.000 25.000 0.555 0.104 13.886 424.730
Example of SED
Example of SED
#lambda Flux Flux_err Aperture
#(um) (Jy) (Jy) (arcsec)
70. 50.079 2.702 12.87
160. 177.209 6.974 15.85
350. 138.412 2.873 27.21
351. 87.937 0. 20.0
500 53.909 2.337 27.21
870. 24.06 3.85 64.14
Example of log from envelope_plummer
Example of log from envelope_plummer
_______________________________________________________________________________
Envelope model, Plummer-like density profile.
Version 2014/11. Robert Estalella
Date: 2014/12/03 Time: 17:55:43
Reading parameters from file test_plummer.par
name= test_source
dist_pc= 1980.000
freq_op_GHz= 345.000
rhmb_op= 1.000
hpmbw_op_arcsec= 22.000
rheb_op= 0.000
hpebw_op_arcsec= 0.000
scale_error_op= 0.100
chopthrow_arcsec= 0.000
Reading radial profile data from file test_source-ellint-15oct-east.dat
freq_op_GHz= 856.550
rhmb_op= 1.000
hpmbw_op_arcsec= 9.600
rheb_op= 0.000
hpebw_op_arcsec= 0.000
scale_error_op= 0.500
chopthrow_arcsec= 0.000
Reading radial profile data from file test_source-ellint-15oct-east.dat
scale_error_os= 0.200
Reading SED data from file test_source_sed_Fnu_lam_scaled.dat
freq_opac_GHz= 230.000
dustopac= 0.009
beta= 1.500
temp0= 300.000
dens0= 2.000E-16
rc0= 2.000
qd= 1.800
T_min= 12.000
dn_min= 5000.000
flux_disk_mJy= 0.000
freq_disk_GHz= 43.310
beta_disk= 1.000
r_disk_au= 100.000
freqf_GHz= 8.300
fluxf_mJy= 0.000
spindf= 1.900
Reading fit parameters from file model_plummer_fitpar.dat
Source name= test_source
Distance (pc)= 1.980E+03
Resolution (arcsec)= 6.000E-01
(au)= 1.188E+03
Mod max radius (arcsec)= 1.200E+02
(au)= 2.376E+05
Dust opacity index beta= 9.140E-01
rho_0 (r=1kau)(g cm^-3)= 4.062E-16
rho_c (r=r_c) (g cm^-3)= 1.662E-18
r_c (kau)= 1.540E+01
rho power-law index qd= 2.011E+00
n_min (cm^-3)= 5.000E+03
rho_min (g cm^-3)= 2.343E-20
r_env= r(rho_min) (au)= 1.273E+05
(pc)= 6.171E-01
(arcsec)= 6.428E+01
r_env/r_c= 8.265E+00
Mass (analytic) (M_sol)= 9.616E+02
Total mass (num)(M_sol)= 8.598E+02
Radius Avg dens Mass
(cm^-3) (M_sol)
0.05 pc 3.084E+05 1.118E+01
0.10 pc 1.730E+05 5.016E+01
0.15 pc 1.171E+05 1.146E+02
0.20 pc 8.153E+04 1.891E+02
r_c 2.291E+05 2.765E+01
r_env 1.262E+04 8.598E+02
T_0= T(1kau) (K)= 8.611E+01
T power-law index qt= 4.070E-01
T_min (K)= 1.200E+01
r(T_min) (au)= 1.267E+05
Radial profile wl (um)= 850
Peak intens. (Jy/beam)= 6.115E+00
Flux density (Jy)= 2.494E+01
Profile fit to 24 points
Fit bias, chi^2, chi_r= -0.026 5.222 0.524
Radial profile wl (um)= 350
Peak intens. (Jy/beam)= 1.761E+01
Flux density (Jy)= 1.366E+02
Profile fit to 23 points
Fit bias, chi^2, chi_r= 0.493 16.175 0.948
SED fit to 6 points
Fit bias, chi^2, chi_r= -0.036 15.077 3.883
Radial profiles and SED: np chi^2 chi_r
Profile at 850 um: 24 5.222 0.524
Profile at 350 um: 23 16.175 0.948
SED: 6 15.077 3.883
Total: 53 36.474 0.872
Model fit
Par: beta T_0 rho_0 rc/r0 qd chi_r
Val: 0.914 86.108 4.062E-16 15.401 2.011 0.872
Rnge 1.500 300.000 1.000E-16 10.000 1.000
Iteration parameters: Niter Nloop LoopFct
200 10 0.700
Model fit
Par: beta T_0 rho_0 rc/r0 qd chi_r
Val: 0.914 86.108 4.062E-16 15.401 2.011 0.872
Rnge 1.500 300.000 1.000E-16 10.000 1.000
Iteration parameters: Niter Nloop LoopFct
20 5 0.700
Model fit
Par: beta T_0 rho_0 rc/r0 qd chi_r
Val: 1.000 100.000 1.000E-16 20.000 2.000 2.941
Rnge 1.500 300.000 1.000E-16 10.000 1.000
Iteration parameters: Niter Nloop LoopFct
20 5 0.700
Model fit
Par: beta T_0 rho_0 rc/r0 qd chi_r
Val: 1.000 100.000 1.000E-16 20.000 2.000 2.941
Rnge 0.500 20.000 1.000E-17 1.000 0.500
Iteration parameters: Niter Nloop LoopFct
20 5 0.700
Par: beta T_0 rho_0 rc/r0 qd chi_r
Loop 1 0.989 119.909 1.082E-16 20.081 1.871 1.675
Rnge 0.500 20.000 1.000E-17 1.000 0.500
Loop 2 0.793 96.424 1.076E-16 19.638 1.706 1.395
Rnge 0.350 14.000 7.000E-18 0.700 0.350
Loop 3 0.691 107.035 1.151E-16 20.767 1.723 1.362
Rnge 0.245 9.800 4.900E-18 0.490 0.245
Loop 4 0.931 97.298 1.117E-16 20.946 1.705 1.320
Rnge 0.171 6.860 3.430E-18 0.343 0.171
Loop 5 0.796 108.920 1.129E-16 20.347 1.747 1.262
Rnge 0.120 4.802 2.401E-18 0.240 0.120
beta T_0 rho_0 rc/r0 qd chi_r
Val: 0.887 102.828 1.118E-16 20.192 1.743 1.253
Err: 0.061 4.659 1.327E-17 5.837 0.023
Err+ 0.058 4.659 1.973E-17 1.483 0.019
Err- 0.064 4.659 6.811E-18 10.190 0.027
Source name= test_source
Distance (pc)= 1.980E+03
Resolution (arcsec)= 6.000E-01
(au)= 1.188E+03
Mod max radius (arcsec)= 1.200E+02
(au)= 2.376E+05
Dust opacity index beta= 8.870E-01
rho_0 (r=1kau)(g cm^-3)= 1.118E-16
rho_c (r=r_c) (g cm^-3)= 5.928E-19
r_c (kau)= 2.019E+01
rho power-law index qd= 1.743E+00
n_min (cm^-3)= 5.000E+03
rho_min (g cm^-3)= 2.343E-20
r_env= r(rho_min) (au)= 1.273E+05
(pc)= 6.169E-01
(arcsec)= 6.427E+01
r_env/r_c= 6.302E+00
Mass (analytic) (M_sol)= 7.827E+02
Total mass (num)(M_sol)= 7.006E+02
Radius Avg dens Mass
(cm^-3) (M_sol)
0.05 pc 1.231E+05 4.462E+00
0.10 pc 8.071E+04 2.340E+01
0.15 pc 6.253E+04 6.120E+01
0.20 pc 4.795E+04 1.112E+02
r_c 8.604E+04 2.340E+01
r_env 1.029E+04 7.006E+02
T_0= T(1kau) (K)= 1.028E+02
T power-law index qt= 4.092E-01
T_min (K)= 1.200E+01
r(T_min) (au)= 1.904E+05
Radial profile wl (um)= 850
Peak intens. (Jy/beam)= 4.493E+00
Flux density (Jy)= 2.493E+01
Profile fit to 24 points
Fit bias, chi^2, chi_r= -0.052 32.308 1.304
Radial profile wl (um)= 350
Peak intens. (Jy/beam)= 1.135E+01
Flux density (Jy)= 1.456E+02
Profile fit to 23 points
Fit bias, chi^2, chi_r= 0.218 22.259 1.112
SED fit to 6 points
Fit bias, chi^2, chi_r= 0.511 20.819 4.563
Radial profiles and SED: np chi^2 chi_r
Profile at 850 um: 24 32.308 1.304
Profile at 350 um: 23 22.259 1.112
SED: 6 20.819 4.563
Total: 53 75.385 1.253
Best fit parameters and errors
Par: beta T_0 rho_0 rc/r0 qd
Val: 0.887 102.828 1.118E-16 20.192 1.743
Err: 0.061 4.659 1.327E-17 5.837 0.023
Derived parameters and errors
rho_c= 5.928E-19 +/- 3.078E-19
Model fit
Par: beta T_0 rho_0 rc/r0 qd chi_r
Val: 0.887 102.828 1.118E-16 20.192 1.743 1.253
Rnge 0.084 3.361 1.681E-18 0.168 0.084
Iteration parameters: Niter Nloop LoopFct
20 5 0.700
Date: 2014/12/04 Time: 16:51:36
List of routines
List of routines
Main, initialization
Main, initialization
initial (subroutine)
- called from: main
- makes calls to: read_profile, read_sed, read_fitpar
mainmenu (subroutine)
- called from: main
- makes calls to: modcalc, shortcalc, modelfit, write_fitpar
Model computation
Model computation
modcalc (subroutine)
- called from: main, mainmenu
- makes calls to: adaptgrid, envelopemass, inigraf, intensity, profile, graficpro, sedcalc, flux_rprofile, write_sed_obs, compare_sed, graficsed
shortcalc (subroutine)
- passed as rms_calcul to: lhalton_seed_fit, uncertainty, rms_estimate
- called from: mainmenu, halton_seed_fit, rms_estimate
- makes calls to: adaptgrid, intensity, profile, flux_rprofile, compare_sed
plotcalc (subroutine)
- passed as best_calcul to: halton_seed_fit
- called from: halton_seed_fit
- makes calls to: adaptgrid, intensity, profile, flux_rprofile, compare_sed
Parameters fit
Parameters fit
modelfit (subroutine)
- called from: mainmenu
- makes calls to: lhalton_seed_fit, uncertainty
halton_seed_fit (subroutine)
- called from: modelfit
- makes calls to: halton, rms_calcul (shortcalc), best_calcul (plotcalc)
uncertainty (subroutine)
- called from: modelfit
- makes calls to: rms_estimate
rms_estimate (subroutine)
- called from: uncertainty
- makes calls to: rms_calcul (shortcalc)
halton (function)
- called from: fitall
Intensity radial profile
Intensity radial profile
intensity (subroutine)
- called from: shortcalc, modcalc, sedcalc
column_density (subroutine)
- called from: write_herschel
profile (subroutine)
- called from: shortcalc, modcalc
- makes calls to: profiletomap, beamsmooth, maptoprofile, beamchop, flux_rprofile, beamprofile, compare_profile, write_profile, profiletomap, exportfits,
profiletomap (subroutine)
- called from: profile
beamprofile (subroutine)
- called from: profile
maptoprofile (subroutine)
- called from: profile
beamchop (subroutine)
- called from: profile
beamsmooth (subroutine)
- called from: profile
- makes calls to: conv2g
compare_profile (subroutine)
- called from: profile
SED
SED
sedcalc (subroutine)
- called from: modcalc
- makes calls to: intensity, flux_rprofile
compare_sed (subroutine)
- called from: shortcalc, modcalc
Small calculations
Small calculations
rho (function)
- called from: intensity, column_density
write_herschel (subroutine)
- called from: modcalc
- makes calls to: column_density
flux_rprofile (subroutine)
- called from: shortcalc, modcalc, profile, sedcalc
envelopemass (subroutine)
- called from: modcalc
beam_solidangle (function)
- called from: initial, read_profile, profile
freefree (function)
- called from: shortcalc, modcalc, sedcalc
Utilities, special functions
Utilities, special functions
jwl_um (function)
- called from: profile, write_profile, graficpro, exportfits
adaptgrid (subroutine)
- called from: shortcalc, modcalc
expmu (function)
- called from: intensity
rounded (function)
- called from: adaptgrid
gammln (function)
- called from: modcalc
Reading/writing files
Reading/writing files
read_fitpar (subroutine)
- called from: initial
write_fitpar (subroutine)
- called from: mainmenu
read_profile (subroutine)
- called from: initial
write_profile (subroutine)
- called from: profile
read_sed (subroutine)
- called from: initial
write_sed_obs (subroutine)
- called from: modcalc
write_herschel (subroutine)
- called from: modcalc
Graphic routines
Graphic routines
inigraf (subroutine)
- called from: main, modcalc
graficsed (subroutine)
- called from: modcalc
graficpro (subroutine)
- called from: modcalc
2D Gaussian convolution
2D Gaussian convolution
conv2g (subroutine)
- called from: beamsmooth
- makes calls to: rfft2, ftgauss, hfft2
ftgauss (subroutine)
- called from: conv2g
Fourier transform
Fourier transform
rfft2 (subroutine)
- called from: conv2g
- makes calls to: four1, fxrl1
hfft2 (subroutine)
- called from: conv2g
- makes calls to: four1, fxrl1
four1 (subroutine)
- called from: rfft2, hfft2
fxrl1 (subroutine)
- called from: rfft2, hfft2
FITS file writing
FITS file writing
exportfits (subroutine)
- called from: profile
- makes calls to: headerfits, writefits
headerfits (subroutine)
- called from: exportfits
writefits (subroutine)
- called from: exportfits
- makes calls to: testendian, swap
testendian (subroutine)
- called from: writefits
swap (subroutine)
- called from: writefits
Report abuse