function vac_to_air,vac1 ;--------------------------------------------------------------------------- ; name: ; vac_to_air ; purpose: ; returns an air wavelength given a vacuum wavelength in nm ; comments: ; see Morton (1991, ApJS, 77, 119) ;--------------------------------------------------------------------------- vac=double(vac1) res=1d-1/(1.0002735182d-1/vac+131.4182d-3/vac^3+2.76249d3/vac^5) ; nm return,res end function dfdvac,vac ;--------------------------------------------------------------------------- ; name: ; dfdvac ; purpose: ; analytic derivative of vac_to_air with respect to the argument vac in nm ; comments: ; see Morton (1991, ApJS, 77, 119) ;--------------------------------------------------------------------------- a=1.0002735182d-1 b=131.4182d-3 c=2.76249d3 res=1d-1*(a/vac^2+3d0*b/vac^4+5d0*c/vac^6)/(a/vac+b/vac^3+c/vac^5)^2 return,res end function air_to_vac,air1 ;--------------------------------------------------------------------------- ; name: ; air_to_vac ; purpose: ; returns a vacuum wavelength given an air wavelength in nm ; comments: ; see Morton (1991, ApJS, 77, 119) ; uses Newtown-Raphson method ;--------------------------------------------------------------------------- air=double(air1) nair=n_elements(air) ; first guess vac=air for iair=0,nair-1 do begin while 1 do begin vac0=vac[iair] vac[iair]=vac0-(air[iair]-vac_to_air(vac0))/(-dfdvac(vac0)) if (abs(vac0/vac[iair]-1d0) lt 1d-9) then break end endfor return,vac end function hunt,grid,ngrid,val,count,linear=linear ;--------------------------------------------------------------------------- ; name: ; hunt ; purpose: ; find up to two nodes to the left and up to two nodes ; to the right of the desired point ; comments: ; none yet ;--------------------------------------------------------------------------- if (n_elements(linear) eq 0l) then is_linear=0l else is_linear=linear distance=min(abs(grid-val),ilo) distance=grid[ilo]-val if (distance eq 0d) then begin result=[ilo] count=1l return,result endif if (distance gt 0d) then ilo-=1l ihi=ilo+1l if (is_linear) then result=[ilo,ihi] else result=[ilo-1l,ilo,ihi,ihi+1l] index=where(result ge 0l,/null,count) if count eq 0l then stop result=result[index] index=where(result le ngrid-1l,/null,count) if count eq 0l then stop result=result[index] count=long(n_elements(result)) return,result end function getflux,teff_in,lgg_in,feh_in,vsini_in,$ vbroad_in,wl_air_in,line_in,grid=grid ;--------------------------------------------------------------------------- ; name: ; getflux ; description: ; returns rotationally-broadened, gaussian-broadened, ; normalised flux profiles for Halpha/beta/gamma ; example: ; see example.pro ;--------------------------------------------------------------------------- ; check the grid if (n_elements(grid) eq 0l) then begin ; no grid provided ; try to read it from current directory file=file_search('flux_3d.fits',count=nfile) if nfile eq 0l then begin print,'Grid not provided -- stopping!' stop endif else grid=mrdfits('flux_3d.fits',1) endif ; check line line=['hgamma','hbeta','halpha'] iline=where(line eq line_in,/null,count) if count ne 1l then begin print,'Line must be halpha, hbeta, or hgamma -- stopping!' stop endif else iline=iline[0l] ; vacuum wavelengths wl_in=air_to_vac(wl_air_in) ; input stellar parameters star={teff:teff_in,lgg:lgg_in,feh:feh_in,vsini:vsini_in,$ vbroad:vbroad_in,wl:wl_in} ind_teff=hunt(grid.teff,grid.nteff,star.teff,nteff,/linear) ind_lgg=hunt(grid.lgg,grid.nlgg,star.lgg,nlgg,/linear) ind_feh=hunt(grid.feh,grid.nfeh,star.feh,nfeh,/linear) ind_vsini=hunt(grid.vsini,grid.nvsini,star.vsini,nvsini) ind_vbroad=hunt(grid.vbroad,grid.nvbroad,star.vbroad,nvbroad) teff=reform(grid.teff[ind_teff]) lgg=reform(grid.lgg[ind_lgg]) feh=reform(grid.feh[ind_feh]) vsini=reform(grid.vsini[ind_vsini]) vbroad=reform(grid.vbroad[ind_vbroad]) ; space for the result nwl_in=n_elements(star.wl) result=dblarr(nwl_in) for jwl=0l,nwl_in-1l do begin ind_wl=hunt(reform(grid.wl[*,iline]),grid.nwl,star.wl[jwl],nwl) wl=reform(grid.wl[ind_wl,iline]) b=dblarr(nlgg) c=dblarr(nfeh) d=dblarr(nvsini) e=dblarr(nvbroad) f=dblarr(nwl) for iwl=0l,nwl-1l do begin for ivbroad=0l,nvbroad-1l do begin for ivsini=0l,nvsini-1l do begin for ifeh=0l,nfeh-1l do begin for ilgg=0l,nlgg-1l do begin a=grid.flux[*,ind_lgg[ilgg],ind_feh[ifeh],$ ind_vsini[ivsini],ind_vbroad[ivbroad],ind_wl[iwl],iline] a=a[ind_teff] if nteff eq 1l then b[ilgg]=reform(a) else $ b[ilgg]=interpol(a,teff,star.teff) endfor ; lgg if nlgg eq 1l then c[ifeh]=reform(b) else $ c[ifeh]=interpol(b,lgg,star.lgg) endfor ; feh if nfeh eq 1l then d[ivsini]=reform(c) else $ d[ivsini]=interpol(c,feh,star.feh) endfor ; vsini if nvsini eq 1l then e[ivbroad]=reform(d) else begin is_spline=(nvsini ge 4l) e[ivbroad]=interpol(d,vsini,star.vsini,spline=spline) endelse endfor ; vbroad if nvbroad eq 1l then f[iwl]=reform(e) else begin is_spline=(nvbroad ge 4l) f[iwl]=interpol(e,vbroad,star.vbroad,spline=spline) endelse endfor ; wl if nwl eq 1l then result[jwl]=reform(f) else begin is_spline=(nwl ge 4l) result[jwl]=interpol(f,wl,star.wl[jwl],spline=spline) endelse endfor ; input wls return,result end