Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.mrao.cam.ac.uk/~bn204/alma/memo-turb/uvblow.py Дата изменения: Fri Aug 10 20:17:36 2012 Дата индексирования: Tue Oct 2 10:12:12 2012 Кодировка: Поисковые слова: hydrogen

# Bojan Nikolic ,
#
# Initial version October 2007
# Revised 2008
#
# Blow phase screens across UV data

import os
import math
from itertools import izip

import uvfits
import numpy
import numarray
import pyfits

def GenWindVector(fin,
ws_px,
toffset=0.0,
dc="x",
vismask=None):
"""Generate the wind displacement vector list

:param ws_px: windspeed in pixels per second

:param toffset: time offset (in seconds) to add to all times (so we don't
start at begining of screen)

:param dc: direction of wind ("x" or "y")

:return: Two dimensional vector to faciliate other blow directions
in later implementations.
"""

if vismask is None:
ts = uvfits.IntegrationTimeSeries(fin) + toffset
res=numpy.zeros( ( uvfits.NVis(fin), 2))
else :
ts = uvfits.IntegrationTimeSeries(fin)[vismask] + toffset
res=numpy.zeros( ( sum(vismask) , 2))

if dc == "x" :
res[:,0] += ts * ws_px
elif dc == "y" :
res[:,1] += ts * ws_px
else:
raise "unkown direction"

return res


def InterpolateScreen(screenname,
pl):
"""Bi-linearly interpolate the screen"""

d=pyfits.open(screenname)[0].data
nx,ny=d.shape

res= numpy.zeros( len(pl) )

for i,p in enumerate(pl):

xl = int(math.floor(p[0]))
yl = int(math.floor(p[1]))

xh = xl+1
yh = yl+1

if xl == p[0] and yl ==p[1]:
# Excatly on pixel
res[i]=d[xl, yl]
else:
xd = p[0]-xl
yd = p[1]-yl

res[i] = ( xd * yd * d[xh,yh] +
(1-xd) * yd * d[xl, yh] +
xd * (1-yd) * d[xh, yl] +
(1-xd) * (1-yd) * d[xl, yl])

return res

def InterpolateScreen_fast(screenname,
pl):

"""Like InterpolateScreen, but faster

Note that this uses the B. Nikolic's astromap/pyplot library

:param pl: is the position list in *pixels*

"""

import pyplot

screen=pyplot.FitsMapLoad(screenname ,
1)

res= numpy.zeros( len(pl) )

for i,p in enumerate(pl):
bc=pyplot.MkBiLinearCoeffs(p[0],p[1], screen)
# Whoops, can loose pointers here
bc.this.acquire()
res[i]= bc(screen)


return res

def InterpolateScreen_ff(screenname,
pl):

"""Like InterpolateScreen, but faster

This useses both B.Nikolic's astromap and bnlib. It is the fastest
of the low however.
"""

# Make the imports here to avoid global dependency.
import pyplot
import pybnlib

screen=pyplot.FitsMapLoad(screenname ,
1)

n = len(pl)
xv = numarray.array( [xx[0] for xx in pl], numarray.Float64)
yv = numarray.array( [xx[1] for xx in pl], numarray.Float64)
res= numarray.zeros( n, numarray.Float64 )

pyplot.BiLinearMapInterp( pybnlib.doubleCvt( int(xv.__array_data__[0], 16) ),
pybnlib.doubleCvt( int(yv.__array_data__[0], 16) ),
pybnlib.doubleCvt( int(res.__array_data__[0], 16) ),
n,
screen);


return res


def BlowWindV(fin,
windspeed,
pixelscale,
toffset,
dc,
vismask=None):

"""Return wind vector suitable for simple blowing screens"""

return GenWindVector(fin, float(windspeed) / float(pixelscale),
toffset=toffset,
dc=dc,
vismask=vismask)

def WVToAntennaPosV(fin,
wv,
pixelscale,
vismask=None,
offset_px=[0,0]):
"""Convert a wind vector to antenna positions vector

Use offset_px to move away a little from the edgesp (in pixels)
"""

