大脸猫爱吃鱼

模板光谱库的处理

astronomy
字数统计: 1.1k阅读时长: 5 min
2018/12/16 Share

这是我的毕业论文的内容之1，这个内容一共分为两个部分，这是第一个部分，第二个部分是我毕业论文的主体。这一篇只是将我对得到的数据进行的一些预处理。
模板谱来自陈燕梅老师2012年发表在MNRAS上的文章。
一共有3个版本，如下：

1.“lib_spec”:the library of 25000 models in ASCIItext modewhich is generated by gene_lib_tim_cont_truncate.pro

  • first row：wavelengh (6900 points)
  • rest rows：flux of each model spectrum
  1. “lib_spec_flt.fits”： the library of model spectra in FITS image mode, single precision
  2. “lib_spec_dbl.fits”： the library of model spectra in FITS image mode, double precision

以上的文件内容均一致，但是在文件大小上fits文件比起ASCII文件要小很多，lib_spec文件大小为2.08GB，转为fits文件之后，单精度文件的大小变为了658MB，双精度的文件的大小变为了1.28GB。并且在使用IDL读取ASCII与fits文件上，fits文件的读取比起ACSII文件的读取要便捷很多，但是需要在函数库上加入天文学的函数库。
还得到了陈老师生成模板谱的IDL程序：gene_lib_tim_cont_truncate.pro

  • 通过阅读发现陈老师的这个程序产生的模板谱并没有卷积速度弥散，所以需要对数据进行再次处理。
    但是由于这个不是我自己写的，所以也就不贴上来了。
  • 进行卷积速度弥散之前有一些对数据的操作。
    下面是对数据进行插值，并且取出数据中的3550Å-5500Å的数据输出至文件"wave.txt",“spec.dat”
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
pro new
  niw = 6900
  nmod = 25000
  wave = make_array(niw,/float)
  flux = make_array(niw,nmod,/float)

  openr,lun,'lib_spec',/get_lun
  readf,lun,wave
  readf,lun,flux
  free_lun,lun
  
  dlogwv=1.d-4
  waveint0=alog10(3500.)
  waveint=waveint0+indgen(2050)*dlogwv
  waveint=10^waveint
  fluxint = make_array(2050,nmod,/float)
  
  for ii = 0L,nmod-1 do begin
  linterp,wave,flux[*,ii],waveint,fluxaa
  fluxint[*,ii]=fluxaa
  endfor
  
  ind = where(waveint ge 3550 and waveint le 5550,ct)
  covflux = make_array(ct,nmod,/float)
  covwave = waveint[ind]
  
  for ii=0,nmod-1 do begin
    covflux[*,ii]=fluxint[ind,ii]
  end
  print,ct
  ;print,covwave
  openw,lun,'wave.txt',/get_lun
  printf,lun,covwave,format='(1941(e12.5,1x))'
  free_lun,lun
  help,covwave
  
  openw,lun,'spec.txt',/get_lun
  for ii = 0L,nmod-1 do begin
    printf,lun,covflux[*,ii],format='(1941(e12.5,1x))'
  end
  free_lun,lun
  help,covflux
  
end

在输出了插值后的数据后，便需要对数据进行卷积速度弥散，下面是对数据进行卷积速度弥散的fortran程序。输出文件:“spec_vdisp.txt”
文件"spec_vdisp.txt"为卷积了速度弥散后的数据，"wave.txt"为包含了波长点的文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
program vdisp
    implicit none
    integer :: m,n,i,j
    parameter(m=25000)
    parameter(n=1941)
    real :: wave(n),vdisp_add,bc03_pix ,bc03_vdisp
    real :: spec(n),vd,sigma_pix
    !read file
    !read vdisp
    open(11,file='wave.txt')
        read(11,*) wave
    close(11)
    bc03_pix = 70.0    
    bc03_vdisp =  75.0
100 format(1941E15.5E2)
    open(12,file='vd.dat')
    open(13,file='spec.txt')
    open(14,file='spec_vdisp.txt')
    do i=1,m
        read(12,*) vd
        read(13,*) spec
        call GAUSSIAN_V_DISP(wave,spec,n,vd)
        write(14,100) spec
    enddo
    close(14)
    close(13)
    close(12)
end
     
    SUBROUTINE GAUSSIAN_V_DISP(X,Y,N,SIGMA)
!	Applies a gaussian filter to the sed in (X,Y) to reproduce
!	the effect of the velocity dispersion of stars in a galaxy
!	Array declaration
	real x(n),y(n),z(10000),u(1000),g(1000)
!	Enter sigma in Km/s
!	Speed of light in Km/s
	c=3.00E5
!	Number of sigmas to deviate from v=0
	m=6
!	Copy array y to array z
	do i=1,n
	z(i)=y(i)
	enddo
	do i=1,n
	xmax=c*x(i)/(c-m*sigma)
	call locate(x,n,xmax,i1)
	m2=i1+1
	m1=2*i-m2
	if (m1.lt.1) m1=1
	if (m2.gt.n) m2=n
	k=0
	do j=m2,m1,-1
	v=(x(i)/x(j)-1.)*c
	k=k+1
	u(k)=v
	g(k)=z(j)*gauss(v,0.,sigma)
	enddo
	y(i)=trapz1(u,g,k)
	enddo
	return
	end

      SUBROUTINE locate(xx,n,x,j)
      INTEGER j,n
      REAL x,xx(n)
      INTEGER jl,jm,ju
      jl=0
      ju=n+1
10    if(ju-jl.gt.1)then
        jm=(ju+jl)/2
        if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then
          jl=jm
        else
          ju=jm
        endif
      goto 10
      endif
      j=jl
      return
      END
!  (C) Copr. 1986-92 Numerical Recipes Software '5s!61"25:5<.



	FUNCTION GAUSS(X,X0,SIGMA)
!	Returns value of gaussina function normalized to area = 1
!	Normalization constant
	data a/0./
	if (a.eq.0) then
		pi=4.*atan(1.)
		a=1./sqrt(2.*pi)
	endif
        c=a/sigma
        gauss=c*exp(-(x-x0)**2/sigma**2/2)
        return
        end
        
	FUNCTION TRAPZ1 (X,Y,N)
	REAL X(N),Y(N)
	TRAPZ1=0.
	IF (N.LE.1) RETURN
!	IF (N.EQ.1) GOTO 2
	DO 1 J=2,N
1	TRAPZ1= TRAPZ1 + ABS(X(J)-X(J-1))*(Y(J)+Y(J-1))/2.
	RETURN
2	TRAPZ1=Y(1)*X(1)/2.
	RETURN
	END

我已经是一只废喵了
我已经是一只废喵了

原文作者：Miao cc

原文链接：http://www.bilibiliwiki.com/模板光谱库的处理/

发表日期：December 16th 2018, 11:19:10 am

更新日期：January 18th 2019, 8:31:03 pm

版权声明：本文采用知识共享署名-非商业性使用 4.0 国际许可协议进行许可

CATALOG
Total : 15
2019
2018
astronomy Linux python IDL
Linux