模板光谱库的处理

这是我的毕业论文的内容之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
2. “lib_spec_flt.fits”： the library of model spectra in FITS image mode, single precision
3. “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”

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,covfluxend

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

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.0100 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+110    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,N1   TRAPZ1= TRAPZ1 + ABS(X(J)-X(J-1))*(Y(J)+Y(J-1))/2.    RETURN2   TRAPZ1=Y(1)*X(1)/2.    RETURN    END

我已经是一只废喵了