|
一、参考
1、IDL
Reference
Coyote's Guide to IDL Programming
https://www.harrisgeospatial.com/docs/routines-1.html
https://ips.ucsd.edu/help/smei/SMEI_list.html
(Fortran)
https://ips.ucsd.edu/help/ssw/SSW_list.html
(IDL)
https://www.cnblogs.com/alex-bn-lee/category/725300.html?page=1
IDL控制语句：https://blog.csdn.net/weixin_43955546/article/details/105661962
2、Formatting IDL Graphics Symbols and Lines
https://www.harrisgeospatial.com/docs/formattingsymsandlines.html#Color
clip=[0,255]:最暗到最亮 loadct: 0~31 cgloadct:
0~40 /brewer:颜色内外相反 /reverse,
/brewer的情况下，cgloadct: 0~31
plot和oplot绘图是，用16进位颜色顺序是'蓝绿红'，color='0000ff'xL（red）；color='0050ff'xL（orange）；color='00ffff'xL（yellow）；color='00ff00'xL（green）；color='00ff00'xL（blue）；color='ffff00'xL（cyan）；color='ff00ff'xL（pink）……当然，可以用数字或字母使得颜色亮度逐渐增加，比如：'0000ff'xL为全红最亮，'000099'xL为红色较暗，'000050'xL为红色更暗，'0000ff'xL为全黑了；'505050'xL到'707070'xL到'909090'xL则为灰色越来越亮；'004050'xL到'007000'xL到'009000'xL则为深绿越来越亮。
FORTRAN-Style Format Codes
https://www.l3harrisgeospatial.com/docs/format_codes_fortran.html
fortran格式代码的语法如下:
[n] C[+][-][width]
n：格式的重复数(n≥1)。若不指定n，则为1。
C：格式代码。
+：可选，指定输出的正数应该带有“+”前缀。“+”标志仅对数字格式代码有效。通常，只有负数输出带有符号前缀。非十进制数字代码(B、O和Z)允许“+”标志，但忽略它。
-：可选，指定字符串或数字值在输出时应左对齐。通常，输出是右对齐的。
width：宽度规范(可选)。宽度规范和默认值是特定于格式代码的。
3、PLOT
https://www.harrisgeospatial.com/docs/PLOT_Procedure.html
Over-plot
data points with accompanying X or Y error bars.
oploterror, [ x,] y, [xerr], yerr, [ /NOHAT,
HATLENGTH= , ERRTHICK =, ERRSTYLE=, ERRCOLOR =,
/LOBAR, /HIBAR, NSKIP = , NSUM = , /ADDCMD, ...
OPLOT keywords ]
https://www.l3harrisgeospatial.com/docs/oploterror.html
4、PLOTHIST
https://www.harrisgeospatial.com/docs/plothist.html
5、cgContour
http://www.idlcoyote.com/idldoc/cg/cgcontour.html
6、WHERE
Result =
WHERE( Array_Expression [, Count] [,
COMPLEMENT=variable] [, /L64] [, NCOMPLEMENT=variable]
)
这里Count为符合条件的数量；COMPLEMENT为不符合条件的脚码；NCOMPLEMENT为不符合条件的数量；/L64在默认条件下返回的是32位整数，加了/L64则返回64位整数。
可以简写，比如简写为：result=where(array, c, com=a, ncom=b)，或简写为：result=where(array,
c, co=a, nc=b)
https://www.l3harrisgeospatial.com/docs/WHERE.html
https://northstar-www.dartmouth.edu/doc/idl/html_6.2/WHERE.html
7、REBIN——重新调整数组或向量的尺寸，注意：调整后的维数必须为原数组维数的整数倍
Result =
REBIN( Array, D1 [, ..., D8] [, /SAMPLE] )
IDL>
a=[1.,2.,3.,4.]
IDL> b=rebin(a,12)
IDL> print,b
1.00000 1.33333 1.66667 2.00000 2.33333 2.66667
3.00000 3.33333 3.66667 4.00000 4.00000 4.00000
IDL> b=rebin(a,8)
IDL> print,b
1.00000 1.50000 2.00000 2.50000 3.00000 3.50000
4.00000 4.00000
IDL> c=rebin(b,4)
IDL> print,c
1.25000 2.25000 3.25000 4.00000
https://northstar-www.dartmouth.edu/doc/idl/html_6.2/REBIN.html
8、LINTERP——从一个网格线性插值一维数据到一个新的网格。LINTERP的结果在数值上与IDL内禀结果是等价的！
LINTERP,
Xtab, Ytab, Xint, Yint, [MISSING =, /NoInterp ]
IDL>y1=fltarr(10)
IDL>x=[0,1.45,2.37,3.56,4.28,5.10,6.89,7.31,8.55,9.61]
IDL>y=x^2
IDL>x1=findgen(10)
IDL>linterp,x,y,x1,y1
IDL>print,x1,y1
0.00000 1.00000 2.00000 3.00000 4.00000 5.00000
6.00000 7.00000 8.00000 9.00000
0.00000 1.45000 4.20350 9.35280 16.1232 25.0720
36.8010 49.0341 64.3795 81.2745
https://www.l3harrisgeospatial.com/docs/linterp.html
9、波长数组和取log后的对应数组
λ1、λ2、λ3、……、λn
（波长间隔为1Å）
y1、y2、y3、……、yn
（对应的取log后的值）
y2
＝ y1+lgλ2－lgλ1
y3
＝ y2＋lgλ3－lgλ2
＝ y1+lgλ2－lgλ1＋lgλ3－lgλ2＝y1+lgλ3－lgλ1
y4
＝ y3＋lgλ4－lgλ3
＝ y1+lgλ3－lgλ1＋lgλ4－lgλ3＝y1+lgλ4－lgλ1
…………
yn
＝y1+lgλn－lgλ1
例如λ从3400Å到8500Å，求其取log后对应的波长
wave
= alog10(3400) + dindgen(4001) * (alog10(8500) -
alog10(3400))/4d3
flux
= dindgen(4001) * 0.
10、添加lts_fits的几个PRO：
http://www-astro.physics.ox.ac.uk/~mxc/software/#lts
https://pypi.org/project/ltsfit/
解压后文件夹放在/usr/local/idl71/lib/lts_fits，然后在IDL_PATH中添加lts_fits路径目录，IDL_PATH=/usr/local/idl71/lib:/usr/local/idl71/lib/astro/pro:/usr/local/idl71/lib/lts_fits
然后IDL>
print, file_which('lts_linefit.pro')就找到了lts_linefit.pro所在的路径：
/usr/local/idl71/lib/lts_fits/lts_linefit.pro
简单地：IDL>which('lts_linefit.pro')即可
CALLING
SEQUENCE:
lts_linefit, x, y, sigx, sigy, par, sig_par,
chi_sq, PLOT=plot, RMS=rms, BAD=bad, /BAYES,
CLIP=clip, EPSY=epsy, FRAC=frac, GOOD=good,
OVERPLOT=overplot, PIVOT=pivot, TEXT=text
"LTS_LINEFIT
program described in Cappellari et al. (2013,
MNRAS, 432, 1709), which combines the Least
Trimmed Squares robust technique of Rousseeuw &
van Driessen (2006) into a least-squares fitting
algorithm which allows for errors in both
variables and intrinsic scatter."
11、计算光度距离——lumdist(z)，单位：Mpc
result =
lumdist(z, [H0 = 70 , k = 0, Omega_M = 0.3,
Lambda0 = 0.7, q0 = -0.55, /SILENT])
方括号中为默认值，如果lumdist(z)，则其它参数默认为上面的值
，光度距离结果以Mpc为单位。
http://www.harrisgeospatial.com/docs/lumdist.html
https://idlastro.gsfc.nasa.gov/ftp/pro/astro/lumdist.pro
如：计算SDSS目标光度
L = Flux·4πD2 = F×10-17
erg/cm2/s ×4π×(lumdist(z)×106×3.26×365×24×3600×3×1010)2
cm2
= F×10-17
erg/cm2/s ×4π×(lumdist(z)×3.086×1024)2
cm2 =···erg/s = ···×10-7
J/s = ···W/(3.828×1026 W) = ···L⊙
其中：L⊙
= 3.828×1026 W
12、半峰全宽（FWHM，Full Width at Half
Maximum）
FWHM≈2.3548σ
参考：https://mathworld.wolfram.com/GaussianFunction.html
13、等值宽度（Equivalent Width，EW）
等值宽度，即与吸收（或发射）谱线轮廓和连续谱之间所包围的面积相当的高度为1的矩形的宽度。
参考：https://astronomy.swin.edu.au/cosmos/E/Equivalent+Width
14、关于赤道坐标系转换
（1）RADEC——To convert RA and Dec
from decimal to sexagesimal units.
radec, ra, dec, ihr, imin, xsec, ideg, imn, xsc,
[/HOURS}
ra - Right ascension, scalar or vector, in
DEGREES unless the /HOURS keyword is set
dec - declination in decimal DEGREES, scalar or
vector, same number of elements as RA
/HOURS - if set, then the input righ ascension
should be specified in hours instead of degrees.
https://www.l3harrisgeospatial.com/docs/radec.html
（2） JPRECESS——Precess astronomical coordinates
from B1950 to J2000
jprecess, ra, dec, ra_2000, dec_2000, [ MU_RADEC
= , PARALLAX = RAD_VEL =, EPOCH = ]
Example：
RA(1950) = 13h 39m 44.526s Dec(1950) = 8d 38'
28.63'' Mu(RA) = -.0259 s/yr Mu(Dec) = -.093
''/yr
IDL> mu_radec = 100D* [ -15D*.0259, -0.093 ]
IDL> ra = ten(13,39,44.526)*15.D
IDL> dec = ten(8,38,28.63)
IDL> jprecess, ra, dec, ra2000, dec2000,
mu_radec = mu_radec
IDL> print, adstring(ra2000, dec2000,2)
===> 13h 42m 12.740s +08d 23' 17.69"
https://www.l3harrisgeospatial.com/docs/jprecess.html
15、K-correction
If you want to convert apparent
magnitudes in band R to absolute magnitudes in
band Q, you need to calculate the K-correction,
which is defined by the equation:
mR = MQ + DM(z) + KQR(z),
where mR is the apparent magnitude, MQ
is the absolute magnitude, DM(z) is the distance
modulus, accounting for the angular diameter
distance and cosmological surface-brightness
dimming, and KQR(z) is the
K-correction. By absolute magnitude we mean: the
apparent magnitude in band Q the object would
have if it were observed at rest, 10 pc away,
using an aperture that contains its total flux.
The distance modulus accounts for the difference
between an object's actual distance and 10 pc.
The K-correction accounts for the fact that you
observed a redshifted galaxy in band R but the
absolute magnitude requires a rest-frame
observation in band Q. Obviously the difference
between the fluxes observed in different
bandpasses is fully determined by the galaxy SED
and the description of the bandpasses.
https://arxiv.org/pdf/astro-ph/0210394.pdf
http://kcorrect.org/
二、拟合
1、linfit线性拟合（适应于只有X、Y，但没有误差值的两组数据）：
;
有一对n个元素的向量：
X = [-3.20, 4.49, -1.66, 0.64, -2.43, -0.89,
-0.12, 1.41,
2.95, 2.18, 3.72, 5.26]
Y = [-7.14, -1.30, -4.26, -1.90, -6.19, -3.98,
-2.87, -1.66,
-0.78, -2.61, 0.31, 1.74]
; 定义其中一个向量的泊松测量误差为：
measure_errors = SQRT(ABS(Y))
; 计算模型参数A和B，并输出结果为：
result = LINFIT(X, Y, MEASURE_ERRORS=error,
sigma=sig, yfit=yfit, chisq=chi2, prob=prob)
PRINT, result
IDL prints:
-3.16574 0.829856
IDL>
print, sig
0.282226 0.0991508
y = (a±σa)
+ (b±σb)x = (res[0]±sig[0]) + (res[1]±sig[1])x
https://www.harrisgeospatial.com/docs/LINFIT.html
2、fitexy线性拟合（适应于包含x、y的误差
的数据）：
Linear
Least-squares approximation in one-dimension (y
= a + b*x), when both x and y data have errors
FITEXY, x, y, A, B, X_SIG=xerr, Y_SIG=yerr, [sigma_A_B,
chi_sq, q, TOL=]
为了避免X和Y由于数量级差别太大带来的问题，通常需要对不同数量级的数据进行归一化处理后再进行拟合。（参考Press和Teukolsky的"Numerical
Recipes"一书中的"Numerical Recipes"一栏）
归一化过程如下：
xx = (x - xm) / xs, sigx = x_sigma / xs，式中：xm =
mean(x)， xs = stddev(x)
yy = (y - ym) / ys, sigy = y_sigma / ys，式中：ym =
mean(y)， ys = stddev(y)
https://www.l3harrisgeospatial.com/docs/fitexy.html
3、ladfit最小绝对偏差拟合：
; 有一对n个元素的向量：
X = [-3.20, 4.49, -1.66, 0.64, -2.43, -0.89,
-0.12, 1.41,
2.95, 2.18, 3.72, 5.26]
Y = [-7.14, -1.30, -4.26, -1.90, -6.19, -3.98,
-2.87, -1.66,
-0.78, -2.61, 0.31, 1.74]
; 将X值进行升序排序，并对Y值也进行升序排序以匹配X中元素的新顺序：
XX = X[SORT(X)]
YY = Y[SORT(X)]
; 计算模型参数A和B:
PRINT, LADFIT(XX, YY)
IDL prints:
-3.15301 0.930440
https://www.harrisgeospatial.com/docs/LADFIT.html
4、poly_fit线性拟合：
; Define
an 11-element vector of independent variable
data:
X = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
0.8, 0.9, 1.0]
; Define an 11-element vector of dependent
variable data:
Y = [0.25, 0.16, 0.09, 0.04, 0.01, 0.00, 0.01,
0.04, 0.09,
0.16, 0.25]
; Define a vector of measurement errors:
measure_errors = REPLICATE(0.01, 11)
; Compute the second degree polynomial fit to
the data:
result = POLY_FIT(X, Y, 2, MEASURE_ERRORS=measure_errors,
$
SIGMA=sigma)
; Print the coefficients:
PRINT, 'Coefficients: ', result
PRINT, 'Standard errors: ', sigma
IDL prints:
Coefficients: 0.250000 -1.00000 1.00000
Standard errors: 0.00761853 0.0354459 0.0341395
https://www.harrisgeospatial.com/docs/POLY_FIT.html
5、高斯拟合：
（1）Gauss1——计算给出x值、平均值（中心值μ）、σ、面积的高斯曲线y值
形式：yvals =
gauss1(xvals, [mean, sigma, area], skew=skew)
举例：
p = [2.2D,
1.4D, 3000.D]
x = dindgen(200)*0.1 - 10.
y = gauss1(x, p)
计算等间距(间距为0.1)情况下高斯函数的值。高斯函数的中心值为2.2，标准差为1.4，总面积为3000。
http://www.harrisgeospatial.com/docs/gauss1.html
（2）Gauss1p——计算给出x值、平均值（中心值μ）、σ、面积的高斯曲线y值
形式：gauss1,
xvals, [mean, sigma, area], yvals, skew=skew
举例：
p = [2.2D,
1.4D, 3000.D]
x = dindgen(200)*0.1 - 10.
gauss1p, x, p, y
http://www.harrisgeospatial.com/docs/gauss1p.html
（3）mpfitfun——对IDL函数执行Levenberg-Marquardt最小二乘拟合
形式：parms = mpfitfun(myfunct, x, y, err,
start_params, ...)
用户自定义函数应按如下方式：
function myfunct, x, p
; The independent variable is x
; Parameter values are passed in "p"
ymod = ... computed model values at X ...
return, ymod
END
举例：
; First, generate some synthetic data
npts = 200
x = dindgen(npts) * 0.1 - 10. ; Independent
variable
yi = gauss1(x, [2.2D, 1.4, 3000.]) ; "Ideal" Y
variable
y = yi + randomn(seed, npts) * sqrt(1000. + yi);
Measured, w/ noise
sy = sqrt(1000.D + y) ; Poisson errors
; Now fit a Gaussian to see how well we can
recover
p0 = [1.D, 1., 1000.] ; Initial guess (cent,
width, area)
p = mpfitfun('GAUSS1', x, y, sy, p0) ; Fit a
function
print, p
生成具有高斯峰和泊松统计不确定度的数据集合，然后将相同的函数拟合到数据(具有不同的初始参数)，以视得到的数据接近程度。
http://www.harrisgeospatial.com/docs/mpfitfun.html
https://pages.physics.wisc.edu/~craigm/idl/mpfittut.html#download
http://www.eg.bucknell.edu/physics/ASTR201/IDLTutorial/tutorial_05.html
附：结构体http://www.weixueyuan.net/a/93.html
|