# De-convolution of a Gaussian beam from planetary observationsΒΆ

From `astromap` version 1.5b3, the OOF software supports deconvolution of an observation of a planet to retrieve the assumed Gaussian parameters of the beam. This deconvolution can be useful for empirical monitoring of telescope beams over time over which the planets used as sources vary significantly in size.

Here is a simple script showing how to carry out this deconvolution:

```# Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>
# Initial version 2005.
#
"""
Example of how to deconvolve for a Gaussian beam, assuming underlying
source is a planet represetned by a tapered disk
"""

import pickle
import numpy

import localsetup

import pybnlib
import pyplot
import pybnmin1
import bnmin1utils

def setPlanet(m,
Cm=0):
"""
Set the map to with a planet-like shape
"""
pfn=pybnlib.TaperedTopHatDD()
pfn.Cm = Cm
pyplot.WorldSet(m,
pfn)

def mkPlanet(npix,
mapsize,
Cm=0):
"""
Create a map with planet model

:param npix: Size of map in pixels per dimension

:param mapsize: Angular size to assign to each map dimension

Example:

>>> m=mkPlanet(256, 20.0, 500.0)

"""
m=pyplot.Map(npix, npix)
pyplot.MkRectCS(m,
mapsize*0.5,
mapsize*0.5)
setPlanet(m,
Cm=Cm)
return m

def mkMockObs(npix,
mapsize,
noise,
params=[1, 0,  0,  1.0, 0, 0]):
"""
Create a mock observation

:param parms: Parameters for the Gaussian. The order is amplitude,
x0, y0, sigma, rho and diff
"""
o=mkPlanet(npix,
mapsize)
# First parameter is not used since we never fit
fm=pyplot.GaussConvMap(o,
o)
bnmin1utils.SetIC(fm,
params)
o.mult(0)
fm.eval(o)
mnoise=pyplot.Map(npix, npix)
pyplot.NormDist(mnoise,
1.0*noise)
return o

def fitObs(o,
mapsize,
ic=[1, 0,  0,  1.0, 0, 0],
fit=["amp", "x0", "y0", "sigma", "diff", "rho" ],
Cm=0):
"""
Fit the observations.

:param mapsize: Angular size of the map

:param ic: Initial condition from which to start the minimiser

:param fit: List of parameters to fit for -- default is all
parameters

:returns: List of [best fitting model (*after convolution*),
best fitting Gaussian parameters, final chi-sq]

Example:
>>> m, r, ch=fitObs(obsmap, 1.0, 50.0, ic=[0, 0, 0, 1.0, 0, 0])

"""
npix=o.nx
pm=mkPlanet(npix,
mapsize,
Cm=Cm)
fm=pyplot.GaussConvMap(o,
pm)
lmm=pybnmin1.LMMin(fm)
for p in fit:
lmm.getbyname(p).dofit=1
m1=pybnmin1.ChiSqMonitor()
m1.thisown=0
lmm.solve()
chisq=lmm.ChiSquared()
bestmodel=pyplot.Map(npix,
npix)
pyplot.MkRectCS(bestmodel,
mapsize*0.5,
mapsize*0.5)
fm.eval(bestmodel)
return bestmodel, [lmm.getbyname(p).getp() for p in  fit], chisq

def npToMap(a,
mapsize):
"""
Conver a numpy array to a map.

:param mapsize: Angluar size of one of the dimensions of the map.

"""
m=pyplot.Map(*a.shape)
for i in range(a.shape[0]):
for j in range(a.shape[1]):
m.set(i, j, a[i,j])
pyplot.MkRectCS(m,
mapsize*0.5,
mapsize*0.5)
return m
```

#### Previous topic

Building the OOF package: worked solutions

#### Next topic

Aperture-plane phase change due to defocus of a Gregorian telescope