Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.mrao.cam.ac.uk/~eic22/ALMAmemos/memo590/dispersivephase.py Дата изменения: Wed Nov 11 16:36:38 2009 Дата индексирования: Tue Oct 2 12:40:04 2012 Кодировка:

"""
Emily Curtis

Script for working with ATM dispersive phase for Memo 590
"""

import os
import string
import numpy
from atmtools import *

"""System variable parameters:"""
class uservar:
outputdir = 'Newdata' #Users' output directory
dcval = 0.01 #Change in water vapour to compute change in path /mm
fmin = 0 #Minimum frequency to compute phase /GHz
fmax = 1000 #Maximum frequency to compute phase /GHz
fstep = 1 #Frequency step /GHz

"""Define median atmosphere parameters"""
class medatmos:
pwv = 1.22 #Median precipitable water vapour /mm
pg = 560 #Median ground pressure /mb
Tg = 270 #Median ground temperature/K
tlr = -7.28 #Median tropospheric lapse rate/K/km
hzero = 1.16 #Water vapour scale height/km

findATMPath()

_atmexecpath=findATMPath()

def keepcolumns(x,col1,col2):
"""Function to keep columns col1 and col2 of numpy array"""
y = numpy.zeros((x.shape[0],x.shape[1]-2),float)
y[:,0]=x[:,col1]
y[:,1]=x[:,col2]
return y

def subcolumns(x,y,a):
"""Function to subtract the 2nd column of two matrices x y, first column is the same and multiply each by a"""
z = numpy.zeros(x.shape,float)
z[:,0]=x[:,0]
z[:,1]=numpy.subtract(y[:,1]*a,x[:,1]*a)
return z

def dividecolumns(x,ref):
z = numpy.zeros((x.shape[0],3),float)
z[:,0] = x[:,0]
z[:,1] = x[:,1]
z[:,2] = numpy.divide(x[:,1],ref[:,1])
return z

def finddispcoefficients(f,deltac,params):
"""Function to calculate the dispersive phase correction coefficients """
print 'Frequency range... ',f[0],':',f[1],' GHz in steps of',f[2],' GHz'
print 'Model parameters...'
print 'pwv=',params[0],'mm pg=',params[1],'mb Tg=',params[2],'K TLR=',params[3],'K/km h0=',params[4],'km'
r1 = keepcolumns(parsefgrid(runatm("dispersive", {"fmin":f[0], "fmax":f[1], "fstep":f[2], "pwv":params[0], "gpress":params[1], "gtemp":params[2],"tlr":params[3], "wvl":params[4]})),0,2)
print ''
print 'Calculating numerical derivative, offset water vapour:',deltac,'mm'

r2 = keepcolumns(parsefgrid(runatm("dispersive", {"fmin":f[0], "fmax":f[1], "fstep":f[2], "pwv":params[0]+deltac, "gpress":params[1], "gtemp":params[2],"tlr":params[3], "wvl":params[4]})),0,2)

z = subcolumns(r1,r2,1000.0/deltac)

return z

def textout(filename,array):
file = open(filename, 'w')
for line in array:
dataline = string.join(map(str, line))
dataline = dataline+'\n'
file.write(dataline)
file.close()

cval = [0.44, 0.69, 2.56, 5.45] # Values of pwv to compute dL/dc
pval = [520,540,580,600]
Tval = [262,265,276,281]
tlrval = [-9.71,-8.83,-5.69,-4.80]
hval = [0.97,1.06,1.29,1.54]


"""Compute median atmosphere parameters"""

print '==============================================================='
print 'Computing phase for median atmosphere'
print '==============================================================='
refcoeffs = finddispcoefficients([uservar.fmin,uservar.fmax,uservar.fstep],uservar.dcval,[medatmos.pwv,medatmos.pg,medatmos.Tg,medatmos.tlr,medatmos.hzero])
textout(uservar.outputdir+'/pwv/dldc_median.dat',refcoeffs)


print ''
print '==============================================================='
print 'Computing phase for variations in c'
print '==============================================================='

for c in cval:
newcoeffs = finddispcoefficients([uservar.fmin,uservar.fmax,uservar.fstep],uservar.dcval,[c,medatmos.pg,medatmos.Tg,medatmos.tlr,medatmos.hzero])
division = dividecolumns(newcoeffs,refcoeffs)
textout(uservar.outputdir+'/pwv/dldc_'+str(int(c*1000))+'um.dat',division)

print ''
print '==============================================================='
print 'Computing phase for variations in T'
print '==============================================================='

for temp in Tval:
newcoeffs = finddispcoefficients([uservar.fmin,uservar.fmax,uservar.fstep],uservar.dcval,[medatmos.pwv,medatmos.pg,temp,medatmos.tlr,medatmos.hzero])
division = dividecolumns(newcoeffs,refcoeffs)
textout(uservar.outputdir+'/T/dldc_'+str(temp)+'K.dat',division)

print ''
print '==============================================================='
print 'Computing phase for variations in ground pressure'
print '==============================================================='

for p in pval:
newcoeffs = finddispcoefficients([uservar.fmin,uservar.fmax,uservar.fstep],uservar.dcval,[medatmos.pwv,p,medatmos.Tg,medatmos.tlr,medatmos.hzero])
division = dividecolumns(newcoeffs,refcoeffs)
textout(uservar.outputdir+'/p/dldc_'+str(p)+'mbar.dat',division)

print ''
print '==============================================================='
print 'Computing phase for variations in lapse rate'
print '==============================================================='

for tlapse in tlrval:
newcoeffs = finddispcoefficients([uservar.fmin,uservar.fmax,uservar.fstep],uservar.dcval,[medatmos.pwv,medatmos.pg,medatmos.Tg,tlapse,medatmos.hzero])
division = dividecolumns(newcoeffs,refcoeffs)
textout(uservar.outputdir+'/tlr/dldc_'+str(tlapse)+'kkms.dat',division)

print ''
print '==============================================================='
print 'Computing phase for variations in water vapour scale height'
print '==============================================================='

for scaleheight in hval:
newcoeffs = finddispcoefficients([uservar.fmin,uservar.fmax,uservar.fstep],uservar.dcval,[medatmos.pwv,medatmos.pg,medatmos.Tg,medatmos.tlr,scaleheight])
division = dividecolumns(newcoeffs,refcoeffs)
textout(uservar.outputdir+'/h/dldc_'+str(scaleheight)+'km.dat',division)


"""
the following compute the data for Figs. 9 & 10

"""
tcons = [271,269,272,268]
for T in tcons:
newcoeffs = finddispcoefficients([uservar.fmin,uservar.fmax,uservar.fstep],uservar.dcval,[medatmos.pwv,medatmos.pg,T,medatmos.tlr,medatmos.hzero])
division = dividecolumns(newcoeffs,refcoeffs)
textout(uservar.outputdir+'/tempconstraints/dldc_'+str(T)+'K.dat',division)

tlrcons = [-8.08,-6.48,-8.78,-5.78]

for tlapse in tlrcons:
newcoeffs = finddispcoefficients([uservar.fmin,uservar.fmax,uservar.fstep],uservar.dcval,[medatmos.pwv,medatmos.pg,medatmos.Tg,tlapse,medatmos.hzero])
division = dividecolumns(newcoeffs,refcoeffs)
textout(uservar.outputdir+'/tempconstraints/dldc_'+str(tlapse*-1)+'Kkm.dat',division)