Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.mrao.cam.ac.uk/~bn204/alma/memo-infer-587/air.py Дата изменения: Fri Aug 10 20:17:35 2012 Дата индексирования: Tue Oct 2 10:25:50 2012 Кодировка: Поисковые слова: m-a 17

"""
Bojan Nikolic ,
Initial version December 2008
Revised 2009

The driver module for all libAIR functionality and retrieval of
atmospheric parameters.

Do not put into this file:
- Testing routines
- Input/Output or plotting
"""
import math

import numpy

import setup

setup.pyair_minV("0.8")
import pyair

setup.pybnmin1_minV("1.4")
import pybnmin1

from bnmin1utils import SetIC, ChainToArray

MK_weak_pr=( ("n", 0, 10),
("T", 200, 320),
("P", 100, 650))

MK_flat_pr=( ("n", 0, 5),
("T", 260, 280),
("P", 530, 610))

MK_tight_pr=( ("n", 0, 5),
("T", 260, 280),
("P", 570, 590))

def applyPriors(obsmod,
prior=None,
log_pr=False):
"""
Apply prior information to a model -- supports flat or log-flat
priors

:param obsmod: The base model

:param prior: List of lists defining the prior ranges

:param log_pr: If true, use log-priors

"""
if prior is None:
return obsmod
else:
obsmod.thisown=False
if log_pr :
po=pybnmin1.LogFlatPriors(obsmod)
else:
po=pybnmin1.IndependentFlatPriors(obsmod)
for pname, plow, phigh in prior:
po.AddPrior(pname, plow, phigh)
return po

def mkModel(coupling=None,
TSpill=275,
za=0):
"""
Make a libAIR model based on user inputs

:param za: Zenith angle for the plane-parrelel model
"""
wm=pyair.SingleLayerWater(pyair.ALMADickeProto,
pyair.PartTable,
pyair.AirCont)
wm.thisown=False
if coupling is not None:
wm=pyair.CouplingModel(wm)
wm.setSpill(coupling,
TSpill)
wm.thisown=False
wm=pyair.PPDipModel(wm)
wm.setZA(za)
return wm

def AbsRetrieve(err=1.0,
d=None,
ic=[1.0, 270, 600],
nsample=2000000,
prior=None,
coupling=None,
TSpill=275,
el=None,
PSigma=0.1,
TSigma=0.1):
"""
Absolute retrieval from brightness temperatures only. For testing
& comparison to the joint analysis.

:param err: Gaussian error on each measurement

:param d: The observed data. If not supplied (is None) then the
model for the initial condition point (ic) is used. The list is in
[ch1, ch2, ch3, ch4] format.

:param ic: The initial condition to start the chain from

:param nsample: Number of proposals to make

:param coupling: The coupling of the WVR to the sky

:param el: Elevation in radian. If None assume pointing straight
up

"""
if el is None:
za=0
else:
za=math.pi*0.5-el
wm=mkModel(coupling=coupling,
TSpill=TSpill,
za=za)
ll=pyair.AbsNormMeasure(wm)
ll.thermNoise=pybnmin1.DoubleVector([err]*4)
SetIC(ll, ic)
if d is None:
ll.modelObs()
else:
ll.obs= pybnmin1.DoubleVector(d)
print "Simulated points is ", list(ll.obs)
ll=applyPriors(ll,
prior=prior,
log_pr=True)
mm=pybnmin1.MetropolisMCMC(ll,
[0.001, TSigma, PSigma],
33)
if coupling is not None:
mm.getbyname("coupling").dofit=0

return ChainToArray(mm.sample(nsample))

def dipRetrieve(elv,
tskyv,
err=1.0,
ic=[1.0, 270, 600, 0.99, 280],
nsample=200000,
fitTTerm=False
):
"""
Retrieve model parameters from a sky-dip measurement

:param err: Random error on each brighntess measured, in Kelvin,
assumed independent and gaussian
"""
wm=mkModel(coupling=ic[3],
TSpill=ic[4],
za=0)
ll=pyair.DipNormMeasure(wm)
SetIC(ll, ic)
for el,tsky in zip(elv,tskyv):
ll.addObs(math.pi/2 - el,
pyair.doubleV(tsky))
ll.thermNoise=pybnmin1.DoubleVector([err]*4)
sigmas=[0.001, 0.1, 0.1, 0.01]
if fitTTerm:
sigmas += [0.1]
mm=pybnmin1.MetropolisMCMC(ll,
sigmas,
33)
if fitTTerm:
mm.getbyname("TTerm").dofit=1
return ChainToArray(mm.sample(nsample))






def calcdTdL(pchain,
coupling=None,
TSpill=275,
el=math.pi*0.5):
"""
Calculate dTdL from a chain of model parameters

:param pchain: Chain of model parameter values

:returns: Chain of the four dTdL coefficients
"""
wm=mkModel(coupling=coupling,
TSpill=TSpill,
za=math.pi*0.5-el)
res=[]
for p in pchain:
SetIC(wm,
p)
v=pybnmin1.DoubleVector([0,0,0,0])
wm.dTdL_ND(v)
res.append(list(v))
return numpy.array(res)