! 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
! BURE (or PDBI), PICOVELETA (or VELETA) ALMA, APEX,
! ATF, CSO, FCRAO, JCMT, KITTPEAK, KOSMA, SEST, SMA,
! SMT EFFELSBERG, VLA LASILLA, MAUNAKEA, PARANAL (or VLT),
! IOTA, PTI, GI2T
! 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
atm
! write to file
say ‘freq’ ‘tau_tot’
! add incr to the frequency and
! restart the loop
let freq = freq+incr
next
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.ion()
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
pl.figure(1,(11,6))
# plot the graphs
pl.plot(data0_3[0], transm_percent(data0_3[1]))
pl.xlabel(‘Frequency [GHz]’)
pl.ylabel(‘Transmission (%)’)
pl.title(‘Atmospheric Transmission at PdBI’)
What you should end up with is the following
Leave a Reply