Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.naic.edu/~rminchin/idl/masveltrans.pro Дата изменения: Wed Jun 18 18:13:05 2014 Дата индексирования: Sun Apr 10 06:12:12 2016 Кодировка:

pro masveltrans,b,bout,bary=bary,lsr=lsr,topo=topo,resamp=resamp,spline=spline,lsquadratic=lsquadratic,quadratic=quadratic,errcode=errcode,rest=rest,redshift=redshift,vel=vel,exact=exact

; Transform frequency from topocentric to another frame
; Uses formulae from http://www.gb.nrao.edu/~fghigo/gbtdoc/doppler.html

; Use '/exact' to change channel widths as well as velocities.

nrows=n_elements(b)
crval1 = b.h.crval1
cdelt1 = b.h.cdelt1
specsys=b.h.specsys
c=299792458D
bout=b
if (n_elements(resamp) eq 0) then resamp = 1
if (n_elements(bary) eq 0) then bary = 0
if (n_elements(lsr) eq 0) then lsr = 0
if (n_elements(topo) eq 0) then topo = 0
if (n_elements(rest) eq 0) then rest = 0
if (n_elements(vel) eq 0) then vel = 0
if (n_elements(exact) eq 0) then exact = 0
errcode = 0



if (rest) then begin
myvel=vel
if (n_elements(redshift) ne 1) then begin
redshift = b.h.req_vel
veltype = b.h.req_vel_type
if (strcmp(veltype[0],'v',1,/fold_case)) then vel = 1 else vel = 0
endif
if (myvel) then redshift = redshift * 1000
if (vel) then redshift = redshift/c
restfactor = 1.0 + redshift
endif else restfactor = 1.0



if ((topo + bary + lsr) EQ 0) then begin
print, 'Error: no target coordinate system specified!'
errcode = 1
endif else if ((topo + bary + lsr) GT 1) then begin
print, 'Error: only one target coordinate system may be specified!'
errcode = 1
endif



; Look for input specsys
; If LSR or Barycent, set topo to convert to topocentric before
; converting back out

if (strcmp(b[0].h.specsys,'TOPOCENT',/fold_case)) then begin
oldsys = 'topo'
if (topo) then begin
print, 'Error: Requested transform is from topocentric to topocentric!'
errcode = 1
endif
endif else if (strcmp(b[0].h.specsys,'BARYCENT',/fold_case)) then begin
oldsys = 'bary'
if (bary) then begin
print, 'Error: Requested transform is from barycentric to barycentric!'
errcode = 1
endif
topo = 1
endif else if (strcmp(b[0].h.specsys,'LSRK',/fold_case)) then begin
oldsys = 'lsr'
if (lsr) then begin
print, 'Error: Requested transform is from LSR to LSR!'
errcode = 1
endif
topo = 1
endif else begin
print,'Error: specsys not recognized: ',b[i].h.specsys
errcode = 1
endelse

if (errcode EQ 0) then begin
for i = 0, nrows-1 do begin
freqs = masfreq(b[i].h) * 1E6
nchans = b[i].nchan
oldcrval1 = b[i].h.crval1
relbary = b[i].h.vel_bary/c
rarad = b[i].h.crval2*!pi/180.d
decrad = b[i].h.crval3*!pi/180.d
radec3vec = anglestovec3(rarad,decrad)
rellsr = vellsrproj(radec3vec,relbary)


if (topo) then begin
if (oldsys EQ 'bary') then begin
relold = relbary
endif else if (oldsys EQ 'lsr') then begin
relold = rellsr
endif else begin
print, 'something weird has happened!'
errcode = 1
relold = 0
endelse

tfreq = freqs * sqrt(1 - relold^2)/(1 - relold)
tcrval1 = oldcrval1 * sqrt(1 - relold^2)/(1 - relold)
relframe = 0
newspecsys = 'TOPOCENT'

endif else begin
tfreq = freqs
tcrval1 = oldcrval1

endelse

if (bary) then begin
relframe = relbary
newspecsys='BARYCENT'
endif else if (lsr) then begin
relframe = rellsr
newspecsys='LSRK'
endif

bfreq = tfreq * (1 - relframe) / sqrt(1 - relframe^2)
newcrval1 = tcrval1 * (1-relframe)/sqrt(1-relframe^2)
specsys[i] = newspecsys

if (resamp) then begin
bfreqresamp = freqs - (oldcrval1 - newcrval1)
npols = b[i].npol
resampled = fltarr(nchans,npols)
for j = 0,npols - 1 do begin
resampled[*,j] = interpol(b[i].d[*,j],bfreq,bfreqresamp,spline=spline,lsquadratic=lsquadratic,quadratic=quadratic)
endfor
bout[i].d = resampled
endif
crval1[i] = newcrval1
endfor
endif


bout.h.crval1 = crval1*restfactor

if (exact) then begin
bout.h.cdelt1 = cdelt1*restfactor
endif else begin
bout.h.cdelt1 = cdelt1
endelse
bout.h.specsys = specsys

return
end