Identifier syn_phot Purpose To produce synthetic photometry on IS0 SWS/LWS spectra or other spectra derived within the ISAP environment. This routine will allow comparison of integrated or averaged SEDs to other measurements such as from ISOPHOT, KAO observations, ground-based telescope observations, IRAS, or theoretical predictions. Synopsis struc_out = syn_phot(struc_in, filter, faar, saar, flux_int, err_stat, mode=mode, units_out_in=units_out_in, units_out=units_out) Arguments Name 1/0 Type Description struc_in I YAAAR Input structure (spectrum) containing the x-scale (wavelength) and the y-scale (flux density in Jy, W/m^2/um, or W/cm^2/um) that is to be used for synthetic photometry. filter I Name of filter which contains transmittance vs. wavelength. This name defines either an ASCII file or else a structure. faar O Filter structure in AAR format (possibly rebinned, see Comment 3 below) saar O Source Spectrum in AAR Format (possibly rebinned, see Comment 3 below) flux_int When mode = 'ave_nufnu' this is the integrated in-band flux returned by the procedure in W/cm^2 err_stat I Error status: 1 for fatal, 0 otherwise mode I The three choices are: 1. "ave" where the filter is used as a weighting function to compute the average flux density (making no assumptions about the spectral shape of the source) 2. "ave_nufnu" where the filter is used to compute the in-band flux and a flux density is returned assuming the source shape is nu x f_nu = const , or equivalently lambda f_lambda = constant 3. "int" where where the filter is used to compute the in-band flux default mode = "ave" units_in I allowed input units_out for the Source Spectrum are: `jy', `wcm2um', `wm2um', units_out default = `jy' units_out I For the "ave" and "ave_nufnu" modes the allowed output units_out are `jy', `wcm2um', `mag'. For the "int" mode the allowed units_out are w/cm2, w/m2 if units_out = `mag' the flux density is computed. A zero point value is found in the filter function directory area. The flux density is converted to magnitudes using the zero point and mag = 2.5 * LOG (zero_point/flux_density) units_out_out default = `jy' struc_out O struct YAAAR Structure containing the single flux or flux density or mag measurements vs. the single effective wavelength Returns A valid YAAAR structure containing a single flux, flux-density or a magnitude measurement at a single wavelength that represents the filter center. It returns the filter and input spectrum, with one or these rebinned onto the x axis or the other. For the ave_nufnu mode, the in-band flux is returned. Description Performs synthetic photometry on IS0 SWS or SWS spectra. Multiplies and integrates SWS and LWS spectra with filter functions to compute simulated in-band flux, or flux density, or magnitude Comments - The routine will operate on SWS, LWS, and PHT-S AARs (once PHT-S AAR are read into the YAAAR format). - This version will not make any adjustments for beam size and Source/Beam coupling. - The user may supply their own filter or else read a filter function out of the filter library using "read filter.pro", which will also check for the presence of zero-point values. The user must set up an environment variable (this is done in the standard ISAP start-up script) "SAP_FILTERS", which points to the directory containing filter files. - The input filter y-axis values are unitless transmissions which we call T'. In the case of the four IRAS filters, T' are the "Relative System Responsivities" - See The IRAS Explanatory Suppl. Section II.C.5. The input spectrum struct_in has a y-axis flux density value called S'. The valid flux densities are Jy, W/cm2/um, and W/m2/um. Values are internally converted to W/cm2/um. The routine does not yet propagate uncertainty values in its flux determination, although this is a high priority requirement Details of processing: 1) The input SWS/LWS spectrum or modeled SED spectrum is called S' and is also a function of lambda. Its flux density units_out are converted to W/cm2/um by the routine. 2) The filter is identified by name. If the name is the name of a structure, the routine uses this structure. Otherwise, the routine looks for the filter in an ISAP directory defined by an environment variable "SAP_FILTERS" The routine uses filter.pro to locate the filter and make YAAAR which is a function of lambda and has in place of its flux values, the filter's transmission. T'. 3) The routine examines the sampling of the wavelengths in the S' and filter T' arrays. It records the maximum delta_lambda in T' and the maximum delta_lambda in the wavelenguh range of S' The array has the smaller maximum delta_lambda is defined as the reference array for rebinning the other. The rebin is done. Now either T' is rebinned to T or S' is rebinned to S. We will hereafter refer to the two arrays that are on the same wavelength scale as T and S. In the rebinning process there is no extrapolation: beyond their range, each array is assumed to have y-values equal to zero. 4) Verify each S and filter are single-valued spectra with chk_multix using epsilon of 1e-5 microns. 5) Convert the S array to W/cm2/um. If the S array does not have units_in of flux density, exit with an error message. 6) When mode = "ave", normalize the T-values of the (possibly rebinned) filter function so that summation(1 to N) ( T x delta_lambda ) = 1 (See Summation Details Below) The routine will give a warning if any T < 0 on input and will set these values to 0 7) If mode = `ave' we compute the flux density by averaging S between the start and end wavelengths of the filter, weighted by T: flux density [W/cm2/um] = summation(1 to N) ( T x delta_lambda x S The mean wavelength in this mode is mean_lambda = summation(1 to N) ( T x delta_lambda x lambda) (See Summation Details Below.) 8) If mode = `int' we compute the flux treating the filter transmission values as actual transmission (e.g. a filter with all T values equal to 1 can be used integrate with unity transmission flux density over a given wavelength range to compute total flux.) the in-band flux is: Flux = summation(1 to N) (T delta_lambda S), W cm-2 or W m-2 The routine gives a warning if T < 0 or T > 1 on input. the mean wavelength in this mode is weighted by the transmission. mean_lambda = summation(1 to N) (T x wave) / summation (1 to N) (T) 9) If mode= `ave_nufnu' follow step 8 above to arrive at the flux in W/cm2. After this, we convert the in-band flux to flux density under the assumption that the source shape is of the form nu x F_nu = constant. To do so, we handle the four IRAS filters and all other filters slightly differently. We do this in order to remain faithful to published IRAS filter effective widths: (13.48, 5.16, 2.58, 1.00) x 10^12 Hz, and effective wavelengths (12, 25, 60, 100um). For other filters we derive delta_lambda_eff and lambda_eff. 10) For the IRAF filters, eg. filter = 'iras_12.dat', our mode='ave_nufnu' procedure is to convert the flux to W/m2 and then to Jy by dividing the filter widths - the delta_nu's in Hz listed above, and finally multiplying by 10^26. By definition the effective wavelengths of the IRAS filters are lambda_eff = 12, 25, 60, or 100um If the user has requested output units_out of W/m2/um, we convert to W/m2/um by multiplying Jy by 10^-12 / lambda_eff^2 If the user has requested magnitudes, we use the flux density in Jy and mag = 2.5 log (zero_point/flux_density) where the zero_point value (in Jy) is found in the SAP_FILTERS path. 11) If the filter is not an IRAS filter, we first compute the following Two quantities: lambda_eff = summation(1 to N) (t delta_lambda] / summation(1 to N) (t delta_lambda/lambda] delta_lambda_eff= summation(1 to N) (t delta_lambda) For the relation of t to T, see the summation description below. For the derivation of these two formulae see the Derivation below. We divide the in-band flux (step 8) by delta_lambda_eff to get W/m2/um or W/cm2/um. If the user has requested Jy as the preferred output then we convert W/cm2/um to Jy by multiplying by lambda_eff^2/2.9979e-16. For output in magnitudes we again use mag = 2.5 log (zero_point/flux_density). We report the wavelength as lambda_eff as above. 12) Derivation of lambda_eff and delta_lambda_eff for Non-IRAS Filters ------------------------------------------------------------------- The idea here is that we have an arbitrary filter, and we are using it to measure the flux of an SED of arbitrary shape. Call this SED the "actual" SED. We will report a flux density at a wavelength at the filter's center (center as defined as lambda_eff). If a hypothetical source had an SED with flux density of the form: lambda f_lambda = constant and had the flux_density which we report at lambda_effective, then such a source would have the same total flux within the filter which is observed from the _actual_ SED. We first derive lambda_eff: lambda_eff = summation(1 to N) lambda weight / summation(1 to N) weight. The weight we use at each point is the product of the transmission t, the wavelength interval, delta_lambda, and a hypothetical SED with lambda f_lambda = C, or f_lambda = C/lambda. Thus, weight = delta_lambda t C/lambda. So lambda_eff = summation(1 to N) (lambda t delta_lambda C/lambda) / summation(1 to N) ( t delta_lambda C/lambda) which reduces to: lambda_eff = summation(1 to N) (t delta_lambda) / summation(1 to N) (t delta_lambda/lambda) Now, what is the effective width of the filter, delta_lambda_eff? Assume lambda_eff has been computed as above. delta_lambda_eff is defined for an SED of the form f_lambda=C/lambda, such that: f_lambda(lambda_eff)= {summation(1 to N) (t delta_lambda f_lambda)} / delta_lambda_eff Substituting in f_lambda=C/lambda, C/lambda_eff = {summation(1 to N) (t delta_lambda C/lambda)} /delta_lambda_eff and noting that lambda_eff is given by the ratio of summations above, we arrive at the simple result: delta_lambda_eff = summation(1 to N) (t delta_lambda) 13) Summation Details ----------------- There are several summations refered to above. The source and filter wavelength scale are not necessarily uniformly spaced so we have established a specification as to how the summations will be made to handle this. Filters will be assumed to end at the last specified lambda. We will take summations between the first and last lambda using mid-point lambda and Y values: S S ----s---- S S --s-- -s- S S S S ------------------------------- ----+------+-------+------ l d l d l d l If the S array (after rebinning) looks as shown above on the left, and we wish to integrate over it, we create the "s" array which is has the mean S values "s" at the midpoint (eventually weighted by errorbars), and an associated set of delta_lambdas with midpoints at locations "d". The T array is similarly prepared in t, d pairs. The sums mentioned above become as follows: (T delta_lambda S) becomes (t delta_lambda s), (T delta_lambda wave) becomes (t delta_lambda d). 14) STRUC_OUT (YAAAR) HEADER/HISTORY : The primary header of output YAAAR structure contains the following keywords. SIMPLE = T / Standard FITS format BITPIX = 8 / number of bits per data pixel NAXIS = 0 / number of data axes EXTEND = T / FITS dataset may contain extension ORIGIN = 'ISAP/IPAC' / Data written by MAKEYAAARHEAD COMMENT = 'ISO YAAAR' END The secondary header of output YAAAR structure contains the following keywords. XTENSION= 'BINTABLE' / Binary table FITS extension BITPIX = 8 / 8 BITS character format NAXIS = 2 / Tables are 2-d character array NAXIS1 = 16 / Characters in a row NAXIS2 = 1 / Number of rows in a table PCOUNT = 0 / Parameter count always 0 GCOUNT = 1 / Group count always 1 TFIELDS = 4 / No of columns in table TFORM1 = '1E ' / TTYPE1 = 'WAVE ' / wavelength of data point TUNIT1 = 'um ' / TFORM2 = '1E ' / TTYPE2 = 'FLUX ' / flux TUNIT2 = 'units_out ' / TFORM3 = '1E ' / TTYPE3 = 'STDEV ' / standard deviation TUNIT3 = 'units_out ' / TFORM4 = ''1J ' TTYPE4 = ''FLAG ' / masking flag TUNIT4 = '' ' END The history tag will contain the history from input structure struc_in, history from the filter YAAAR (read_filter.pro) , and from syn_phot itself. 15) Error (Uncertainty) propagation - TBD [We need a lot of different examples still. I would like to replace the one below because the IRAS filters will not typically be used like this.] Examples 1) out = syn_phot(sws01_aar, "IRAS_12um", units_out='jy') Multiply the IRAS 12um filter transmittance curve with an SWS01 spectrum and estimate the flux under the spectrum in this band-pass. User will get F_band in W/cm^2/um at a mean wave near 12 um. The Jy's at this wave are Jy= F_band x wave ^ 2 /2.99792E-16 Dependencies sap_rebin, read_filter, environment variable "SAP_FILTERS" pointing to filters directory. CALLED FROM: IA, GUI Category ISAP Filename syn_phot.pro Author Iffat Khan (irk@ipac.caltech.edu) Version 0.1 History 27. 7.95 requirements outline D. Lutz 16. 1.96 draft requirements J. Mazz, S. Lord 14. 2 96 more req.s (utypes) S. Lord 26. 2 96 even more req.s (units_out) S. Lord 06. 3 96 Design and Code I. Khan 06.05 96 Revision I. Khan 27.06 96 made an IRAS nufnu mode S. Lord 19.07 96 delivered code for v_2.1 irk 22.07 96 add an additinal param "stat" irk to write_hdrkey call irk