# To calculate the time averaging interval dt to achieve time-smearing
# at radius theta of no worse than R where R = apparent/true amplitude
# of an unresolved source, for synthesised beamwidth thetaB
#
# thetaB can be approximated as (wavelength/longest baseline) but for
# short tracks or low-elevation sources, the natural, synthesised
# beam-width should be used. This expression becomes less accurate for
# non-circular uv tracks; using the minor beam axis will give a
# conservative estimate of dt.
#
# The phase-rate dphi/dt may also cause smearing at large theta.
#

#
R = 0.95               # lowest acceptable dimunition of signal (Iobs/Ireal)
thetaB = 1.56          # synthesised beam size (arcsec)
theta = 31.            #  outermost source, arcsec from centre

###############################################################
#
# This form is useful for predictions; it could also be modified to
# run inside CASA on a given MS to extract frequency and other
# parameters and apply smoothing automatically.
#
#################################################################
# Calculate time smearing
#
# Other effects e.g. primary beam, phase rate, may provide worse
# distortions at theta.
# Confusing as well as wanted sources should be considered
#
# Based on Bridle & Schwab (B&S) ch 18 in Synthesis Imaging II
# (Taylor, Carilli & Perley 1999)
# Average intensity loss for a circumpolar point source is given by
# R = 1-const (theta/thetaB)^2 dt^2
#
# where const = 1.08E-09 for uniform, circular uv coverage (possibly
# suitable for the most compact configurations) and
# C = 1.22E-09 for a Gaussian distribution of uv coverage (more
# suitable for most, especially extended configurations).
# So
# dt = sqrt[(1-R)/const] (thetaB/theta) 
######################################################
from glob import glob
from matplotlib import *
from pylab import *
import math
import matplotlib.pyplot as plt
import numpy as np
import re
C = [1.08E-09, 1.22E-09]

dtU = (math.sqrt((1.-R)/C[0]))*thetaB/theta
dtG = (math.sqrt((1.-R)/C[1]))*thetaB/theta

# Phase rate df = dphi/dt = 2 pi (theta/thetaB)/24/3600 (rad/sec)
df = 2.*math.pi*(theta/thetaB)/86400.
# leading to an amplitude reduction of Rf = sinc(df*dtG/2.) 
# (Perley, ch 13 ibid, eq 13-15).
# But, (theta/thetaB) cancels, giving
Rf = (sin(1.041*math.sqrt(1-R)))/(1.041*math.sqrt(1-R))

# total phase change in dtG
tf = df*dtG

# time interval for phase change of pi/6
dtF = math.pi/6./df

t = min(dtG, dtF) # tightest constraint on averaging time

# Dynamic range in 1 hr of t-sec integrations, 50 antennas, pi/6 phase
# change. This assumes that each integration is independent which may
# not be true for small t

DdtF50 = sqrt(3600./t)*50./(math.pi/6.)
DdtF16 = sqrt(3600./t)*16./(math.pi/6.)  # 16 antennas

# time for dynamic range of 1000
t1000 = (3600.*50.**2./(df*1000.)**2)**(1./3.)

dtU = '%3.f' % (dtU)
dtG = '%3.f' % (dtG)

df = '%4.2f' % (df*180./math.pi)
Rf =  '%4.2f' % (Rf)
tf = '%3.f' % (tf*180./math.pi)
dtF = '%3.f' % (dtF)
t = '%3.f' % (t)
t1000 = '%3.f' % (t1000)
DdtF50 = '%4.f' % (DdtF50)
DdtF16 = '%4.f' % (DdtF16)

print '\nSources at a radius of ',theta, 'arcsec, observed using a synthesised beam width \nof', thetaB, 'arcsec will show an an apparent to actual flux density reduction \nratio ~', R,'due to time-smearing in averaged integration time',dtU,'to',dtG,'sec\nfor very well-filled to more extended array configurations.\n'
print 'The phase rate at',theta, 'arcsec is ',df,' deg/s, leading to a reduction of \namplitude to',Rf,'in',dtG,'s. To avoid calibration and imaging errors a \nphase change < 30 deg and so an averaging time of <',dtF,'s is required. \nThe shortest of these averaging times, ',t,'s limits the dynamic range to \n<',DdtF50,'in 1 hr on-source, 50 antennas, or  <',DdtF16,' with 16 antennas. \nAn averaging time <',t1000,'s is needed for dynamic range 1000, 50 antennas.\nThis may be further reduced by noise and calibration errors.' 

###############################################################
#
# Time smearing at primary beam FWHM ES compact config ~1190 s
dtU = (math.sqrt((1.-R)/C[0]))*0.175
# Time smearing at primary beam FWHM ES extended config ~330 s
dtG = (math.sqrt((1.-R)/C[1]))*0.051

#############################################################
# Plot time smearing interval against baseline length
# for imaging full field of view
#
# and limit for phase-rate smearing using
# dphi/dt = 2 pi (theta/thetaB)/24./3600.
clf()

PB = 12.    # Antenna diameter
baselines=np.arange(100., 16000.,10.)
bc = 0.85   # fudge factor to make primary beam/baseline ratio
            # consistent with FoV/resolution in documentation

#dtU = (math.sqrt((1.-R)/C[0]))*2.*bc*PB/baselines

a1=axes([0.1,0.1,0.7,0.7])
R=0.95
dtG = (math.sqrt((1.-R)/C[1]))*2.*bc*PB/baselines
a1.plot(baselines, dtG,'-',color='b',linewidth=2)
R=0.8
dtG = (math.sqrt((1.-R)/C[1]))*2.*bc*PB/baselines
a1.plot(baselines, dtG,'--',color='b',linewidth=2)

# Phase rate at primary beam FWHM
dfPB = 2.*math.pi*(baselines/(bc*PB))/86400.
# max averaging time to reach a dynamic range of 1000 in 1 hr, 50 ants.
D=100.
Nant=16
tPB = (3600.*Nant**2./(dfPB*D)**2.)**(1./3.)
a1.plot(baselines, tPB,'--',color='r',linewidth=2)
D=1000.
Nant=50
tPB = (3600.*Nant**2./(dfPB*D)**2.)**(1./3.)
a1.plot(baselines, tPB,'-',color='r',linewidth=2)


ylabel('Averaged integration time (s)')
xlabel('Baseline length (m)')

a1.set_xscale('log')
a1.set_yscale('log')
title('Time smearing and dynamic range limit at primary beam FWHM',size=11)
a1.xaxis.grid(True, which='majorminor')
a1.yaxis.grid(True, which='majorminor')

a1.plot([100,170],[9,9],'-',color='r',linewidth=2)
a1.plot([100,170],[11,11],'--',color='r',linewidth=2)

a1.plot([100,170],[14,14],'-',color='b',linewidth=2)
a1.plot([100,170],[17,17],'--',color='b',linewidth=2)

text(200,8,'Dynamic range 1000, 50 antennas',color='r')
text(200,10,'Dynamic range 100, 16 antennas',color='r')

text(200,21,'Time smearing reduction',color='b')

text(200,13,'95% of peak',color='b')
text(200,16,'80% of peak',color='b')


a1.set_xlim(80.,19000.)
a1.set_ylim(7,600.)

#show()
#savefig('ALMA_timesmear.eps')
#savefig('ALMA_timesmear.png')