This page was created by the IDL library routine
mk_html_help. For more information on
this routine, refer to the IDL Online Help Navigator
or type:
? mk_html_help
at the IDL command line prompt.
Last modified: Fri Apr 18 12:56:58 2008.
function addrsp combine two response matrices to make a composite RSP and return a structure in the same format as RD_OGIP_RMF() steps through each row of an RMF and tacks on the array from a 2nd RMF, updating the N_GRP, N_CHAN, F_CHAN elements as appropriate. syntax rsp=addrsp(rm1,rm2,ea1=ea1,ea2=ea2,eps=eps,verbose=verbose) parameters rm1 [INPUT; required] response matrix structure as read in by RD_OGIP_RMF rm2 [INPUT; required] response matrix structure to be added to RM1 * RM1 and RM2 are expected to have the following fields: NNRG,ELO,EHI,NCHAN,EMN,EMX,N_GRP,F_CHAN,N_CHAN,MATRIX,FIRSTCHAN keywords ea1 [INPUT] if set to a scalar or size matches RM1.ELO, then multiplies RM1 by EA1 prior to addition ea2 [INPUT] as EA1, but for RM2 * default for EA1 and EA2 is 1.0 eps [INPUT] a small number, below which to set everything to zero * default is 1e-6 morecol [INPUT] designed as a buffer so that a matrix won't overflow the bounds automatically set. by default, the new matrix is defined at an early stage to be max(RM1.N_CHAN)+max(RM2.N_CHAN)+100+MORECOL. * note that MORECOL can even be negative. * use wisely. verbose [INPUT] controls chatter ea3 [JUNK] here just to prevent user error ea4 [JUNK] here just to prevent user error ea5 [JUNK] here just to prevent user error ea6 [JUNK] here just to prevent user error * if EA3, EA4, EA5, or EA6 are set, the program will complain and quit. this is done to avoid the obvious user error of trying to specify EA3=EA3, etc. if adding a third RMF/ARF to a preexisting one. _extra [JUNK] here only to avoid crashing the program subroutines KILROY OGIPZIP history vinay kashyap (Nov2001) bug correction re max possible matrix elements when NGROUP>>1; added dummy keywords EA3, EA4, EA5, EA6 and extra checks against tricky inputs (VK; Aug02) bug correction: wasn't handling N_GRP=1 (VK; Aug04)
(See pro/util/addrsp.pro)
procedure adjustie adjust the parameters to account for any ties among them, and return the corrected set syntax adjustie,a,ties=ties,vname=vname parameters a [I/O; required] parameters to be updated keywords ties [INPUT] string array describing how a parameter is tied to other parameters. * must be a mathematical expression of the form "a0 = a1/10.", "a2=a3+a4", "a5 = a5 < 0", etc. vname [INPUT; default='a'] in case TIES contains variable names such as "p0=p2^2", etc, then set VNAME='p' NOTE: the following keywords are generated in situ the first time the function is called. on subsequent calls they are used instead of EXECUTEing TIES in order to improve speed. NOTE2: actually, they are not implemented yet. someday. ifit [I/O] integer array describing which parameters are frozen (0's) and which are not (1's). updated upon exit to include information given in TIES itie [I/O] integer array showing which parameters are tied. tieto [I/O] linked list showing which parameters each of the ITIE are tied to which parameter tieops [I/O] as TIETO, but containing information on which operator is used, e.g., "+", "-", "*", "/", "<", ">", etc. _extra [JUNK] here only to prevent crashing the program history vinay kashyap (Oct98) converted to IDL5 (VK; OctMM)
(See pro/fitting/adjustie.pro)
procedure apedance remove or correct the abundance dependance of APED emissivities APED emissivities as stored on disk include both ion balance (Mazzotta et al.) and abundances (Anders & Grevesse). This can be inconvenient in some PINTofALE tasks since other line emissivity databases in PINTofALE do not include either factor. It is not possible to remove the effects of the included ion balance, but removing the abundance dependance is simply a matter of dividing the emissivity of each line by the abundance appropriate to the element producing that line. syntax apedance,line,Z,abund=abund,apabref=apabref,verbose=verbose parameters line [I/O] APED line emissivities * emissivity structure out of RD_LINE() * if 2D array, assumed to be an array of size (T,Wvl) Z [INPUT] atomic numbers * required if LINE is an array and not a structure * must match size of 2nd dimension in LINE keywords abund [INPUT] if given, updates LINE by multiplying by the ratio of ABUND/(Anders & Grevesse) * the default is to just divide by Anders & Grevesse abundances, effectively removing the abundance dependance from the emissivities apabref [INPUT] if set to a reference other than Anders & Grevesse, it is assumed that the input LINE emissivities need to be corrected according to the specified abundances. * WARNING: do not set this unless you know exactly what you are doing verbose [INPUT] controls chatter _extra [JUNK] here only to avoid crashing the program warning no checks are made to verify that the emissivities have have already had the abundances taken out, or that the abundances are corrected according to some other baseline. USE CAREFULLY! history vinay kashyap (Apr2004)
(See pro/apedance.pro)
function arithtogram return the result of an arithmetical operation on the frequency distributions of two lists syntax hx=arithtogram(x1,x2,xout,operator,/plus,/minus,/divide,$ w1=w1,w2=w2,nbin=nbin,xmin=xmin,xmax=xmax,/nonorm) parameters x1 [INPUT; required] list of numbers whose frequency distribution must be multiplied with that of X2 x2 [INPUT; required] list of numbers whose frequency distribution must be multiplied with that of X1 * it is assumed that X1 and X2 both have the same units because otherwise this kind of multiplication makes no sense xout [OUTPUT; required] the list of numbers at which the output is tabulated * by default, this is the sorted unique set of [X1,X2] * if NBIN is set, generates a logarithmic or linear grid of that size oper [INPUT] must be one of '+', '-', '/', '*': the arithmetical operation to carry out on the frequency distributions of X1 and X2, i.e., result = X1 OPER X2 * if not given, assumed to be '*', unless one of keywords PLUS, MINUS, or DIVIDE are set keywords plus [INPUT] if set, OPER is assumed to be '+' minus [INPUT] if set, OPER is assumed to be '-' divide [INPUT] if set, OPER is assumed to be '/' * DIVIDE takes precedence over MINUS takes precedence over PLUS w1 [INPUT] optional weight to be applied to histogram(X1) w2 [INPUT] optional weight to be applied to histogram(X2) * if not given, W1 and W2 are assumed to be 1 nbin [INPUT] number of bins in the output * if not set, will be the sorted unique set of [X1,X2] * if -ve, will produce logarithmic gridding xmin [INPUT] minimum value in the output grid * by default, uses min(X1)>min(X2) xmax [INPUT] maximum value in the output grid * by default, uses max(X1)<max(X2) nonorm [INPUT] if set, does not make a correction for the width of the bins in XOUT. _extra [JUNK] here only to prevent crashing the program description this program does not create histograms as intermediate products but rather produces an operation at the best available data resolution. this is of great use in Monte Carlo calculations or in creating probability distributions in MCMC. the algorithm is straightforward: build up a cdf for each list, and interpolate onto a common grid, and multiply, divide, add, or subtract the d(cdf)'s, suitably normalized for the bin widths and number of elements. example x1=randomn(seed,10000L) & x2=randomn(seed,10000L)+2. xmin=-2. & xmax=3. & nbin=101L hh=arithtogram(x1,x2,xout,'*',xmin=xmin,xmax=xmax,nbin=nbin,/nonorm) hp=arithtogram(x1,x2,xout,'+',xmin=xmin,xmax=xmax,nbin=nbin,/nonorm) hd=arithtogram(x1,x2,xout,'/',xmin=xmin,xmax=xmax,nbin=nbin,/nonorm) hm=arithtogram(x1,x2,xout,'-',xmin=xmin,xmax=xmax,nbin=nbin,/nonorm) h1=histogram(x1,min=xmin,max=xmax,binsize=median(xout[1:*]-xout)) h2=histogram(x2,min=xmin,max=xmax,binsize=median(xout[1:*]-xout)) plot,xout,hh,psym=10 & oplot,xout,h1,col=2 & oplot,xout,h2,col=3 oplot,xout,h1*h2,psym=10,col=4 plot,xout,hp,psym=10 & oplot,xout,h1,col=2 & oplot,xout,h2,col=3 oplot,xout,h1+h2,psym=10,col=4 plot,xout,hd,psym=10 & oplot,xout,h1,col=2 & oplot,xout,h2,col=3 oplot,xout,float(h1)/float(h2),psym=10,col=4 plot,xout,hm,psym=10 & oplot,xout,h1,col=2 & oplot,xout,h2,col=3 oplot,xout,h1-h2,psym=10,col=4 history vinay kashyap (Jun2005)
(See pro/util/arithtogram.pro)
function arrayeq checks whether two input arrays are equal, and returns 1: if they are identical 0: if at least one element differs past tolerance -1: if sizes differ, but one array is a simple subset of the other syntax AeqB=arrayeq(A,B,tol=tol) parameters A [INPUT; required] first array B [INPUT; required] second array * arrays must be of some >numerical< type keywords tol [INPUT; default: 1e-4] tolerance within which two numbers are said to be identical _extra [JUNK] here only to avoid crashing the program history vinay kashyap (Nov98) converted to IDL5 notation (VK; OctMM)
(See pro/util/arrayeq.pro)
function bamabs return photoelectric absorption cross-sections in energy range 30 eV - 10 keV using the polynomial fit coefficients determined by Monika Balucinska-Church and Dan McCammon (Balucinska-Church, M.\ and McCammon, D.\ 1992, ApJ 400, 699). This is an update to Morrison & McCammon 1983. syntax sigabs=bamabs(w,abund=abund,/ikeV,/Fano,verbose=v) parameters w [INPUT; required] photon energy values at which to compute the absorption * default units are [Ang], unless IKEV is set, when they are assumed to be [keV] keywords abund [INPUT] abundances relative to H=1 * default is to use Anders & Grevesse 1982 ikeV [INPUT] if set, assumes that W are in units of keV Fano [INPUT] if set, the 4 strongest auto-ionizing resonances of He I are included; numbers come from Oza (1986, Phys Rev A, 33, 824), and Fernley et al. (1987, J.Phys B 20, 6457). noHeH [INPUT] if set to any number other than 1 or 2, excludes H and He from the cross-sections * if set to 1, only excludes H * if set to 2, only excludes He verbose [INPUT] controls chatter _extra [JUNK] here only to avoid crashing the program subroutines GETABUND INICON restrictions works only in energy range 30 eV - 10 keV be warned that the edges are idealized and thus unrealistic for high resolution spectra because they do not include any resonance structure history vinay kashyap (OctMM; based on TOTLXS.FOR and XSCTNS.FOR of Balucinska-Church & McCammon, obtained from ADC catalog VI/62A, ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/6/6062A/ ) bug corrections (Brian Kern, 5Dec2001): Ne cross-section was missing "^3"; Fano calc had EPSj[j] (VK; Dec2001) force calculations only in valid energy range; added keyword NOHEH (VK; Jun03) bug correction: "exclude=0" means "include" (VK; Jul03) bug correction: lambda selection for He (Andy Ptak; VK Jun05) edge lambda correction: O, N, C changed acc. to numbers in Atomic Data Booklet (2nd Edition), Thompson, A., et al., 2001, LBL/PUB-490 Rev. 2 (Jeremy Drake; VK Sep2006)
(See pro/bamabs.pro)
function BEHR_hug
IDL wrapper to Bayesian Estimate of Hardness Ratios (BEHR)
returns a structure containing the relevant outputs neatly
summarized into fields for IDL consumption
references:
"Bayesian Estimation of Hardness Ratios: Modeling and Computations",
Park, T., Kashyap, V.L., Siemiginowska, A., van Dyk, D., Zezas, A.,
Heinke, C., and Wargelin, B., 2006, ApJ, 652, 610
"BEHR: Bayesian Estimation of Hardness Ratios", Park, T., Kashyap, V.,
Zezas, A., Siemiginowska, A., van Dyk, D., Connors, A., and Heinke, C.,
2005, at Six Years of Science with Chandra Symposium, Nov 2-5, 2005, #5.6
"Computing Hardness Ratios with Poissonian Errors",
van Dyk, D.A., Park, T., Kashyap, V.L., \& Zezas, A.,
2004, AAS-HEAD #8, 16.27
http://www.aas.org/publications/baas/v36n3/head2004/137.htm
syntax
behr=behr_hug(softsrc,hardsrc,softbkg,hardbkg,softarea,hardarea,$
softeff=softeff,hardeff=hardeff,softidx=softidx,hardidx=hardidx,$
softscl=softscl,hardscl=hardscl,softtbl=softtbl,hardtbl=hardtbl,/post,$
level=level,algo=algo,details=details,nsim=nsim,nburnin=nburnin,$
/hpd,nbins=nbins,outputf=outputf,outputR=outputR,outputHR=outputHR,$
outputC=outputC,outputMC=outputMC,outputPr=outputPr,BEHRdir=BEHRdir,$
verbose=verbose)
parameters
softsrc [INPUT; required] counts in source region in soft (S) band
hardsrc [INPUT; required] counts in source region in hard (H) band
softbkg [INPUT; required] counts in background region in S band
hardbkg [INPUT; required] counts in background region in H band
softarea[INPUT; required] background scaling factor in S band
hardarea[INPUT; required] background scaling factor in H band
* (background region area)/(source region area)
* can also include differences in exposure time into
the ratio, in the same manner as geometric area
keywords
softeff [INPUT] effective area in soft (S) band
* if not set, then set equal to HARDEFF if given, or
1 otherwise
hardeff [INPUT] effective area in hard (H) band
* if not set, then set equal to SOFTEFF if given, or
1 otherwise
* can also be the effective area relative to some
special point on the detector (e.g., aimpoint)
or even some specific detector (e.g., ACIS-I v/s ACIS-S)
softidx [INPUT] index of prior on S (range = 0+)
* if not set, then set equal to HARDIDX if given, or
0.5 otherwise
hardidx [INPUT] index of prior on H (range = 0+)
* if not set, then set equal to SOFTIDX if given, or
0.5 otherwise
* similar to AGAMMA of PPD_SRC()
softscl [INPUT] scaling index of prior on S
* if not set, then set to HARDSCL if given, or
0 otherwise
hardscl [INPUT] scaling index of prior on H
* if not set, then set to SOFTSCL if given, or
0 otherwise
* similar to BGAMMA of PPD_SRC()
softtbl [INPUT] filename containing a tabulated prior for S
hardtbl [INPUT] filename containing a tabulated prior for H
* the table prior must be an ascii file with the following format:
line 1: number of entries, say NLIN
line 2: labels for the columns, ignored
lines 3..NLIN+2: two whitespace separated columns of numbers,
with each row containing the source intensity and the posterior
density, in that order
* the default filenames are "./tblprior_{soft|hard}.txt"
* the default filenames are used iff SOFTTBL and HARDTBL are set
but do not appear to exist
* WARNING: if regex is used in the filename specification, only the
first file from the list will be used. furthermore, if specified,
the table priors are applied to _all_ SOFTSRC and HARDSRC
post [INPUT] if set, suggests the values of (SOFTIDX,SOFTSCL)
and (HARDIDX,HARDSCL) going forward, i.e., what you should
set the priors to in your next calculation for the same
source
level [INPUT] percentage confidence level at which to report error
(default = 68)
details [INPUT] compute various ratios (true/false)?
(default = true)
algo [INPUT] calculation method, GIBBS (default) or QUAD
* if set to AUTO or AUTO=N, then uses GIBBS for any
case where {SOFT|HARD}SRC > N, and QUAD below that
unless one of {SOFT|HARD}SRC > 99, in which case
GIBBS is set automatically
- if N is not set, assumed to be 15
- N can be set to any integer less than 100
nsim [INPUT] number of draws if algo=gibbs (default=10000)
nburnin [INPUT] number of burn-in draws if algo=gibbs
(default=5000 or NSIM/2, whichever is smaller)
nbins [INPUT] number of bins in integration if algo=quad (default=500)
hpd [INPUT] if set, computes the highest posterior density interval,
which also is the smallest interval that includes the mode.
* this is set by default
* only has an effect if ALGO=QUAD
outputF [INPUT] root of filename in which to place output
(default = 'none')
* output will be placed in the file, OUTPUTF.txt
* unless OUTPUTF='none', in which case nothing is written out,
unless one of OUTPUT{R|C|HR|MC|PR} are set, in which case
the corresponding output is put in file BEHR_{R|C|HR|MC|PR}.txt
outputR [INPUT] if set, writes output for R in OUTPUTF_R.txt
outputC [INPUT] if set, writes output for C in OUTPUTF_C.txt
outputHR[INPUT] if set, writes output for HR in OUTPUTF_HR.txt
outputMC[INPUT] if set, writes Monte Carlo draws for lamS and lamH
to OUTPUTF_draws.txt when algo=gibbs
outputPr[INPUT] if set, writes the probability distributions for
R, HR, C, lamS, and lamH to OUTPUTF_prob.txt when algo=quad
BEHRdir [INPUT] full path to directory where BEHR executable resides
(default = '/data/fubar/kashyap/AstroStat/BEHR')
verbose [INPUT] controls chatter
requirements
BEHR executable must be installed in BEHRDIR
BEHR should be executable under the shell via SPAWN
BEHR output assumed to be compatible with 12-19-2005 version
side-effects
potentially creates numerous ascii files in $cwd
history
vinay kashyap (SepMMV)
bugfix: output file format was overflowing; force NBURNIN < NSIM;
outputF for 1st iteration; added more references (VK; NovMMV)
added keywords SOFTSCL,HARDSCL,POST,HPD; changed DETAILS default
to true, OUTPUTF default to none; added extra fields to output
structure; added ALGO option AUTO (VK; DecMMV)
altered behavior of ALGO=AUTO to use GIBBS if counts in one of the
bands exceeds 100 (VK; OctMMVI)
added keywords SOFTTBL and HARDTBL (VK; FebMMVIII)
(See pro/stat/behr_hug.pro)
function binErsp
returns position indices of input energ(y/ies) appropriately
binned into the specified binning scheme
returns -1 where E is outside the bin range
syntax
iE=binersp(nrg,binmn,maxbin,bstr=bstr,nbin=nbin)
parameters
nrg [INPUT; required] photon energies for which to find
the binning numbers
binmn [INPUT] array of bin minimum values
* overrides all keywords
* if not given, see description below on how it is
derived from BSTR.ELO, and NBIN & NRG.
binmx [INPUT] maximum of the bin max values.
* if not a scalar, then the last element of array is used
* see description below for how it is determined if
not given.
keywords
bstr [INPUT] structure containing the following pieces of
information, as obtained from an OGIP-compliant response
matrix (see RDRESP.PRO/RDARF.PRO): {NNRG, ELO, EHI}
nbin [INPUT; default=101] number of equally spaced bins into
which to split the input array of photon energies
* if -ve, binning will be logarithmic!
_extra [JUNK] here only to prevent crashing the program
description
if BINMN is not given, then BINMN=BSTR.ELO
if BSTR is also not given (or is incorrect), then
BINMN=FINDGEN(ABS(NBIN))*DE+EMIN
where DE=(EMAX-EMIN)/(ABS(NBIN)-1), EMIN=MIN(EE,MAX=EMAX),
and EE=NRG if NBIN>0, EE=ALOG10(NRG) if NBIN<0
if BINMX < MAX(BINMN), BINMN is truncated
if BINMX < MIN(BINMN), BINMX is reset to MIN(BINMN)
if BINMX is not given, then BINMX=MAX(BINMN)+BINMX(LAST)-BINMX(LAST-1)
if BINMN is not given, BINMX=MAX(BSTR.EHI)
if BSTR is also not given (or is incorrect), then BINMX=EMAX
history
vinay kashyap (Jul97)
converted to IDL5 notation (VK; OctMM)
(See pro/util/binersp.pro)
function bland returns a "blending factor" for each line in the emissivity list, defined as the ratio of the flux in the line to that within the given range. (i.e., small values mean heavily blended lines) syntax b=bland(flstr,dwvl,flx=flx,wflx=wflx,arstr=arstr,$ dem=dem,abund=abund,effar=effar,wvlar=wvlar,$ /temp,/noph,/ikev,/regrid) parameters flstr [INPUT; required] structure containing emissivity and wavelength information. see RD_LINE for details. * ion-balance assumed to be included * if FLX and WFLX are set on input, FLSTR is ignored -- still required, but ignored. dwvl [INPUT; default=median(delta(FLSTR.WVL))] the spectral "resolution" -- any line nearer than this value at any given wavelength is considered to be blended in to this. * if array of size FLSTR.WVL (or twice that size), set to different values (with possible asymmetric ranges) at each wavelength. keywords arstr [INPUT] structure of structures containing effective area information (see ARIA or GRATFLX) * if given in the right format, overrides call to LINEFLX with call to GRATFLX flx [I/O] contains the calculated flux due to each line * may be given on input and if WFLX is also set, FLSTR and ARSTR are ignored wflx [INPUT] if size matches that of FLX, assumed to be the wavelengths of the lines. * if FLX and WFLX are set on input, FLSTR and ARSTR are ignored. _extra [INPUT] pass defined keywords to GRATFLX (DEM, ABUND, TEMP, NOPH, IKEV, REGRID) LINEFLX (DEM, ABUND, EFFAR, WVLAR, TEMP, NOPH, IKEV, REGRID) subroutines ARIA GRATFLX LINEFLX history vinay kashyap (Sep98) bug fix for when FLSTR is ignored (VK; Oct98) added keyword ARSTR, added call to GRATFLX (VK; Nov98) converted to IDL5 notation (VK; OctMM) bug fix: crashed if FLX and WFLX set, but DWVL not set (VK;JanMMI) bug fix: FLSTR can have more than 8 tags (JJD; MayMMVII)
(See pro/bland.pro)
function cat_id
returns a concatenated structure of Line IDs
syntax
C=cat_id(A,B,pick=pick,ask=ask,/okev,comm=comm,flst=flst,dbdir=dbdir,$
fluxes=fluxes,fluxerr=fluxerr)
parameters
A [INPUT; required] output of LINEID
B [INPUT] output of LINEID to be merged with A
* if any wavelengths are common, the IDs in A take
precedence (but see keyword ASK)
* if not given, then the program *prints* a summary
of the IDs to the screen
keywords
pick [INPUT; default: all] position indices of elements in A
to be selected (or not).
* if scalar, returns all elements EXCEPT the specified one
* has no effect on B
* use these two to delete entries, for e.g.
ask [INPUT] if set, asks for confirmation before throwing
away an ID from B that also exists in A. allows overriding
the default action.
* ASK='n' will reverse the default by throwing away the
ID from A instead of the one from B
* ASK='r' *reverses* the default, in that if there is a
component in B that's also in A, the one in A is >replaced<
by the one in B (as opposed to if ASK='n', then the one in
A is deleted, and B is just appended at the end)
* ASK='k' overrides the default by keeping BOTH
okev [INPUT] if set, converts wavelengths [Ang] to energy [keV]
comm [INPUT] appends a comment to each line while printing
* takes the comments from A.LABL(1,i)+A.LABL(0,i)
* truncates the comment to a maximum of 30 characters
unless COMM is set to a higher number
* for this to work when FLST is set, make sure DBDIR is
also defined
flst [INPUT] if set, prints out the line list in the form that
RD_LIST likes
* if FLST is a string, writes to file FLST
mplet [OUTPUT] if set, will contain on output a 2 dimensional MPLET
string array as prescribed by MIXIE()
dbdir [INPUT] if set and FLST is set, appends DBDIR to IDs
* 1 => $SCAR , 2 => $SPEX , 3 => $CHIANTI
fluxes [INPUT] if set and matches either the number of components
or the number of matches, then renormalizes/overwrites the
line fluxes in the output.
* ONLY if IDs are being printed
fluxerr [INPUT] if set and matches either the number of components
or the number of matches, then appends an error to output
* ignored if FLST is not set
* if FLUXERR is given for each match, square-adds the errors
for the components
* if single element, but there are more than 1 component to
ID structure, this is assumed to be the fractional error
(if < 1) or S/N (if > 1) for each match
_extra [JUNK] here only to prevent crashing the program
restrictions
* works ONLY on the outputs of LINEID
* requires subroutines:
-- INICON
-- CREATE_STRUCT
-- LINEID
history
vinay kashyap (Feb97)
bug -- was breaking down in simplest case (VK; Mar97)
was breaking down if ID was unknown (VK; Jul97)
added keyword FLST (VK; Oct98)
added call to INITSTUFF; added Tmax to output; added keywords DBDIR,
FLUXES, FLUXERR; corrected bug re undefined output when FLST was set;
increased number of significant digits in FLST output; generalized
ASK to "Yes/No/Replace/Keep" (VK; Nov98)
updated FLUXERR behavior to account for new field (VK; Dec98)
modified to allow for field NOTES to be printed out (VK; Mar99)
changed call to INITSTUFF to INICON (VK; 99May)
added FLUXERR to output (VK; MarMM)
converted to IDL5 array notation (VK; OctMM)
changed keyword KEV to OKEV (VK; JanMMI)
now /COMM also works when FLST (and DBDIR) are defined (VK;Nov01)
added MPLET keyword (LL; Feb04)
(See pro/cat_id.pro)
function cat_ln returns a concatenated structure of Line emissivities syntax C=cat_ln(A,B,pick=pick,ask=ask,comm=comm,reord=reord,/okev,$ abund=abund,/noph,effar=effar,wvlar=wvlar,dem=dem,/temp,/ikev,/regrid) parameters A [INPUT; required] output of RD_LINE B [INPUT] output of RD_LINE to be merged with A * if any wavelengths are common, the IDs in A take precedence (but see keyword ASK) * if the temperature grids do not match, the emissivities in B are interpolated to the same grid as in A * if not given, then the program *prints* a summary of the lines to the screen keywords pick [INPUT; default: all] position indices of elements in A to be selected (or not). * if scalar, returns all elements EXCEPT the specified one * has no effect on B * use these two to delete entries, for e.g. ask [INPUT] specifies how to handle duplicate lines * if not set, doesn't care about duplicates and keeps all * if set (see exceptions below), will discard duplicates from B * ASK='n' will reverse the normal by throwing away the line from A instead of the one from B * ASK='r' will >replace< the emissivity, etc. in A with that of B * ASK='k' will act just as though ASK has not been set okeV [INPUT] if set, converts wavelengths [Ang] to energy [keV] in the screen listing comm [INPUT] appends a comment to each line while printing * takes the comments from A.DESIG[1,i]+A.DESIG[0,i] * truncates the comment to a maximum of 30 characters unless COMM is set to a higher number flst [INPUT] if set, prints out the line list in the form that RD_LIST likes * if FLST is a string, writes to file FLST reord [INPUT] if set, rearrange the order of the output * +ve: descending order of fluxes * -ve: ascending order of fluxes * if abs(REORD) .NE. 1, limit output to first REORD entries _extra [INPUT ONLY] use to pass defined keywords to LINEFLX: ABUND, NOPH, EFFAR, WVLAR, DEM, TEMP, REGRID restrictions * works ONLY on the outputs of RD_LINE * requires subroutines: -- LINEFLX [GETABUND, WHEE] -- RD_LINE -- CREATE_STRUCT -- INICON history vinay kashyap (Jun98, based on CAT_ID) changed keyword LNLST to FLST, added ABUND and REORD (VK; Jul98) delete keyword ABUND, and added call to LINEFLX (VK; Aug98) corrected bug that was skipping B because of roundoff error in first element; bug in handling ECONF (VK; Nov98) changed keyword KEV to OKEV to avoid conflict with LINEFLX keyword of same name (VK; 99Apr) added support for structure field JON (VK; 99May) allowed merging of mismatched LOGT grids (VK; 99Jul) converted to IDL5 array notation (VK; OctMM) added full RD_LIST compatible printout for FLST; changed output file wavelength format from f8.4 to f10.5 (VK; JanMMI) bug -- crashing on <empty> for FLST -- fixed (VK; Jul04) added extra columns (FLUX, logTMAX) to output of FLST (VK; Feb08)
(See pro/cat_ln.pro)
function chanarf compute the average effective area over each channel, as a weighted average of the effective area from all the energies that contribute to that channel. makes a difference mainly for low and medium spectral resolution data. for high-res grating data, not so much. syntax carf=chanarf(rmfstr,effar,spec=spec) parameters rmfstr [INPUT; required] response matrix structure, in the same format as returned from RD_OGIP_RMF() effar [INPUT] effective area as a function of energy [cm^2] * assumed to be on the same grid as the input energies of RMFSTR, i.e., RMFSTR.ELO and RMFSTR.EHI * ignored if size doesn't match RMFSTR keywords spec [INPUT] a spectrum to further weight the contribution of EFFAR * assumed to be on the same grid as RMFSTR.ELO and RMFSTR.EHI * ignored if size doesn't match EFFAR * ideally should be in units of [ph/keV/...] chnrg [OUTPUT] the average energy of a photon in this channel * it should be similar to 0.5*(RMFSTR.EMN+RMFSTR.EMX), and serves as a check on the gain calibration verbose [INPUT] controls chatter _extra [JUNK] here only to prevent crashing history vinay kashyap (MMVI.IX)
(See pro/chanarf.pro)
PROJECT : CHIANTI
CHIANTI is an Atomic Database Package for Spectroscopic Diagnostics of
Astrophysical Plasmas. It is a collaborative project involving the Naval
Research Laboratory (USA), the University of Florence (Italy), the
University of Cambridge and the Rutherford Appleton Laboratory (UK).
NAME : CH_GET_FILE
PURPOSE : to select a file from either a selected directory or the working
directory, having an extension.
EXPLANATION : a file in either a selected directory or the working
directory, having an extension can be selected using a
widget. Note that both directory and extension have to be
supplied. If no file is found, an empty string is returned. ;
USE : IDL> name = ch_get_file( '~/', '.pro', tit=' Select a procedure ')
EXAMPLES : dir= concat_dir(!xuvtop),'dem')
dem_name=ch_get_file(path=dir,filter='*.dem',title='Select DEM File')
INPUTS : directory, extension
OPT. INPUTS :
OUTPUTS : the file name
OPT. OUTPUTS:
KEYWORDS : title
CALLS : findfile, break_file
COMMON : co
RESTRICTIONS: both directory and extension have to be
supplied.
SIDE EFFECTS:
CATEGORY :
PREV. HIST. : extracted from CDS/CHIANTI routines.
WRITTEN :
Giulio Del Zanna (GDZ),
DAMTP (University of Cambridge, UK)
MODIFIED : Version 1, GDZ 10-Oct-2000
V.2, GDZ, corrected a typo at the end of the file.
V.3, GDZ, generalized directory concatenation to work for
Unix, Windows and VMS.
V. 4, 19-July-2002, GDZ
Added the option to select files also with the standard IDL
dialaog_pickfile, and changed a few things...
V.5, 2-Aug-02, GDZ
reduced the size of the widget.
V.6, 12-Aug-02, GDZ
corrected for a bug in the directory output.
V.7, 3-Nov-03 GDZ
Fixed a bug when using Windows, the returned path was not
correct.
VERSION : 7, 3-Nov-03
(See pro/external/ch_get_file.pro)
function colbehr
computes hardness ratios based on 3-band data by invoking
BEHR (Bayesian Estimate of Hardness Ratios) twice, and
returns the inputs as well as the following quantities
and their credible ranges in a structure:
{S, M, H, lS=log10(S), lM=log10(M), lH=log10(H),
SpM=S+M, MpH=M+H, SpH=S+H, T=S+M+H,
SmM=S-M, MmH=M-H, SmH=S-H,
R1=S/M, R2=M/H, R3=S/H,
C1=C_SM, C2=C_MH, C3=C_SH,
HR1=(S-M)/(S+M), HR2=(M-H)/(M+H), HR3=(S-H)/(S+H),
HRA=(S-M)/(S+M+H), HRB=(M-H)/(S+M+H), HRC=(S-H)/(S+M+H)}
Warning: Because of the necessity of combining two separate
runs of BEHR, only the MCMC option is used. Thus, if the
length of the chain is small, the computed values may be
subject to computational instability.
Warning: BEHR assumes that the passbands in question
do not overlap, and that the counts input to the program
are statistically independent. It is up to the users to
ensure the validity of this assumption. No checks are
made to verify it.
Reference:
"Bayesian Estimation of Hardness Ratios: Modeling and Computations",
Park, T., Kashyap, V.L., Siemiginowska, A., van Dyk, D., Zezas, A.,
Heinke, C., and Wargelin, B., 2006, ApJ, 652, 610
syntax
behr3=colbehr(Ssrc,Msrc,Hsrc,Sbkg=Sbkg,Mbkg=Mbkg,Hbkg=Hbkg,$
Sarea=Sarea,Marea=Marea,Harea=Harea,$
Seff=Seff,Meff=Meff,Heff=Heff,$
Sidx=Sidx,Midx=Midx,Hidx=Hidx,$
Sscl=Sscl,Mscl=Mscl,Hscl=Hscl,$
Stable=Stable,Mtable=Mtable,Htable=Htable,$
/post,level=level,nsim=nsim,nburnin=nburnin,hpd=hpd,$
outputf=outputf,BEHRdir=BEHRdir, verbose=verbose)
parameters
Ssrc [INPUT; required] counts in source region in the soft (S) band
Msrc [INPUT; required] counts in source region in the medium (M) band
Hsrc [INPUT; required] counts in source region in the hard (H) band
* can be arrays; if so, the array with the most elements
determines the size of the output and the shortfalls in the
others, if any, are made up by replicating the first elements
keywords
Sbkg [INPUT] counts in background region in the S band
Mbkg [INPUT] counts in background region in the M band
Hbkg [INPUT] counts in background region in the H band
* if not given, assumed to be 0
* if size smaller than Xsrc, first element gets replicated
Sarea [INPUT] background scaling factor in the S band
Marea [INPUT] background scaling factor in the M band
Harea [INPUT] background scaling factor in the H band
* (background region area)/(source region area)
* if not given, assumed to be 1
* can also include differences in exposure time into
the ratio, in the same manner as geometric area
* if size smaller than Xsrc, first element gets replicated
Seff [INPUT] effective area in S band
Meff [INPUT] effective area in M band
Heff [INPUT] effective area in H band
* if none are set, all are assumed to be 1,
else if one is set, all are assumed to be equal to that one,
else if two are set and unequal, third is assumed to be 1
* can also be the effective area relative to some
special point on the detector (e.g., aimpoint)
or even some specific detector (e.g., ACIS-I v/s ACIS-S)
* if size smaller than Xsrc, first element gets replicated
Sidx [INPUT] index of prior on S (range = 0+)
Midx [INPUT] index of prior on M (range = 0+)
Hidx [INPUT] index of prior on H (range = 0+)
* if none are set, all are assumed to be 0.5,
else if one is set, all are assumed to be equal to that one
else if two are set and are unequal, third is assumed to be 0.5
* if size smaller than Xsrc, first element gets replicated
* similar to AGAMMA of PPD_SRC()
Sscl [INPUT] scaling index of prior on Ssrc
Mscl [INPUT] scaling index of prior on Msrc
Hscl [INPUT] scaling index of prior on Hsrc
* if none are set, all are assumed to be 0
else if one is set, all are assumed to be equal to that one
else if two are set and are unequal, third is assumed to be 0
* if size smaller than Xsrc, first element gets replicated
* similar to BGAMMA of PPD_SRC()
Stable [INPUT] filename containing a tabulated prior for Ssrc
Mtable [INPUT] filename containing a tabulated prior for Msrc
Htable [INPUT] filename containing a tabulated prior for Hsrc
* the table prior must be an ascii file with the following format:
line 1: number of entries, say NLIN
line 2: labels for the columns, ignored
lines 3..NLIN+2: two whitespace separated columns of numbers,
with each row containing the source intensity and the posterior
density, in that order
* the default filenames are "./tblprior_{soft|med|hard}.txt"
* the default filenames are used iff Stable, Mtable, and Htable are set
but are not found
* WARNING: if regex is used in the filename specification, only the
first file from the list will be used. furthermore, if specified,
the table priors are applied to _all_ SSRC, MSRC, and HSRC
post [INPUT] if set, suggests the values of (Sidx,Sscl), (Midx,Mscl),
and (Hidx,Hscl) going forward, i.e., what you should set the
priors to in your next calculation for the same source -- the
suggested values are stored in the output structure
level [INPUT] percentage confidence level at which to report error
(default = 68)
details [INPUT] compute various ratios (true/false)?
(default = true)
nsim [INPUT] number of draws if algo=gibbs (default=10000)
nburnin [INPUT] number of burn-in draws if algo=gibbs
(default=5000 or NSIM/2, whichever is smaller)
outputF [INPUT] root of filename in which to place output
(default = 'none')
* output will be placed in the files OUTPUTF.txt and OUTPUTF_draws.txt
* NOTE: if OUTPUTF='none', then MC draws will be in BEHR_draws.txt
BEHRdir [INPUT] full path to directory where BEHR executable resides
(default = '/data/fubar/kashyap/AstroStat/BEHR')
verbose [INPUT] controls chatter
_extra [JUNK] here only to prevent crashing the program
requirements
uses subroutines HIPD_INTERVAL() and MODALPOINT()
BEHR executable must be installed in BEHRDIR
BEHR should be executable under the shell via SPAWN
BEHR output assumed to be compatible with 12-19-2005 version
side-effects
potentially creates numerous ascii files in $cwd or `basedir OUTPUTF`
history
vinay kashyap (Mar07; based on behr_hug.pro)
bug correction with NaNs not being caught in some cases
(VK; Mar07)
added keywords Stable,Mtable,Htable (VK; Feb08)
etymology
getting color-color diagrams using BEHR
(that's my story and I'm sticking to it)
(See pro/stat/colbehr.pro)
PROJECT: CHIANTI
PROJECT: CHIANTI
CHIANTI is an atomic database package for the calculation of
continuum and emission line spectra from astrophysical plasmas. It is a
collaborative project involving the Naval Research Laboratory
(Washington DC, USA), the Arcetri Observatory (Firenze, Italy), and the
Cambridge University (United Kingdom).
NAME:
convertname
PURPOSE:
Ion names as character strings are converted into
numerical values (note c_2 is C II or C^+1
in spectroscopic or atomic notation)
CATEGORY:
naming utility
CALLING SEQUENCE:
CONVERTNAME,Name,Iz,Ion
INPUTS:
Name: such as 'c_2'
OUTPUTS:
Iz: nuclear charge Z (6 for 'c_2', the equivalent of C II)
Ion: ionization stage: (2 for 'c_2')
OPTIONAL OUTPUTS
DIELECTRONIC Set to 1 if NAME has a 'd' appended to it
(indicating dielectronic recombination data) else
set to 0
EXAMPLE:
> convertname,'c_2',iz,ion
> print,iz,ion
> 6,2
> convertname,'o_6d',iz,ion
> print,iz,ion
> 8,6
MODIFICATION HISTORY:
Written by: Ken Dere
March 1996: Version 2.0
October 1999: Version 3. by kpd
Ver.4, 11-Dec-01, Peter Young
Revised routine, removing ch_repstr call.
Added DIELECTRONIC optional output.
(See pro/external/convertname.pro)
NAME:
CONVERT_TO_TYPE
PURPOSE:
Converts its input argument to a specified data type.
AUTHOR:
FANNING SOFTWARE CONSULTING
David Fanning, Ph.D.
1645 Sheely Drive
Fort Collins, CO 80526 USA
Phone: 970-221-0438
E-mail: davidf@dfanning.com
Coyote's Guide to IDL Programming: http://www.dfanning.com
CATEGORY:
Utilities
CALLING SEQUENCE:
result = Convert_To_Type(input, type)
INPUT_PARAMETERS:
input: The input data to be converted.
type: The data type. Accepts values as given by Size(var, /TNAME) or Size(var, /TYPE).
OUTPUT_PARAMETERS:
result: The input data is converted to specified data type.
KEYWORDS:
None.
RESTRICTIONS:
Data types STRUCT, POINTER, and OBJREF are not allowed.
MODIFICATION HISTORY:
Written by David W. Fanning, 19 February 2006.
(See pro/external/convert_to_type.pro)
procedure conv_rmf convolves input energy spectrum with the response matrix unlike the combination of RDRESP and FOLD_RESP, this one goes easy on memory at the expense of speed. syntax conv_rmf,nrg,flx,chan,spec,rmf,effar=effar,nrgar=nrgar,$ rmfstr=rmfstr,fchcol=fchcol,/shift1,verbose=verbose parameters nrg [INPUT; required] mid-bin values at which input FLX is defined. * usually [keV], depends on RMF file * if size = N(FLX)+1, assumed to be bin-boundaries flx [INPUT; required] fluxes in [ph/bin/(s/cm^2/sr/...)] chan [OUTPUT; required] bin boundary values at which output SPEC is defined * same units as NRG spec [OUTPUT; required] output, convolved spectrum rmf [INPUT; required] name of OGIP-compliant FITS file containing the response matrix * may also be a structure, such as the output of RD_OGIP_RMF() keywords effar [INPUT] effective areas [cm^2] nrgar [INPUT] energies at which EFFAR are defined * same units as NRG * sizes of EFFAR and NRGAR must match * range of NRGAR must overlap NRG * if EFFAR and NRGAR are legal, FLX is multiplied by EFFAR (over intersection of NRGAR and NRG; zeroed elsewhere) before convolving with response matrix * of course, this multiplication makes sense only if units on FLX are [ph/bin/cm^2/...] rmfstr [OUTPUT] structure containing info of response matrix, the output of RD_OGIP_RMF() verbose [INPUT] controls chatter _extra [INPUT ONLY] pass defined keywords to subroutines: RD_OGIP_RMF: FCHCOL, SHIFT1 REBINW : SLOWOK restrictions requires IDLAstro library requires RD_OGIP_RMF() requires REBINW() requires BINERSP() history vinay kashyap (Apr99) added keyword SHIFT1 (VK; JanMMI) bug: incorrect usage of WVLAR instead of NRGAR in effar block (Erica R; Aug01) deleted keyword SHIFT1, cleaned up and speeded up (VK; Nov2001) bug: when spectrum size < RMF size and RMF has multiple groups, was using only the first group (VK; Nov'02) bug: was crashing when NRG was outside the bounds of the RMF in cases where the RMF was at a higher resolution (LL/VK; Apr'03) added keyword VERBOSE (VK; MarMMV) bug: when input spectrum energy range was smaller than RMF range, was assuming RMF was at higher resolution even if it wasn't (VK; Aug'07) bug: wasn't checking for frequency beating between input spectrum grid and input RMF grid in above case; now if it does, and calls REBINW() and forces input spectrum to the right grid (VK; Aug'07) bug: EFFAR multiplication was being ignored when input spectrum had to be rebinned to match RMF grid (VK; Oct'07)
(See pro/util/conv_rmf.pro)
NAME:
CURVEFIT
PURPOSE:
Non-linear least squares fit to a function of an arbitrary
number of parameters. The function may be any non-linear
function. If available, partial derivatives can be calculated by
the user function, else this routine will estimate partial derivatives
with a forward difference approximation.
CATEGORY:
E2 - Curve and Surface Fitting.
CALLING SEQUENCE:
Result = CURVE_FIT(X, Y, Weights, A, SIGMA, FUNCTION_NAME = name, $
ITMAX=ITMAX, ITER=ITER, TOL=TOL, /NODERIVATIVE)
INPUTS:
X: A row vector of independent variables. This routine does
not manipulate or use values in X, it simply passes X
to the user-written function.
Y: A row vector containing the dependent variable.
Weights: A row vector of weights, the same length as Y.
For no weighting,
Weights(i) = 1.0.
For instrumental (Gaussian) weighting,
Weights(i)=1.0/sigma(i)^2
For statistical (Poisson) weighting,
Weights(i) = 1.0/y(i), etc.
A: A vector, with as many elements as the number of terms, that
contains the initial estimate for each parameter. IF A is double-
precision, calculations are performed in double precision,
otherwise they are performed in single precision. Fitted parameters
are returned in A.
KEYWORDS:
FUNCTION_NAME: The name of the function (actually, a procedure) to
fit. IF omitted, "FUNCT" is used. The procedure must be written as
described under RESTRICTIONS, below.
ITMAX: Maximum number of iterations. Default = 20.
ITER: The actual number of iterations which were performed
TOL: The convergence tolerance. The routine returns when the
relative decrease in chi-squared is less than TOL in an
interation. Default = 1.e-3.
CHI2: The value of chi-squared on exit (obselete)
CHISQ: The value of reduced chi-squared on exit
NODERIVATIVE: IF this keyword is set THEN the user procedure will not
be requested to provide partial derivatives. The partial
derivatives will be estimated in CURVEFIT using forward
differences. IF analytical derivatives are available they
should always be used.
OUTPUTS:
Returns a vector of calculated values.
A: A vector of parameters containing fit.
OPTIONAL OUTPUT PARAMETERS:
Sigma: A vector of standard deviations for the parameters in A.
COMMON BLOCKS:
NONE.
SIDE EFFECTS:
None.
RESTRICTIONS:
The function to be fit must be defined and called FUNCT,
unless the FUNCTION_NAME keyword is supplied. This function,
(actually written as a procedure) must accept values of
X (the independent variable), and A (the fitted function's
parameter values), and return F (the function's value at
X), and PDER (a 2D array of partial derivatives).
For an example, see FUNCT in the IDL User's Libaray.
A call to FUNCT is entered as:
FUNCT, X, A, F, PDER
where:
X = Variable passed into CURVEFIT. It is the job of the user-written
function to interpret this variable.
A = Vector of NTERMS function parameters, input.
F = Vector of NPOINT values of function, y(i) = funct(x), output.
PDER = Array, (NPOINT, NTERMS), of partial derivatives of funct.
PDER(I,J) = DErivative of function at ith point with
respect to jth parameter. Optional output parameter.
PDER should not be calculated IF the parameter is not
supplied in call. IF the /NODERIVATIVE keyword is set in the
call to CURVEFIT THEN the user routine will never need to
calculate PDER.
PROCEDURE:
Copied from "CURFIT", least squares fit to a non-linear
function, pages 237-239, Bevington, Data Reduction and Error
Analysis for the Physical Sciences. This is adapted from:
Marquardt, "An Algorithm for Least-Squares Estimation of Nonlinear
Parameters", J. Soc. Ind. Appl. Math., Vol 11, no. 2, pp. 431-441,
June, 1963.
"This method is the Gradient-expansion algorithm which
combines the best features of the gradient search with
the method of linearizing the fitting function."
Iterations are performed until the chi square changes by
only TOL or until ITMAX iterations have been performed.
The initial guess of the parameter values should be
as close to the actual values as possible or the solution
may not converge.
EXAMPLE: Fit a function of the form f(x) = a * exp(b*x) + c to
sample pairs contained in x and y.
In this example, a=a(0), b=a(1) and c=a(2).
The partials are easily computed symbolicaly:
df/da = exp(b*x), df/db = a * x * exp(b*x), and df/dc = 1.0
Here is the user-written procedure to return F(x) and
the partials, given x:
pro gfunct, x, a, f, pder ; Function + partials
bx = exp(a(1) * x)
f= a(0) * bx + a(2) ;Evaluate the function
IF N_PARAMS() ge 4 THEN $ ;Return partials?
pder= [[bx], [a(0) * x * bx], [replicate(1.0, N_ELEMENTS(f))]]
end
x=findgen(10) ;Define indep & dep variables.
y=[12.0, 11.0,10.2,9.4,8.7,8.1,7.5,6.9,6.5,6.1]
Weights=1.0/y ;Weights
a=[10.0,-0.1,2.0] ;Initial guess
yfit=curvefit(x,y,Weights,a,sigma,function_name='gfunct')
print, 'Function parameters: ', a
print, yfit
end
MODIFICATION HISTORY:
Written, DMS, RSI, September, 1982.
Does not iterate IF the first guess is good. DMS, Oct, 1990.
Added CALL_PROCEDURE to make the function's name a parameter.
(Nov 1990)
12/14/92 - modified to reflect the changes in the 1991
edition of Bevington (eq. II-27) (jiy-suggested by CreaSo)
Mark Rivers, U of Chicago, Feb. 12, 1995
- Added following keywords: ITMAX, ITER, TOL, CHI2, NODERIVATIVE
These make the routine much more generally useful.
- Removed Oct. 1990 modification so the routine does one iteration
even IF first guess is good. Required to get meaningful output
for errors.
- Added forward difference derivative calculations required for
NODERIVATIVE keyword.
- Fixed a bug: PDER was passed to user's procedure on first call,
but was not defined. Thus, user's procedure might not calculate
it, but the result was THEN used.
Steve Penton, RSI, June 1996.
- Changed SIGMAA to SIGMA to be consistant with other fitting
routines.
- Changed CHI2 to CHISQ to be consistant with other fitting
routines.
- Changed W to Weights to be consistant with other fitting
routines.
_ Updated docs regarding weighing.
Vinay Kashyap, CfA, Dec 1998
- Changed name to CURVE_FIT.
- Added keyword _EXTRA=E to calling sequence and to all
calls to CALL_PROCEDURE.
(See pro/external/curve_fit.pro)
function cut_id add or delete IDs to an existing ID structure syntax newid=cut_id(idstr,idx,delet=delet,addon=addon,verbose=verbose,$ /incieq,dbdir=dbdir,sep=sep,pres=pres,logP=logP,n_e=n_e,$ chifil=chifil,chidir=chidir,eqfile=eqfile) warning will not keep track of relative fluxes or their errors properly, though the total flux will be conserved. parameters idstr [INPUT; required] an ID structure (see LINEID.PRO) containing line identifications of spectral features idx [INPUT; required] a zero-based index of the feature being (re)edited * must be a scalar * nothing happens if IDX is outside the legal range defined by IDSTR keywords delet [INPUT] set to a zero-based index to IDSTR.(IDX+1) * if float, assumed to refer to IDSTR.(IDX+1).WVL * if string, assumed to be in same format as understood by RD_LIST: "Z ION <sep> WAVE <sep> SOURCE <sep> DESCRIPTION" (RD_LIST will in fact be called in order to decipher it) * may be an array addon [INPUT] string describing the line to be added as an ID to IDSTR.(IDX+1) * must be in same format as that understood by RD_LIST: "Z ION <sep> WAVE <sep> SOURCE <sep> DESCRIPTION" * may be an array eps [INPUT] a small number * default is 1e-5 verbose [INPUT] controls chatter _extra [INPUT ONLY] pass defined keywords to RD_LIST: INCIEQ,DBDIR,SEP,PRES,LOGP,N_E,CHIFIL,CHIDIR,EQFILE LINEFLX: DEM,ABUND,NOPH,EFFAR,WVLAR restrictions requires subroutines: RD_LIST, RD_LINE, FOLD_IONEQ, READ_IONEQ, RD_IONEQ, LINEFLX, GETABUND, WHEE, SYMB2ZION, LAT2ARAB, CAT_LN requires IDL 5.3+ (because of use of STRSPLIT, STRJOIN, STRCMP) history vinay kashyap (JanMMI)
(See pro/cut_id.pro)
Purpose: checks input data for type, ndimension, etc
(uses IDL size function results)
Keyword Parameters:
type - if set, return idl data type (0,1,2..8) from size function
ndimen - if set, return number dimensions (size(0))
nimages - if set, return #images (1->2D, NI->3D, else 0)
xsize - if set, return X size (size(data))(1)
nx - synonym for xsize
ysize - if set, return Y size (size(data))(2)
ny - synonym for ysize
orr - if set, return value is OR of all boolean flags (def=AND)
string/struct/undefined - if set, return true if type matches
scalar, vector - if set, true if data of specified variety
Calling Examples:
if (data_chk(p1,/type) eq data_chk(p2,/type)) then...
case data_chk(data,/type) of...
if data_chk(data,/string,/scalar) then ...
if data_chk(data,/string,/struct,/undef,/orr)
case data_chk(maybe_cube,/nimages) of...
History:
27-Apr-1993 (SLF)
21-Mar-1994 (SLF) documentation header
10-oct-1996 (SLF) add SCALAR (synonym for the historical mispell SCALER)
2-oct-1997 (SLF) add NIMAGES
15-nov-1997 (SLF) add XSIZE and YSIZE keyword and function
Restrictions:
some keywords are mutually exclusive - for self-documenting code
and reduction of code duplicataion
(See pro/external/data_chk.pro)
function demacs a DEM editor that allows interactively "tweaking" a specified DEM syntax dem=demacs(logt,dem0=dem0,logt0=logt0,norm=norm,slope=slope,$ group=group,igroup=igroup, PLOT_KEYWORDS) parameters logt [INPUT; required] log10(Temperature [K]) at which DEM is to be determined keywords dem0 [INPUT] initial DEM (just for guidance). * if size does not match LOGT0, gets stretched to suit. * overrides SLOPE, but not NORM * if max(DEM0)<100, assumed to be in log form -- i.e., plots are in "linear" form and +vity is not enforced. logt0 [INPUT] logT at which DEM0 are defined. if not given, assumed to be LOGT. norm [INPUT] if set, first element of initial DEM is set to NORM slope [INPUT; default=1] if given, generates initial DEM=NORM*T0^(SLOPE) * if DEM0 is specified, SLOPE is ignored group [OUTPUT] long-int array of same size as LOGT containing grouping information igroup [OUTPUT] array of same size as GROUP, containing index of "representative element" of each group _extra [INPUT] allows passing defined keywords (e.g., YRANGE) to PLOT. restrictions requires X-windows, display, and mouse capability. history vinay kashyap (Feb97) corrected log behavior, added _extra to PLOT (VK; Apr97) cleaned up log behavior, added keystrokes *,/,v,^,z (VK; Aug98)
(See pro/demacs.pro)
EXPLANATION
This routine descales all types of spline fits into upsilons or
rates, i.e., it does both electron upsilons and proton rates,
and both 5-point and 10-point splines. In addition it can
simultaneously descale several temperatures at once.
INPUTS
TEMP Temperature(s), K.
SPLSTR Structure output by read_splups.
INDEX Index of structure.
OUTPUTS
UPS Upsilon value(s) at temperature(s) TEMP.
EXAMPLES
read_splups,splupsfile,splstr
descale_all,[1.e6,2.e6],splstr,5,ups
print,ups
HISTORY
Ver.1, 15-Mar-01, Peter Young
adapted from Ken Dere's descale_ups.pro.
Ver.2, 12-Nov-01, Peter Young
added type 6 transitions (for protons)
(See pro/external/descale_all.pro)
PROJECT: CHIANTI
CHIANTI is an atomic database package for the calculation of
continuum and emission line spectra from astrophysical plasmas. It is a
collaborative project involving the Naval Research Laboratory
(Washington DC, USA), the Arcetri Observatory (Firenze, Italy), and the
Cambridge University (United Kingdom).
NAME:
DESCALE_UPS_VK
(used to be) DESCALE_UPS
PURPOSE:
convert from Burgess-Tully scaling spline fits to Upsilons
CATEGORY:
science.
CALLING SEQUENCE:
DESCALE_UPS,Index,Jndex,xt,upsilion, T_TYPE,C_UPS,SPLUPS
INPUTS:
Index: index of lower energy level (lowest level is 1)
Jndex: index of upper energy level (lowest level is 1)
xt: scaled temperature
T_TYPE: 2 dimensional array contain values of the transition type
C_UPS: 2 dimensional array containing values of the Burgess and Tully
scaling parameter c
SPLUPS: spline fits to the scaled Upsilons
OPTIONAL INPUTS:
None:
KEYWORD PARAMETERS:
None:
OUTPUTS:
Upsilon: the Maxwellian averaged collision strength
COMMON BLOCKS:
NONE:
;common elvlc,l1a,term,conf,ss,ll,jj,ecm,eryd,ecmth,erydth,eref
;common wgfa, wvl,gf,a_value
;common upsilon,t_type,c_ups,splups
PROCEDURE:
see Burgess and Tully, 1992, Astron and Astrophys, 254, 436.
EXAMPLE:
;
MODIFICATION HISTORY:
Written by: Ken Dere
March 1996: Version 2.0
December 1998: Include transition type 5 (kpd)
September 2000: Add parameters T_TYPE,C_UPS,SPLUPS and comment out
all common blocks (Vinay Kashyap)
THIS ROUTINE IS OBSOLETE FOR CHIANTIv4+ (VK; Jun02)
(See pro/external/descale_ups_vk.pro)
function detect_limit compute and return the counts upper limit for detection of a source at a given significance, given the background counts. description compute the cumulative significance of obtaining a specified number of counts given the background, and assume that a source would be considered detected if the counts were to exceed the NSIG threshold. see Pease, Drake, & Kashyap (2006, ApJ 636, 426) for full description. briefly, computes the probability that as many as D counts can be observed for a given background b, p(=<D|b), and thence the probability that a given number of counts can be obtained in an observation simply due to the background. The number of counts required for a detection at a specified probability is the upper limit. Note also that this goes only halfway towards a full upper limit (see Kashyap et al., 2008, AAS-HEAD, 2008.03.03) in that this produces an upper limit in counts space and not in intrinsic flux space. syntax ul=detect_limit(bkg,nsig,asrc=asrc,abkg=abkg,bgalt=bgalt,abgalt=abgalt,$ nsim=nsim,ulsig=ulsig,ulsim=ulsim,/gaussy,nxbin=nxbin,verbose=verbose) parameters bkg [INPUT; required] counts in the background nsig [INPUT] the Gaussian-equivalent sigma at which to compute upper limit to detection * if not given, assumed to be 3 keywords asrc [INPUT; default=1] area in which source counts are collected abkg [INPUT; default=1] area in which background counts are collected bgalt [INPUT; default=0] alternative set of background contamination, say from a different area of the instrument, or from a model, or from an extended source, etc. * may be an array * NOTE: the _same_ BGALT is appended to _all_ elements of BKG abgalt [INPUT; default=asrc] area in which BGALT is "collected" * size must match BGALT * if size < size(BGALT), first element is assumed to be the default nsim [INPUT; default=0] number of Monte Carlo simulations to run to account for error in background * setting this results in computing the upper limit for a number of realizations of BKG; the resulting 1-sigma range in the value of computed upper limits is reported in ULSIG, and a conservative upper limit based on combining all of the simulations is returned in ULSIM[*,0] ulsig [OUTPUT] the 1-sigma error on the upper limit, estimated by bootstrapping BKG ulsim [OUTPUT] a 2D array of size (NBKG,NSIM+1) which contains all the simulated limits * ULSIM[*,0] is identical to the primary output if NSIM=0 * for NSIM>0, ULSIM[*,0] is the conservative limit that is derived from the coadded probability distributions that take into account the variations in the background gaussy [INPUT] if set, computes the limit corresponding to the significance matching the location of the NSIG-sigma *intercept* of a Gaussian, rather than matching the total area under the curve. nxbin [INPUT] number of bins to use in the integration * the integration is carrid out over a range of 0..5*E(bg) or 20, whichever is greater. by default, the number of bins is set by the step size, which is set to 1 count, -- unless E(bg) < 1, in which case a bin width of 0.05 is used by default * changing NXBIN does not change the range, but only changes the bin size. * a hard lower limit of 20 is set -- cannot use a bin width larger than 1 count verbose [INPUT] controls chatter subroutines LNPOISSON() KILROY history vinay kashyap (Apr2004) added keyword ULSIG; changed name from PUPLIM to DETECT_LIMIT (VK; May2004) added keyword NXBIN (VK; Sep2004) added keywords BGALT and ABGALT (VK; Dec2004) modified output behavior of ULSIM[*,0]; now BGALT can be 0 (VK; Mar2005)
(See pro/stat/detect_limit.pro)
procedure did2emis2em for a specified line, read in emissivities at various densities, compute fluxes, and invert to get emission measure estimates at each density. syntax did2emis2em,idstr,edens,fluxes=fluxes,/deblend,DEM=DEM,logT=logT,$ ldir=ldir,NH=NH,EMs=EMs,FXs=FXs,xtitle=xtit,ytitle=ytit,title=tit,$ verbose=v,dWVL=dWVL,eps=eps,/incieq,mapping=mapping,chifil=chifil,$ chidir=chidir,eqfile=eqfile,abund=abund,/noph,effar=effar,wvlar=wvlar,$ /kev,tol=tol,fH2=fH2,He1=He1,HeII=HeII,/Fano,/wam,/bam,/mam,/Zbeda,$ /Ibeda,Wbeda=Wbeda,/Lbeku,wform=wform,wstyle=wstyle,ziform=ziform,$ /nuthin,/noday,/notime,/nouser,/nopack,stacol=stacol,stasiz=stasiz,$ stathk=stathk,charsize=charsize,align=align parameters idstr [INPUT; required] structure containing ID information (see LINEID for description) edens [I/O] array of electron densities [cm^-3] * default is [1e8,1e14] * if not specified, insufficiently specified, or otherwise meaningless, uses default and overwrites input keywords fluxes [I/O] observed fluxes * if size does not match the number of components or the number of IDs, will instead use IDSTR.(#).FLUX and overwrite the input * see UPDATID() for details deblend [INPUT] if _not_ set (which is the default), squishes multiple IDs of a single feature into one single flux and calculates only a "composite" ID. on the other hand, if set, leaves the multiple IDs alone. DEM [INPUT] differential emission measure to use in computing the fluxes * default is to use 1e12 [cm^-5] at each specified LOGT logT [INPUT] temperatures at which DEM is defined * if not specified, set to 6.5, unless DEM has more elements, in which case is interpolated into the 4..8 range ldir [INPUT] string array of database directories to use to search for the emissivities * default is '$CHIANTI' * if size does not match either number of IDSTR components or IDs, only the first element is used NH [INPUT] H column density [cm^-2] EMs [OUTPUT] array of emission measures EMS[EDEN,WVL] that are plotted. * at any given point, EMS = \int DEM[LOGT]*dLOGT, unless DEM has one or two elements only, in which case EMS = SUM(DEM) * what this implies is that the input DEM is scaled at each density to match the observed flux at each wavelength and this scaled DEM is represented in the plot by the summed measure. This ensures that all wavelengths considered are compared over the same temperature range, and the ambiguity regarding the temperature of peak emissivity, etc. has been removed. FXs [OUTPUT] the predicted fluxes FXS[EDEN,WVL] which are used to scale the EMs above. xtitle [INPUT] passed to PLOT ytitle [INPUT] passed to PLOT title [INPUT] passed to PLOT verbose [INPUT] controls chatter _extra [INPUT] pass defined keyword values to subroutines: ID2EMIS2ID: DWVL RD_LIST: EPS, INCIEQ, MAPPING FOLD_IONEQ: CHIFIL, EQFILE RD_IONEQ: CHIDIR SQUISHEM: ABUND LINEFLX: ABUND, NOPH, EFFAR, WVLAR, KEV ARRAYEQ: TOL ISMTAU: fH2, He1, HeII, FANO, WAM, BAM, MAM IDLABEL: ZBEDA, IBEDA, WBEDA, LBEKU, WFORM, WSTYLE, ZIFORM STAMPLE: NUTHIN,NODAY,NOTIME,NOUSER,NOPACK,STACOL,STASIZ,STATHK XYOUTS: CHARSIZE,ALIGN n_e [IGNORE] here only to trap boxing gloves restrictions requires IDL5.3+ because of use of STRMATCH, etc. requires subroutines BAMABS CAT_LN FOLD_IONEQ GETABUND ID2EMIS2ID IDLABEL INICON ISMTAU LAT2ARAB LINEFLX RDABUND RD_IONEQ RD_LINE RD_LIST READ_IONEQ SQUISHEM SYMB2ZION SYZE UPDATID WHEE ZION2SYMB side-effects generates a plot which erases any plot that exists beforehand history vinay kashyap (JanMMI; based on DID2EM) improved color-scale setting for 24-bit consoles (VK; FebMMI)
(See pro/did2emis2em.pro)
function dorren
computes spot modulated stellar light curve using Dorren(1987)
syntax
curve = x, A, f, pder,ldstr=ldstr,ldspt=ldspt,rdns=rdns,$
photT=photT,spotT=spotT,wvl=wvl
parameters
x [INPUT;required] absissca points at which to compute spot
model. These are angular displacements with respect to the
initial spot longitude specified in A, in units of degrees.
A [INPUT;required] spot parameters for dorrens model in
radians where:
a(0) = stellar inclination to LOS
a(1) = spot latitude
a(2) = spot size
a(3) = initial spot longitude
pder[OUTPUT] partial derivative with respect to each
parameter in A
keywords
ldstr [INPUT] limb-darkening coefficient for star
default = 0.32 Van Hamme(1994)
ldspt [INPUT] limb-darkening coefficient for spot (umbra)
rdns [INPUT] set this parameter if spot parameters and absissca are input
in radians rather than degrees
photT [INPUT] photospheric temperature in degrees Kelvin
spotT [INPUT] spot temperature in degrees Kelvin
wvl [INPUT] wavelength in nanometers
subroutines
MOD
history
liwei lin (Sep 03) better documentation/functionality than
mod.pro. switch to modular
(Jul 06) removed loop over longitudes and add pder
(Aug 06) BUGFIX: note that T should be defined as:
Pi - ATAN((-SIN(h(tmp)))*(TAN(b1(tmp)))) when B1>=Pi/2
ATAN((SIN(h(tmp)))*(TAN(b1(tmp)))) when B1<Pi/2
to avoid dicontinuities. This differs from Dorren (1987)
formalism only in the '=' sign placement
(See pro/timing/dorren.pro)
function dove_ciclo this function is acts as a wrapper to all the morphological operations that are used to extract pixels that define a loop from an image of the solar corona the general workflow is: - [CONVOL] - enhance contrast by subtracting background - [MORFO_ROTTANGOLI] open with rectangular structure element - [MORPH_CLOSE] close with small square kernel - [CONVOL] enhance contrast again by subtracting background - [MORFO_SOGLIA] threshold to selectively enhance structures - [MORPH_CLOSE] close the image again - [MORFO_SEGMENTO] group contiguous pixels into region blobs and throw out regions deemed too small - [ROI_SELEZIONI] interactively select some blobs for further study - [MORFO_SCHELETRO] create a skeleton for these blobs - [MORFO_POTARE] prune the skeleton to make a cleaner image - [] extract the pixels that form the pruned skeleton and convert them to heliographic coordinates keywords control whether any of the calls to the subroutines should be skipped or repeated. keywords also set up the input parameters to the subroutines. syntax oimg=dove_ciclo(inpimg,$ /xbgdat,icell=icell,$ /xrect,xsize=xsize,ysize=ysize,angle=angle,gray=gray,/readd,/reuse,$ /xrclose,$ /xbgsub,jcell=jcell,$ /xthresh,thrthr=thrthr,nsigma=nsigma,mimg=mimg,gamma=gamma,$ centrale=centrale,bitlev=bitlev,/zeromin,$ /xtclose,$ /xblob,blbthr=blbthr,areas=areas,$ subidx=subidx,/hardcut,deepimg=deepimg,$ /xroi,pixroi=pixroi,roithr=roithr,jitter=jitter,/hadrian,$ /xskel,skthr=skthr,rmin=rmin,$ /xprune,minpix=minpix,$ outx=outx,outy=outy,sunrad=sunrad,pixsiz=pixsiz,$ offsetx=offsetx,offsety=offsety,suncenx=suncenx,sunceny=sunceny,$ loopsav=loopsav, verbose=verbose) parameters inpimg [INPUT; required] image of the corona keywords For keywords with names beginning with "X", when set, the corresponding call to the subroutine is skipped. But if that keyword is set to a negative number -N, the subroutine is called N times in sequence (unless otherwise noted below). xbgdat [INPUT] if set, skips initial background subtraction icell [INPUT; default=1] defines the size of the cell to subtract background - the "source" is measured in a square of size 2*ICELL+1 and the "background" in a surrounding layer of 1pix xrect [INPUT] if set, skips morphological open with rectangular structure element xsize [INPUT; default=1] width of the structure element ysize [INPUT; default=10] height of the structure element angle [INPUT; default=findgen(36)*5] angle of tilt of the rectangle defined by (XSIZE,YSIZE) in [degrees] xrclose [INPUT] if set, skips morphological close apres open xbgsub [INPUT] if set, skips background subtraction of modified image jcell [INPUT; default=ICELL] defines size of cell to subtract background xthresh [INPUT] if set, skips image thresholding thrthr [OUTPUT] value used if histogram-threshold is applied nsigma [INPUT] if given, sets the threshold at mean+NSIGMA*stddev * default is 1 * if mean+NSIGMA*stddev < 0, this is not applied mimg [OUTPUT] the median threshold image, constructed if CENTRALE is set xtclose [INPUT] if set, skips morphological close apres thresholding xblob [INPUT] if set, skips grouping the image into distinct regions * if set to -ve number, gets automatically unset areas [OUTPUT] number of pixels in each labeled region arclev [INPUT; default=0.68] fraction of the area of the largest blob as the threshold below which to discard blobs deepimg [OUTPUT] an image that depicts the depth at which a particular pixel was added to a blob blbthr [INPUT; default=0] threshold value at which to convert grouped image to bitmap xroi [INPUT] if set, skips selecting a subset of regions * if set to -ve number, gets automatically unset pixroi [INPUT] if set, assumed to be a set of pixels around which to pick out the region of interest non-interactively roithr [INPUT; default=0] threshold value at which to ignore pixels xskel [INPUT] if set, skips making skeletons of regions * if set to -ve number, gets automatically unset rmin [INPUT; default=1] radius of circle that acts as structure element skthr [INPUT; default=0] threshold to filter region xprune [INPUT] if set, skips pruning skeletons * if set to -ve number, gets automatically unset minpix [INPUT; default=4] number of pixels above which a given branch must be kept while pruning outx [OUTPUT] the x-pixel indices of the pruned points outy [OUTPUT] the y-pixel indices of the pruned points * if SUNRAD is set, then floating point values defining the locations in heliospheric coordinates are returned sunrad [INPUT] radius of the Sun as seen from the spacecraft, in [arcsec] * if this is set, then OUTX and OUTY are transformed to heliospheric coordinate system * if set to 1, assumed to be (RSun/AU)*(180/!pi)*3600 pixsiz [INPUT; default=0.5] pixel size, in [arcsec] offsetx [INPUT; default=0] offset to be applied to OUTX to bring them in line with SUNCENX offsety [INPUT; default=0] offset to be applied to OUTY to bring them in line with SUNCENY suncenx [INPUT; default=0] the X coordinate distance of the center of INPIMG from Sun center in [arcsec] sunceny [INPUT; default=0] the Y coordinate distance of the center of INPIMG from Sun center in [arcsec] verbose [INPUT] controls chatter loopsav [INPUT] if set to a filename, saves all the intermediate arrays in the named IDL save file _extra [INPUT ONLY] pass defined keywords to subroutines MORFO_ROTTANGOLI: GRAY, READD, REUSE MORFO_SOGLIA: GAMMA, CENTRALE, BITLEV, ZEROMIN MORFO_SEGMENTO: SUBIDX, HARDCUT ROI_SELEZIONI: JITTER, HADRIAN subroutines MORFO_ROTTANGOLI MORFO_SOGLIA [GMASCL,ERROR_MESSAGE,PATH_SEP,SCALE_VECTOR,CONVERT_TO_TYPE,FPUFIX] MORFO_SEGMENTO ROI_SELEZIONI MORFO_SCHELETRO MORFO_POTARE INICON KILROY history vinay kashyap (May2007)
(See pro/ssw/dove_ciclo.pro)
procedure drakopy
set up the environment to make plots per individual preferences
syntax
drakopy,milieu,/HELP,/UNDO,MAGNIFY=MAGNIFY,VERBOSE=VERBOSE
parameters
milieu [INPUT] string specifying what sort of environment to set to.
accepted values are:
-- 'jeremy'
-- 'olivia'
-- 'jjd' (jeremy old preferences)
-- 'drake' (jeremy older preferences)
-- 'vinay' (different depending device at call time)
-- 'poster'
-- 'original'
keywords
help [INPUT] if set, prints out a somewhat verbose help
undo [INPUT] if set, restores original state
verbose [INPUT] controls chatter level
magnify [INPUT] scale all the variables by this factor
_extra [JUNK] here only to prevent crashing the program
commons
drakopy {old_P, old_X, old_Y, old_Z}
description
sets the !P, !X, !Y, and !Z variables to special and specific values.
stores the old values in a common block for later access.
example
set_plot,'ps' & device,file='/tmp/test_drakopy.ps'
plot,findgen(10),title='MILIEU=Original'
drakopy,'jeremy' & plot,findgen(10),title='MILIEU=Jeremy'
drakopy,'olivia' & plot,findgen(10),title='MILIEU=Olivia'
drakopy,'original' & plot,findgen(10),title='MILIEU=Original'
drakopy,/undo & plot,findgen(10),title='MILIEU=previous (Olivia)'
device,/close & set_plot,'x'
history
vinay kashyap (SepMM)
added option "drake" (VK; Aug01)
added option "jjd"; added keyword MAGNIFY (VK; Feb02)
added option "vinay" (VK; May03)
added option "poster" (VK; Jul03)
told the code that "poster" is a valid option (VK; Sep03)
(See pro/util/drakopy.pro)
function drat return the ratio of fluxes in one line relative to another line for a variety of densities syntax frat=drat(line1,line2,edens,DEM,logT,fx1=fx1,fx2=fx2,$ mapping=mapping,chifil=chifil,chidir=chidir,eqfile=eqfile,$ abund=abund,/noph,effar=effar,wvlar=wvlar,NH=NH,fH2=fH2,$ He1=He1,HeII=HeII,/Fano,dbdir=dbdir,sep=sep,prefix=prefix) parameters line1 [INPUT; required] string describing the line/wavelength-range that defines the numerator line2 [INPUT; required] string describing the line/wavelength-range that defines the denominator * LINE1 and LINE2 must be in the format used by RD_LIST, i.e., Z ION <sep> WAVE <sep> DBDIR <sep> DESCRIPTION and WAVE = "WVL"/"WVL +- dW"/"WMIN,WMAX"/"WMIN-WMAX" edens [INPUT; required] electron densities at which to compute the flux ratios DEM [INPUT; required] differential emission measure to use in computing the fluxes logT [INPUT; required] temperatures at which DEM is defined * size must match that of DEM keywords fx1 [OUTPUT] the fluxes due to line 1 at each EDENS fx2 [OUTPUT] the fluxes due to line 2 at each EDENS _extra [INPUT ONLY] pass defined keywords to subroutines RD_LIST: MAPPING,DBDIR,SEP,PREFIX FOLD_IONEQ: CHIFIL,EQFILE RD_IONEQ: CHIDIR LINEFLX: ABUND,NOPH,EFFAR,WVLAR ISMTAU: NH,fH2,He1,HeII,FANO,BAM BAMABS: ABUND subroutines RD_LIST RD_LINE FOLD_IONEQ LINEFLX ISMTAU BAMABS CAT_LN history vinay kashyap (MarMM)
(See pro/drat.pro)
function dummyid
Convert a RD_LINE() style emissivity structure into a dummy
ID structure such as that generated by LINEID()
syntax
idstr=dummyid(linstr,DEM=DEM,/temp,abund=abund,/noph,$
nhne=nhne,effar=effar,wvlar=wvlar,/ikeV,/noabund)
parameters
linstr [INPUT; required] RD_LINE() or RD_LIST() style line
emissivity structure
keywords
_extra [INPUT ONLY] use this to pass defined keywords to subroutines
LINEFLX: DEM,ABUND,EFFAR,WVLAR,NOABUND,NHNE,NOPH,TEMP,IKEV
subroutines
LINEFLX()
history
liwei lin (Feb06) lift parameter check from cat_ln()
cleaned up a bit (VK; Feb06)
(See pro/dummyid.pro)
function eqt_interval computes and returns the double-sided equal-tail interval [lower_bound,upper_bound] at a specified confidence level. interval includes the central mass of the probability distribution unless explicitly required otherwise. syntax hpd=eqt_interval(f,x,/fsample,clev=clev,pdfnorm=pdfnorm,$ xaround=xaround,verbose=verbose) parameters f [INPUT; required] the array for which the confidence interval must be computed * assumed to be a density function unless FSAMPLE is set x [INPUT; optional] abscissae values * if not given, and F is a density function, then taken to be the array indices * ignored if FSAMPLE is set * if N(X).GT.1 but N(X).NE.N(F), X is ignored and FSAMPLE is set automatically keywords fsample [INPUT] if set, assumes that F is a set of samples from a density function, as opposed to being the density function itself clev [INPUT] confidence level at which to compute the intervals * default is 0.68 * if < 0, abs(CLEV) is used * if > 1 and < 100, then assumed to be given as a percentage * if > 100, then 1-1/CLEV is used pdfnorm [INPUT] if set, forces F to integrate to abs(PDFNORM) * if explicitly set to 0, does not normalize F at all * if not set, normalizes to 1 * ignored if FSAMPLE is set * WARNING: do not use this keyword unless you know what you are doing xaround [INPUT] by default, the central part of the distribution is used, i.e., an equal area is left out on both ends. if this keyword is defined, a fraction CLEV/2 of the area is used on either side of XAROUND. * if F is a sample, XAROUND must be in the same units as F * you can return single-sided intervals by setting XAROUND to min(X) or max(X) (min/max of F if sample) * if not set, XAROUND is set to the median verbose [INPUT] controls chatter _extra [INPUT ONLY] pass defined keywords to subroutines example for i=1,20 do print,eqt_interval(randomn(seed,10000L)*i,/fsample) history vinay kashyap (Apr2006)
(See pro/stat/eqt_interval.pro)
procedure erors compute asymmetric error bars on fit parameters by projecting the chi-sq surface onto each parameter axis. what this means is that for each non-frozen parameter, step through a range of values of the parameter, finding the best-fit chi-square calculated by fitting the rest of the parameters, and report that range where the chi-square increases by a set amount from the minimum. this program tries to be a wee bit clever by starting from the best-fit value and then stumbling about on either side until the chi-sq goes above the mark; the length of the strides change according to the projected location of said mark. syntax erors,x,y,a,erru,errl,ysig=ysig,freeze=freeze,dchi=dchi,$ dchi=dchi,algo=algo,maxstep=maxstep,verbose=verbose,$ erra=erra,yfunc=yfunc,$ itmax=itmax,chithr=chithr,/dumb, jumpup=jumpup,jumpdn=jumpdn,$ svdthr=svdthr,funcs=funcs, function_name=function_name,$ type=type, missing=missing,/fwhm,/norm,betap=betap, /poisson parameters x [INPUT; required] data points y [INPUT; required] Y(X) * sizes of X and Y must match a [I/O; required] parameters for user-supplied function * on input, these are assumed to be initial guesses * on output, these contain the best-fit values erru [OUTPUT; required] upper limits of confidence range interval errl [OUTPUT] lower limits of confidence range interval * if ERRL is not specified on input, ERRU will contain the average value of the upper and lower >>deviations<<. keywords ysig [INPUT] standard deviations on Y * default=sqrt(abs(Y)+0.75)+1 * if single element, then sig(Y[*])=YSIG(0) * if -ve, taken to be the fractional error freeze [INPUT] freeze numbered parameters (index starts from 0!) dchi [INPUT] how big a change in chi-square to look for? * default is 2.7 (corresponding to 90% CL) algo [INPUT] fitting algorithm * only the following are implemented: -- LevMarq+SVD (default; calls FIT_LEVMAR) -- IDL-Curvefit (calls CURVE_FIT) maxstep [INPUT] maximum number of steps to take before giving up on actually finding the bounds * default is 100 verbose [INPUT] verbosity level erra [OUTPUT] formal "curvature" errors on the best-fit parameters yfunc [OUTPUT] best-fit Y(X;A) _extra [INPUT] pass defined variables to subroutines:- FIT_LEVMAR: ITMAX, CHITHR, DUMB ADJUSTIE: TIES, VNAME LEVMARQ: JUMPUP, JUMPDN, SVDTHR LMCOEFF: FUNCS, POISSON CURVE_FIT: FUNCTION_NAME note:- FUNCS and FUNCTION_NAME refer to name of user-defined function that takes as input X and A, and returns YMODEL(X;A), and the partial derivatives of the model function wrt A. Any function that was written to work with CURVEFIT or GHRS' WFIT will do. The default for FIT_LEVMAR is X3MODEL. MK_3MODEL: TYPE MK_GAUSS: MISSING, FWHM, NORM MK_LORENTZ: BETAP, MISSING, NORM history vinay kashyap (MM.I) (yes, I _do_ know how to spell "error") added call to ADJUSTIE to handle constraints on ERRU,ERRL (VK; FebMM) what if adjusted ERRA becomes 0? (VK; MarMM) allowed halting with either "q" or "x" (VK; SepMM) force extra iteration after "r"; increased wait time (VK; JanMMI) added confirmation check for too many thawed params (VK; FebMMI) interpol was crashing because it was getting only 1 element arrays; now correctly updates all parameters if better fit is found (VK; Aug01) made various bug fixes that was causing program to go bonkers for some special cases, such as small fluxes (VK; Apr02) added hooks into MPFIT (Liwei Lin/VK; Oct02) bug correction re VERBOSE (LL; Apr03) bug correction, first step was failing sometimes when input A were integers (VK; Mar08)
(See pro/fitting/erors.pro)
NAME:
ERROR_MESSAGE
PURPOSE:
The purpose of this function is to have a device-independent
error messaging function. The error message is reported
to the user by using DIALOG_MESSAGE if widgets are
supported and MESSAGE otherwise.
In general, the ERROR_MESSAGE function is not called directly.
Rather, it is used in a CATCH error handler. Errors are thrown
to ERROR_MESSAGE with the MESSAGE command. A typical CATCH error
handler is shown below.
Catch, theError
IF theError NE 0 THEN BEGIN
Catch, /Cancel
ok = Error_Message()
RETURN
ENDIF
Error messages would get into the ERROR_MESSAGE function by
throwing an error with the MESSAGE command, like this:
IF test NE 1 THEN Message, 'The test failed.'
AUTHOR:
FANNING SOFTWARE CONSULTING
David Fanning, Ph.D.
1645 Sheely Drive
Fort Collins, CO 80526 USA
Phone: 970-221-0438
E-mail: davidf@dfanning.com
Coyote's Guide to IDL Programming: http://www.dfanning.com/
CATEGORY:
Utility.
CALLING SEQUENCE:
ok = Error_Message(the_Error_Message)
INPUTS:
the_Error_Message: This is a string argument containing the error
message you want reported. If undefined, this variable is set
to the string in the !Error_State.Msg system variable.
KEYWORDS:
ERROR: Set this keyword to cause Dialog_Message to use the ERROR
reporting dialog. Note that a bug in IDL causes the ERROR dialog
to be used whether this keyword is set to 0 or 1!
INFORMATIONAL: Set this keyword to cause Dialog_Message to use the
INFORMATION dialog instead of the WARNING dialog. Note that a bug
in IDL causes the ERROR dialog to be used if this keyword is set to 0!
TITLE: Set this keyword to the title of the DIALOG_MESSAGE window. By
default the keyword is set to 'System Error' unless !ERROR_STATE.NAME
equals "IDL_M_USER_ERR", in which case it is set to "Trapped Error'.
TRACEBACK: Setting this keyword results in an error traceback
being printed to standard output with the PRINT command. Set to
1 (ON) by default. Use TRACEBACK=0 to turn this functionality off.
OUTPUTS:
Currently the only output from the function is the string "OK".
RESTRICTIONS:
The WARNING Dialog_Message dialog is used by default.
EXAMPLE:
To handle an undefined variable error:
IF N_Elements(variable) EQ 0 THEN $
ok = Error_Message('Variable is undefined')
MODIFICATION HISTORY:
Written by: David W. Fanning, 27 April 1999.
Added the calling routine's name in the message and NoName keyword. 31 Jan 2000. DWF.
Added _Extra keyword. 10 February 2000. DWF.
Forgot to add _Extra everywhere. Fixed for MAIN errors. 8 AUG 2000. DWF.
Adding call routine's name to Traceback Report. 8 AUG 2000. DWF.
Added ERROR, INFORMATIONAL, and TITLE keywords. 19 SEP 2002. DWF.
Removed the requirement that you use the NONAME keyword with the MESSAGE
command when generating user-trapped errors. 19 SEP 2002. DWF.
Added distinctions between trapped errors (errors generated with the
MESSAGE command) and IDL system errors. Note that if you call ERROR_MESSAGE
directly, then the state of the !ERROR_STATE.NAME variable is set
to the *last* error generated. It is better to access ERROR_MESSAGE
indirectly in a Catch error handler from the MESSAGE command. 19 SEP 2002. DWF.
Change on 19 SEP 2002 to eliminate NONAME requirement did not apply to object methods.
Fixed program to also handle messages from object methods. 30 JULY 2003. DWF.
Removed obsolete STR_SEP and replaced with STRSPLIT. 27 Oct 2004. DWF.
Made a traceback the default case without setting TRACEBACK keyword. 19 Nov 2004. DWF.
(See pro/external/error_message.pro)
procedure exaped EX APED AD IDL extract line and/or continuum emissivity data from an APED file syntax exaped,filroot,emis,tlog,wvl,edens,Z,ion,desig,econf,jon,src,$ /llist,/coco,/tlist,/noerg,/okeV,atomdb=atomdb,verbose=verbose parameters filroot [INPUT; required] APED fits file root * if this includes a ".fits", the suffix is assumed to be already included and nothing extra is added. * the suffix '_linelist.fits' is automatically applied, so don't specify that part. (this is inflexible at this point because we only support reading in this file. sometime later on the other APED format will also be supported.) * looks for ".gz" automatically if the filename specified without it does not exist emis [OUTPUT; required] emissivity array, of the format (NTLOG,NWVL,NEDENS) * APED units are [ph cm^3/s], which will be converted into [1e-23 ergs cm^3/s] unless the NOERG keyword is set tlog [OUTPUT; required] log10(Temperature [K]) grid on which EMIS are placed wvl [OUTPUT; required] wavelength array * units are in [Ang] unless the OKEV keyword is set, in which case output in [keV] edens [OUTPUT; required] electron densities at which EMIS have been calculated Z [OUTPUT; required] atomic numbers of species producing the given line ion [OUTPUT; required] ionic state of species corresponding to the line desig [OUTPUT] (2,NWVL) string array of lower and upper level designation * currently, just the lower and upper levels econf [OUTPUT] (2,NWVL) string array of electron configuration of lower and upper levels of the transition in question * not implemented yet jon [OUTPUT] the ionic state that matters to the ion balance * not implemented yet src [OUTPUT] numeric code specifiying that this is an APED output, set to 3 keywords llist [INPUT] if set, reads from FILROOT_linelist.fits * currently this is the only option supported coco [INPUT] if set, reads from FILROOT_coco.fits * NOT YET IMPLEMENTED tlist [INPUT] if set, reads from FILROOT_line.fits * NOT YET IMPLEMENTED noerg [INPUT] if set, EMIS will be output in same units as in the APED files, i.e., [ph cm^3/s] okeV [INPUT] if set, WVL will be converted from [Ang] to [keV] atomdb [INPUT] directory containing input file * default is !ATOMDB * hardcoded default is /data/atomdb/ * prepended to FILROOT only if FILROOT does not begin with a '/' (UNIX), ':' (MacOS), '\' (Windows), or '[' (VMS?) logT [INPUT] if defined as an array, then interpolates/extrapolates EMIS(NTLOG,NWVL) to a new grid EMIS(NLOGT,NWVL) prior to output. TLOG will be overwritten by LOGT. verbose [INPUT] controls chatter _extra [JUNK] ignore -- here only to prevent crashing the program subroutines SETSYSVAL ARRAYEQ REBINX description The ATOMDB file format is described in http://cxc.harvard.edu/atomdb/features_formats.html In this case, the line list (grouped by line) file is read, and in particular, only the block EMISSIVITY is read. history vinay kashyap (Jul02) bug correction when multiple densities are included (VK; Mar03)
(See pro/mkemis/exaped.pro)
Project : SDAC
Name : EXIST
Purpose : To See if variable Exists
Explanation : So obvious, that explaining it will take more
lines than the code.
Use : A=EXIST(VAR)
Inputs : VAR = variable name
Opt. Inputs : None.
Outputs : 1/0 if variable exists or not
Opt. Outputs: None.
Keywords : None.
Calls : None.
Common : None.
Restrictions: None.
Side effects: None.
Category : Useful stuff
Prev. Hist. : None.
Written : Dominic Zarro (ARC)
Version : Version 1.0, 18 September 1993
(See pro/external/exist.pro)
procedure fidgit digitally (i.e., using Darwin-given fingers) fit a DEM to a given spectrum syntax fidgit,x,y,lstr,cstr,ysigma=ysigma,dem=dem,logt=logt,abund=abund,$ abfrac=abfrac,ldbdir=ldbdir,cdbdir=cdbdir,ceroot=ceroot,lsf=lsf,$ wdem=wdem,wspc=wspc,verbose=verbose,$ pres=pres,logP=logP,n_e=n_e,desig=desig,econf=econf,/allah,$ chidir=chidir,chifil=chifil,chidir=chidir,eqfile=eqfile,$ effar=effar,wvlar=wvlar parameters x [INPUT; required] wavelengths or channels of observed spectrum * if channels, must specify an RMF via the keyword LSF that matches the binning of the spectrum BEWARE: no checks are performed to verify that the supplied RMF is consistent with the input data y [INPUT; required] spectrum lstr [I/O; required] line emissivity structure of the sort read in by RD_LINE() * if not given on input, will call RD_LINE() and generate it cstr [I/O; required] continuum emissivity structure of the sort read in by RD_CONT() * if not given on input, will call RD_LINE() and generate it keywords ysigma [INPUT] errors on Y; default is 1+SQRT(0.75+ABS(Y)) dem [I/O] DEM to start out with, will contain the final result on output logt [I/O] log(T [K]) at which DEM are defined. ignored on input if it doesn't match the size of DEM abund [INPUT] abundances * if not set, assumed to be from Anders & Grevesse 1989 abfrac [I/O] multiplicative factor by which to modify ABUND * if scalar, assumed to be metallicity * if vector, must match size of ABUND, else gets reset to 1 * default is 1 ldbdir [INPUT] path to line emissivity database to read in if LSTR is not given * default is '$CHIANTI' * uses !LDBDIR if set cdbdir [INPUT] path to continuum emissivity database to read in if CSTR is not given * default is '$CONT' * uses !CDBDIR if set ceroot [INPUT] prefix for continuum emissivity files * default is 'cie' * uses !CEROOT if set lsf [INPUT] line spread function * if scalar, the width of the line in bins used as a halfwidth for boxcar smoothing * if vector, assumed to be the kernel with which to smooth * if RMF structure (of the type read in from RD_OGIP_RMF()) then convolves the theoretical spectrum with this RMF wdem [INPUT] window number to display DEM in wspc [INPUT] window number to display spectrum in verbose [INPUT] controls chatter * uses !VERBOSE if set _extra [INPUT] allows specifying defined keywords to subroutines called by this program * RD_LINE: PRES, LOGP, N_E, DESIG, ECONF, ALLAH * RD_CONT: PRES, LOGP, N_E * FOLD_IONEQ: CHIFIL, VERBOSE * RD_IONEQ: CHIDIR, EQFILE * LINEFLX: EFFAR, WVLAR restrictions * requires subroutines -- DEMACS -- WHEE -- RD_LINE [FOLD_IONEQ, RD_IONEQ, READ_IONEQ] -- RD_CONT -- LINEFLX -- CONV_RMF -- INICON history vinay kashyap (Feb97) modified ion balance calcs (VK; 99May) completely rewritten (VK; Aug04)
(See pro/fitting/fidgit.pro)
Name: file_exist
Purpose: returns true(1) if any files match pattern=file
false(0) if no files found
Input Parameters:
file - file, pathname or file filter to check
Optional Keyword Parameters
direct - set if want full check (slower) - use if might be an
empty directory you are looking for
History: written by slf, 11/30/91
4-Sep-92 (MDM) - Modified to handle an array of input files
23-sep-93 (slf) - make sure null file returns false
(findfile count=1! when null requested)
17-Mar-94 (SLF) - use file_stat(/exist) for non wild card cases
April Fools'94 (DMZ) -- added check for :: in file name
(file_stat doesn't work across DECNET)
31-May-94 (DMZ) -- added check for VMS directory name input
(file_stat doesn't seem to work on VMS
dir names like '[ys.atest.soft]')
16-Aug-94 (DMZ) -- added another VMS patch to handle
case when user enters a VMS subdirectory name
without a file name.
6-Aug-97 (Zarro/GSFC) -- added check for nonstring input
30-may-99 (rdb) -- stop use of file_stat with WINDOWS
cause access violation - for wildcard search
8-Jan-1999 - S.L.Freeland - put in IDL version check since 5.3/IRIX
(at least) not backwardly compatible (must use findfile)
10-Jan-1999 - S.L.Freeland - extended 8-jan mod to all UNIX 5.3
15-Feb-2000 - S.L.Freeland - removed 5.3 changes (moved fix to file_stat)
(See pro/external/file_exist.pro)
Project : SOHO - SSW
Name : FILE_STAT()
Purpose : Vector version of FSTAT
Category : Utility, Operating_system
Explanation : Vector version of FSTAT
Syntax : Result = FILE_STAT( FILES )
Examples :
Inputs : FILES = List of files to check.
Opt. Inputs : None.
Outputs : None.
Opt. Outputs: None.
Keywords : EXIST = If set, then the returned values are whether the
files exist or not. This is the default behavior.
SIZE = If set, then the returned values are the sizes of the
files.
Calls : DATA_CHK
Common : None.
Restrictions: None.
Side effects: None.
Prev. Hist. : 11-Mar-1994 (SLF) Wanted faster file_exist function
History : Version 1, 11-Mar-1994, S. Freeland
Version 1.1 05-Jun-1998, J. Newmark changed loop to long
Version 1.2 15-Feb-2000, S.L.Freeland - work around RSI
Version 5.3 non-backwardly compatible change....
Version 1.3 10-Mar-2000, S.L.Freeland - fix 5.3 /SIZE bug
Contact : SFREELAND
Restrictions:
file size returned for directories under +5.3 is not valid
(See pro/external/file_stat.pro)
ROUTINE: findex
PURPOSE: Compute "floating point index" into a table using binary
search. The resulting output may be used with INTERPOLATE.
USEAGE: result = findex(u,v)
INPUT:
u a monitically increasing or decreasing 1-D grid
v a scalor, or array of values
OUTPUT:
result Floating point index. Integer part of RESULT(i) gives
the index into to U such that V(i) is between
U(RESULT(i)) and U(RESULT(i)+1). The fractional part
is the weighting factor
V(i)-U(RESULT(i))
---------------------
U(RESULT(i)+1)-U(RESULT(i))
DISCUSSION:
This routine is used to expedite one dimensional
interpolation on irregular 1-d grids. Using this routine
with INTERPOLATE is much faster then IDL's INTERPOL
procedure because it uses a binary instead of linear
search algorithm. The speedup is even more dramatic when
the same independent variable (V) and grid (U) are used
for several dependent variable interpolations.
EXAMPLE:
; In this example I found the FINDEX + INTERPOLATE combination
; to be about 60 times faster then INTERPOL.
u=randomu(iseed,200000) & u=u(sort(u))
v=randomu(iseed,10) & v=v(sort(v))
y=randomu(iseed,200000) & y=y(sort(y))
t=systime(1) & y1=interpolate(y,findex(u,v)) & print,systime(1)-t
t=systime(1) & y2=interpol(y,u,v) & print,systime(1)-t
print,f='(3(a,10f7.4/))','findex: ',y1,'interpol: ',y2,'diff: ',y1-y2
AUTHOR: Paul Ricchiazzi 21 Feb 97
Institute for Computational Earth System Science
University of California, Santa Barbara
paul@icess.ucsb.edu
REVISIONS:
(See pro/external/findex.pro)
function findscale returns the lengthscale in pixels at each point on the given curve syntax ls=findscale(curve,dim,/crunch,/half,pick=pick,choice=choice,eps=eps) parameters curve [INPUT; required] regularly gridded array of function values to be used to compute the length scales * if scalar, returns 0 * if 2D, -- use DIM to specify primary dimension -- compute lengthscales separately along each projection * if >2D, convert to 1D dim [INPUT; default=1] primary dimension in case of 2D array (e.g., if CURVE=CURVE(NX,NY), DIM=2 returns SCALE=SCALE(NY)) keywords crunch [INPUT] if set, and CURVE is 2D, collapses the array along the secondary dimension to generate 1D curve half [INPUT] if set, returns the half-scale pick [INPUT; default=0] if 2D, specifies how to combine the scales computed at the different cuts 0: pick the smallest scale 1: pick the largest scale 2: get the average choice [INPUT; default=0] what algorithm to use to find the scale? 0: MexicanHat wavelet 1: use inverse of 1st derivative 2: radius of curvature 3: stepped toggle eps [INPUT; default=1e-7] small number _extra [JUNK] ignore. here only to prevent crashing program. subroutines WVLT_SCALE [ROOFN] history vinay kashyap (Apr97) added CHOICE option 3 (VK; Feb03)
(See pro/findscale.pro)
function fitlines
widget-based procedure to fit Gaussians and Lorentzians (or
combinations, or other 3-parameter functions) to lines in a
spectrum and determine line fluxes. returns a structure of
the form
{POS,PERRP,PERRM,PERRC,$
FLX,FERRP,FERRM,FERRC,$
WDT,WERRP,WERRM,WERRC,$
THAW,TYPE,TIES,EPITHET,CONLEVX,CONLEV,CONSIGX,CONSIG,COMMENT}
All these are also available as output individually as keywords.
See keyword descriptions for what the variables mean.
syntax
fxstr=fitlines(x,y,ysig=ysig,funcs=funcs,/intens,dchi=dchi,$
pos=pos,perrp=perrp,perrm=perrm,perrc=perrc,$
flx=flx,ferrp=ferrp,ferrm=ferrm,ferrc=ferrc,$
wdt=wdt,werrp=werrp,werrm=werrm,werrc=werrc,$
thaw=thaw,type=type,epithet=epithet,ties=ties,$
conlev=conlev,consig=consig,comment=comment,$
histerix=histerix,$
/dumb,verbose=verbose,xsize=xsize,ysize=ysize,$
wid=wid,dynrng=dynrng,/posve,/negve,itmax=itmax,chithr=chithr,$
jumpup=jumpup,jumpdn=jumpdn,svdthr=svdthr,missing=missing,$
/noday,/notime,/nouser,/nopack,stacol=stacol,stacol=stacol,$
stathk=stathk,/nuthin)
parameters
x [INPUT; required] absissa, e.g., wavelength or energy
y [INPUT; required] count spectrum Y(X)
* size of Y must match that of X
keywords
ysig [INPUT] errors on Y
* default is sqrt(abs(Y)+0.75)+1.
funcs [INPUT] name of user defined function (actually a procedure)
that takes as input X and A (the model parameters), and returns
Y(X;A) and the partial derivatives dY/dA
* default is set to X3MODEL
* NOTE: if MPFIT is chosen as the optimization algorithm,
then "_f" is automatically appended to the name. This
is because MPFIT requires a _function_, not a procedure.
The program corresponding to X3MODEL is X3MODEL_F() and
that corresponding to LIBMODEL is LIBMODEL_F()
intens [INPUT] if set, Y(X) is assumed to be intensity, i.e.,
(counts|ergs|.../bin_unit/...)
* default is to take Y(X) to be (counts|ergs|.../bin/...)
* passed straight to FITLINES_EVENT
* WARNING: this has >not< been road-tested!
dchi [INPUT] delta-chisq threshold for projection errors
* default is 1.0
* may be changed within FITLINES_EVENT, but those changes
will not be propagated back.
pos [I/O] best-fit line positions
* length of vector shows number of components
perrp [I/O] 1-sided error on POS (POS+PERRP is upper bound)
perrm [I/O] 1-sided error on POS (POS-PERRM is lower bound)
perrc [I/O] DCHI threshold used in computing PERR?
flx [I/O] best-fit fluxes in the lines
ferrp [I/O] 1-sided error on FLX (FLX+FERRP is upper bound)
ferrm [I/O] 1-sided error on FLX (FLX-FERRM is lower bound)
ferrc [I/O] DCHI threshold used in computing FERR?
wdt [I/O] best-fit widths (sigma, core-radius, etc.) of the lines
werrp [I/O] 1-sided error on WDT (WDT+WERRP is upper bound)
werrm [I/O] 1-sided error on WDT (WDT-WERRM is lower bound)
werrc [I/O] DCHI threshold used in computing WERR?
* on input, FLX, WDT, and ?ERR? are forced to match the length
of POS: excess elements are thrown away, and insufficient
lengths are made up by padding with first elements, 1s, etc.
* on output ?ERRM are identical to ?ERRP >unless< the
projected errors have been calculated using ERORS, in
which case the values may be different
* on output, places where ?ERRC contain 0's are those where
the computed errors are 1-sigma formal curvature errors
thaw [I/O] integer array signaling frozen (0) or thawed parameter (1)
* refers to sequences of {POS,WDT,FLX} for each component in
a 3-parameter model (cf. X3MODEL) -- whatever goes for the
appropriate user-defined function.
* length must match 3*length(POS)
* default is to freeze defined input, thaw all new additions
type [I/O] type of model ('gauss', 'lorentz', etc.)
* default is 'gauss'
epithet [I/O] label for each component
* labels, obtained, for example, with IDLABEL
* Merriam-Webster> a characterizing word or phrase accompanying
or occurring in place of the name of a person or thing
ties [I/O] any ties between parameters?
conlev [I/O] the continuum that was taken out of the spectrum
* CONLEV must match the size of X and Y else ignored
consig [I/O] error on CONLEV
* default for CONSIG is sqrt(abs(CONLEV)+0.75)+1.
* NOTE: CONLEV and CONSIG are compressed using SPLAC in
the output that gets returned via the structure. That
structure therefore also has appropriate abscissae
CONLEVX and CONSIGX.
comment [OUTPUT] descriptive string
histerix [OUTPUT] a structure containing the state at each step
of the fitting process
* if explicitly set to 0, will not contain any output
* passed w/o comment to FITLINES_EVENT
oldstr [INPUT] old fitlines structure from which to start with
_extra [INPUT ONLY] use this to pass defined keywords to subroutines
PICKRANGE: XSIZE, YSIZE, WID, DYNRNG
LINEREM: POSve, NEGve
FIT_LEVMAR: ITMAX, CHITHR, DUMB
LEVMARQ: JUMPUP, JUMPDN, SVDTHR
MK_3MODEL: MISSING
MK_SLANT: ANGLE, BETAP
ERORS: VERBOSE
STAMPLE: NUTHIN, NODAY, NOTIME, NOUSER, NOPACK, STACOL,
STASIZ, STATHK
subroutines
FITLINES_EVENT
FMTSTUFF
PICKRANGE
FUNCS
(X3MODEL/LIBMODEL/MPFITFUN)
(MK_3MODEL, MK_GAUSS, MK_LORENTZ, MK_SLANT, MK_ROGAUSS, MK_POWLAM, MK_POLY1D)
LINEREM
SETCONT
ERORS
FIT_LEVMAR
ADJUSTIE
LEVMARQ
LMCOEFF
CURVE_FIT
SPLAC
WHICH
KILROY
WHEE
STAMPLE
PEASECOLR
IS_KEYWORD_SET
known bugs
cannot handle spectra with varying bin sizes
INTENS has not been tested
SPLAC keywords are not gracefully handled
help is available only in UNIX
history
vinay kashyap (Oct98)
added renormalization option (VK; Nov98)
added keywords INTENS,DCHI,PERRL,WERRL,FERRL,PERRC,WERRC,FERRC;
changed output structure format (VK; FebMM)
changed keywords PERR,WERR,FERR to PERRU,WERRU,FERRU; also the
corresponding output structure fields (VK; MarMM)
changed ?ERRU to ?ERRP and ?ERRL to ?ERRM to avoid confusion
with ERR[U,L] of ERORS (VK; MarMM)
moved QUIT to end of row, far from FIT; added ability to save
to disk; changed ?ERR[M,C] to be I/O; added TIES to output
structure (VK; JulMM)
added labeling capability via keyword EPITHET; reorganized
to combine 6th and 7th rows (VK; AugMM)
changed location of doc file to ardb ; allowed call to STAMPLE
(VK; DecMM)
changed default of ?ERRC to 0.0 (VK; JanMMI)
per Antonio Maggio's suggestion, removed keyword MULTIPLE from
call to WIDGET_LIST (VK; FebMMI)
changed suffix of help file from ".doc" to ".hlp" (VK; FebMMI)
added code to make it easier to switch between color tables
(VK; Oct02)
added monte-carlo errors as option (LL; Aug03)
added keyword HISTERIX (VK; Feb04)
added keyword OLDSTR (LL; Jul05)
updated for IDL5.6 keyword_set([0]) behavior change for vectors
(VK; 20Mar2006)
(See pro/fitting/fitlines.pro)
procedure fitlines_event widget event handler subroutine for FITLINES. see that routine for a description of variables, etc. only rudimentary consistency checks are carried out here. uses subroutine FMTSTUFF, which is included in this file. syntax fitlines_event,x,y,pos,wdt,flx,perr,werr,ferr,thaw,type,ties,epithet,$ conlev,consig,chisq,funcs,sigy,widg, perrm=perrm,werrm=werrm,$ ferrm=ferrm,perrc=perrc,werrc=werrc,ferrc=ferrc,/intens,rmfstr=rmfstr,$ histerix=histerix subroutines FMTSTUFF PICKRANGE FUNCS (X3MODEL/LIBMODEL) (MK_3MODEL, MK_GAUSS, MK_LORENTZ, MK_SLANT, MK_ROGAUSS, MK_POWLAM, MK_POLY1D) LINEREM SETCONT ERORS FIT_LEVMAR , CURVE_FIT , MPFITFUN ADJUSTIE LEVMARQ LMCOEFF STAMPLE KILROY WHEE PEASECOLR IS_KEYWORD_SET history vinay kashyap (Oct98) added renormalization option (VK; Nov98) various bug fixes and reorganizations (VK; Dec98) now allows fitting a single parameter (VK; Aug99) enhancements to include call to ERORS, setting YRANGEs, INTENSity units, etc. (VK; FebMM) avoid repeat calls to CW_BGROUP, saving tremendous headaches; bug correction CONTINUUM->AdjustNorm (VK; MarMM) bug correction -- undefined STYPE for new component; rudimentary changing parameter labels (VK; JunMM) moved QUIT to end of row, far away from FIT; added DUMP to disk (VK; JulMM) added error-bar plotting capability via VIEW->SetDefaults; pass x,y-range to setcont; added labeling capability via keyword EPITHET; merged 6th and 7th rows; bug correction with xcont in adjust norm (VK; AugMM) changed location of help file to ardb; added calls to STAMPLE (VK; DecMM) bug fixes: invisible error output, adding model deleted existing projected errors, streamlined saves (VK; JanMMI) per Antonio Maggio's suggestion, removed keyword MULTIPLE from call to WIDGET_LIST; improved color-scale setting for 24-bit consoles; replaced call to WHICH by call to ROUTINE_INFO (VK; FebMMI) various extra info messages (VK; AprMMI) changed y to z on lines 1012 and 1285, because.. y was inappropriate for some reason? (VK; Jun01) was crashing due to missing ZSIG in cont_acpt (VK; Aug01) bug was not passing predefined CONLEV to SETCONT (piecewise) (VK; Apr02) handle color tables in 24-bit displays (VK; Jun02) bug correction -- if all ties were deleted, then TIES was turning into and integer array (VK; Jul02) made changes to plotting so that hardcopy won't come out with annoying blank sheets (VK; Sep02) changed default colors to be compatible with PEASECOLOR (VK); added hooks into MPFIT (Liwei Lin/VK; Oct02) tied the left column of parameter values to changes in the main parameter list window on the right (VK; Apr03) added monte-carlo errors as option and added keyword RMFSTR (LL; Aug03) added keyword HISTERIX (VK; Jul03) bug correction: algo_type undefined unless fit is run first (LL; Aug04) updated for IDL5.6 keyword_set([0]) behavior change for vectors (VK; 20Mar2006) modified some widget instructions to be slightly clearer (VK; Aug06) ********************************************************************* subroutine fmtstuff takes in the model as described, and returns it properly formatted for display in the widget. parameters are named as before. only minimal consistency checks are carried out. usage fmtstuff,pos,wdt,flx,thaw,type,epithet,perr,werr,ferr,$ perrm=perrm,werrm=werrm,ferrm=ferrm,$ perrc=perrc,werrc=werrc,ferrc=ferrc,$ plist=plist,list=list,$ tiep=tiep,tiew=tiew,tief=tief,ties=ties,delcmp=delcmp keywords plist [OUTPUT] appropriately formatted model parameter list list [OUTPUT] denoting frozen/thawed components tiep [I/O] ties specially setting the range on positions tiew [I/O] ties specially setting the range on widths tief [I/O] ties specially setting the range on fluxes ties [I/O] rest of ties delcmp [ACTION] if set to component number, deletes component DELCOMP and all TIES that contain references to parameters in that component _extra [INPUT] pass defined keywords to STAMPLE (NODAY,NOTIME, NOUSER,NOPACK,STACOL)
(See pro/fitting/fitlines_event.pro)
script fitlines_undump restores the variables saved to disk from FITLINES_EVENT and calls FITLINES again. warning this will restore all saved variables from within FITLINES_EVENT. all of your existing variables with the same names will be overwritten. best to always start from a clean environment. save files written in IDL 5.4 are incompatible with earlier versions, but can be restored under certain circumstances using the CMSVLIB package of Craig Markwardt (http://cow.physics.wisc.edu/~craigm/idl/idl.html) This program has a "restore" command right at the start that may fail, but if the variables have already been loaded in, just type .SKIP and .CON, and bob's your uncle. usage fitsavfil='fitlines.save' .run fitlines_undump input FITSAVFIL name of IDL save file containing the variables dumped from within the FITLINES GUI output FITSTR structure containing the fit parameters history vinay kashyap (MMJul) cosmetic surgery (VK; DecMM) included keyword HISTERIX (VK; FebMMIV)
(See pro/scrypt/fitlines_undump.pro)
procedure fit_levmar uses the Levenberg-Marquardt method to find best-fit parameters syntax fit_levmar,x,y,a,yfunc,freeze=freeze,erra=erra,chisq=chisq,$ itmax=itmax,chithr=chithr,/dumb,ties=ties,vname=vname,$ jumpup=jumpup,jumpdn=jumpdn,svdthr=svdthr,funcs=f