Friday, March 22, 2013

Atmospheric transmission curve for Plateau de Bure Interferometer

I tried to find the atmospheric transmission curve for Plateau de Bure Interferometer (PdBI) site the other day. There is no plot/data file showing the atmospheric transmission for the Plateaude Bure Interferometer site! Then I realized that I might be able to plot it in the GILDAS software ASTRO - the observation planning tasks of the bundle. So in this post I will explain how to plot the atmospheric transmission curve at the Plateau de Bure Interferometer (PdBI) using the atmospheric modeling function included in the ASTRO package included in the IRAM GILDAS software bundle, and Python. I guess you can do the same plot in ASTRO using GREG, but I like Python/matplotlib more.

After checking and sending mails back and forth with Jan Martin Winters at IRAM, I realized that there is not direct command to plot the whole transmission curve. What does exist is the command ATMOSPHERE, that calculates τ (the opacity) and a whole bunch of other stuff for a given frequency (running "HELP ATM" gives you the run through of the input/output). This means that I can write a loop that calculates τ for a given frequency. Then I just have to calculate the transmission from I/I0=e-τ and plot it against the frequency.

After each calculation I want to print the frequency and the τ to a file. Thus, the following script produces a file where the first column is frequency and the second the total opacity (τ).

! use the new atm, the 2009 version 
! test that it works with "exam atm%version" 
! after the command below
set atm new
! set observatory, the available names are
! if not satisfied, use coordinates as
! OBS [longitude latitude altitude [radius]]
! radius is the sun avoidance zone in degrees, (def 30) 
obs pdbi                      !<---- change site here!
! which airmass to calculate for 1 = zenith
let airmass 1
! ground temperature in Kelvin
let temperature 273
! create variable freq 
! and the incremental change in frequency
define real freq
let freq 0.1
define real incr
let incr 0.1
! the amount of water vapor in mm 
let water 0.3                ! <---- change water vapor here!
! the name of the output file
sic out pdbi_atm_0_3.dat
for n 0 to 16000
    let freq_sig = freq
    let freq_ima = freq
    ! write to file
    say 'freq' 'tau_tot'
    ! add incr to the frequency and 
    ! restart the loop
    let freq = freq+incr
sic out

Just run this with

astro -nw @script_name.astro

After this you are left with the file pdbi_atm_0_3.dat in this case. Run this script for different amount of water vapor, say 0.3, 0.8 and 1.5. Then you jump into Python and do the following in the same directory as your files

# general imports
from scipy import exp, loadtxt
import matplotlib.pyplot as pl

# various settings
pl.rcParams['font.size'] = 15

# calculate the transmission in percent
transm_percent = lambda tau : exp(-1*tau)*100

# load the data columns
data0_3 = loadtxt('pdbi_atm_0_3.dat').transpose()
data0_8 = loadtxt('pdbi_atm_0_8.dat').transpose()
data1_5 = loadtxt('pdbi_atm_1_5.dat').transpose()

# initialize the figure

# plot the graphs
pl.plot(data0_3[0], transm_percent(data0_3[1]))
pl.plot(data0_8[0], transm_percent(data0_8[1]))
pl.plot(data1_5[0], transm_percent(data1_5[1]))

# get some label action going
pl.legend(['0.3mm', '0.8mm', '1.5mm'])
pl.xlabel('Frequency [GHz]')
pl.ylabel('Transmission (%)')
pl.title('Atmospheric Transmission at PdBI')

# and adjust the margins

What you should end up with is the following

Notable is that it "only" goes up to 1 THz, after this the model does not work, or I could not get it to work at least. The values are not correct. If you figure it out, please let me know. Next I will add the receiver bands and some nice water lines, or something. I have to thank Susanne Wampfler for helping me out with the ASTRO script, programming in SIC is not straight forward...

No comments:

Post a Comment