apos1, apos2= uvfits.BaselinePosList(fin)
if not vismask is None:
apos1=apos1[numpy.array(vismask)]
apos2=apos2[numpy.array(vismask)]

shift=numpy.zeros( (1,2))
for i in range(2):
shift[0,i] = min( [apos1[:,i].min() , apos2[:,i].min() ])
shift[0,i] -= ( offset_px[i] * pixelscale)


ap1_px= (apos1-shift)/pixelscale + wv
ap2_px= (apos2-shift)/pixelscale + wv

return ap1_px, ap2_px

def BlowSimple(fnamein,
screenname,
pixelscale,
windspeed,
fnameout,
phasescale=1.0,
toffset=0.0,
dc="x"):

"""Blow a phase screen across UV data.

:param fnamein: Input UV data
:param screenname: The phase screen as FITS file
:param scale: Scale of each pixel in the phase screen
:param windspeed: Rate at which screen moves (meters/second)
:param fnameout: Output file name
"""

fin=pyfits.open(fnamein)

wv= BlowWindV(fin,
windspeed, pixelscale,
toffset, dc)

ap1_px, ap2_px = WVToAntennaPosV(fin,
wv,
pixelscale)

print "Interpolating"
p1 = InterpolateScreen_ff( screenname, ap1_px)
p2 = InterpolateScreen_ff( screenname, ap2_px)

print "Expanding for IFs, windows, etc"
#scalings for various windows ifs etc....
pdiff = (p1-p2)* phasescale
for l in fin[0].data.field("data").shape[1:-1]:
pdiff= numpy.array( [pdiff for x in range(l) ] )
pdiff=pdiff.transpose()

print "Dephasing"
uvfits.DePhaseData(fnamein,
fnameout,
pdiff)


def BlowMultiSource(fnamein,
screendict,
pixelscale,
windspeed,
fnameout,
phasescale=1.0):

"""Blow phase screens, allow different screens for different sources

screendict is a dictionary connecting source numbers with screen
names.

e.g.:

BlowMultiSource("temp/mid-9-10-12w_test.ms.fits",
{ 1 :"temp/test_screen_9.fits"},
24.0 , 12,
"temp/testcorrupt.fits")
"""

fin=pyfits.open(fnamein)

pdiff=numpy.zeros( uvfits.NVis(fin) )

# For each source separatelly
for srcno in screendict.keys():
screenname=screendict[srcno]
mask = uvfits.SourceVisMask(fin, srcno)

wv= BlowWindV(fin,
windspeed, pixelscale,
0.0,
"x",
vismask=mask)

ap1_px, ap2_px = WVToAntennaPosV(fin,
wv,
pixelscale,
vismask=mask,
offset_px=[2,2])

p1 = InterpolateScreen_ff( screenname, ap1_px)
p2 = InterpolateScreen_ff( screenname, ap2_px)

pdiff[numpy.array(mask)] = (p1-p2)* phasescale

# Expand across IFs, etc
for l in fin[0].data.field("data").shape[1:-1]:
pdiff= numpy.array( [pdiff for x in range(l) ] )
pdiff=pdiff.transpose()

uvfits.DePhaseData(fnamein,
fnameout,
pdiff)



def BaselineRMS(fnamein,
screenname,
pixelscale,
windspeed,
baseline,
phasescale=1.0,
**kwargs):

"Compute simple baseline RMS for reference to imaging sims"

wv=BlowWindV(pyfits.open(fnamein),
windspeed, pixelscale,
0, "x")

sin = pyfits.open(screenname)

hbase= numpy.array( [ 10, sin[0].data.shape[0]/2 + baseline/pixelscale/2] )
lbase= numpy.array( [ 10, sin[0].data.shape[0]/2 - baseline/pixelscale/2] )

ap1_px= hbase + wv
ap2_px= lbase + wv

p1 = InterpolateScreen_ff( screenname, ap1_px)
p2 = InterpolateScreen_ff( screenname, ap2_px)

#scalings for various windows ifs etc....
pdiff = (p1-p2)* phasescale

return numpy.std(pdiff)