git.gag.com Git - debian/gnuradio/blob - gr-radio-astronomy/src/python/local_calibrator.py

   1 #!/usr/bin/env python
   2 #
   3 # Copyright 2003,2004,2005 Free Software Foundation, Inc.
   4 #
   5 # This file is part of GNU Radio
   6 #
   7 # GNU Radio is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation; either version 2, or (at your option)
  10 # any later version.
  11 #
  12 # GNU Radio is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with GNU Radio; see the file COPYING.  If not, write to
  19 # the Free Software Foundation, Inc., 51 Franklin Street,
  20 # Boston, MA 02110-1301, USA.
  21 #
  22
  23 import Numeric
  24 import math
  25 import ephem
  26 import time
  27
  28 #
  29 # Simple class for allowing local definition of a calibration function
  30 #  for raw samples coming from the RA detector chain.  Each observatory
  31 #  is different, and rather than hacking up the main code in usrp_ra_receiver
  32 #  we define the appropriate function here.
  33 #
  34 # For example, one could calibrate the output in Janskys, rather than
  35 #  dB.
  36 #
  37 #
  38
  39 def calib_default_total_power(data):
  40     r = 10.0*math.log10(data)
  41     return(r)
  42
  43 def calib_numogate_ridge_observatory_total_power(data):
  44
  45     me = ephem.Observer()
  46
  47     #
  48     # PyEphem wants lat/long as strings, rather than floats--took me quite
  49     #  a long time to figure that out.  If they don't arrive as strings,
  50     #  the calculations for sidereal time are complete garbage
  51     #
  52     me.long = str(-76.043)
  53     me.lat = str(44.967)
  54
  55     me.date = ephem.now()
  56     sidtime = me.sidereal_time()
  57
  58     foo = time.localtime()
  59     if not "calib_prefix" in globals():
  60         pfx = "./"
  61     else:
  62         pfx = globals()["calib_prefix"]
  63     filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
  64        foo.tm_mon, foo.tm_mday, foo.tm_hour)
  65
  66     numogate_file = open (filenamestr+".tpdat","a")
  67
  68     r = (data / 409.6)
  69     flt = "%6.3f" % r
  70     #r = calib_default_total_power(data)
  71     inter = globals()["calib_decln"]
  72     integ = globals()["calib_integ_setting"]
  73     fc = globals()["calib_freq_setting"]
  74     fc = fc / 1000000
  75     bw = globals()["calib_bw_setting"]
  76     bw = bw / 1000000
  77     ga = globals()["calib_gain_setting"]
  78
  79     now = time.time()
  80
  81     if not "calib_then_tpdat" in globals():
  82         globals()["calib_then_tpdat"] = now
  83
  84     if (now - globals()["calib_then_tpdat"]) >= 20:
  85         globals()["calib_then_tpdat"] = now
  86
  87         numogate_file.write(str(ephem.hours(sidtime))+" "+flt+" Dn="+str(inter)+",")
  88         numogate_file.write("Ti="+str(integ)+",Fc="+str(fc)+",Bw="+str(bw))
  89         numogate_file.write(",Ga="+str(ga)+"\n")
  90     else:
  91         numogate_file.write(str(ephem.hours(sidtime))+" "+flt+"\n")
  92
  93     numogate_file.close()
  94     return(r)
  95
  96 def calib_numogate_ridge_observatory_fft(data,l):
  97
  98     me = ephem.Observer()
  99
 100     #
 101     # PyEphem wants lat/long as strings, rather than floats--took me quite
 102     #  a long time to figure that out.  If they don't arrive as strings,
 103     #  the calculations for sidereal time are complete garbage
 104     #
 105     me.long = str(-76.043)
 106     me.lat = str(44.967)
 107
 108     me.date = ephem.now()
 109     sidtime = me.sidereal_time()
 110
 111     foo = time.localtime()
 112
 113     if not "calib_prefix" in globals():
 114         pfx = "./"
 115     else:
 116         pfx = globals()["calib_prefix"]
 117     filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
 118        foo.tm_mon, foo.tm_mday, foo.tm_hour)
 119
 120     now = time.time()
 121
 122     if not "calib_then" in globals():
 123         globals()["calib_then"] = now
 124
 125     delta = (l/1024)*5
 126
 127     if (now - globals()["calib_then"]) >= delta:
 128
 129         globals()["calib_then"] = now
 130         numogate_file = open (filenamestr+".sdat","a")
 131
 132         r = calib_default_fft(data,l)
 133         inter = globals()["calib_decln"]
 134         fc = globals()["calib_freq_setting"]
 135         fc = fc / 1000000
 136         bw = globals()["calib_bw_setting"]
 137         bw = bw / 1000000
 138         av = globals()["calib_avg_alpha"]
 139         numogate_file.write("data:"+str(ephem.hours(sidtime))+" Dn="+str(inter)+",Fc="+str(fc)+",Bw="+str(bw)+",Av="+str(av))
 140         numogate_file.write(" "+str(r)+"\n")
 141         numogate_file.close()
 142         return(r)
 143
 144     return(data)
 145
 146 def calib_default_fft(db,l):
 147     return(db)
 148
 149 #
 150 # We capture various parameters from the receive chain here, because
 151 #  they can affect the calibration equations.
 152 #
 153 #
 154 def calib_set_gain(gain):
 155     globals()["calib_gain_setting"] = gain
 156     globals()["calib_then_tpdat"] = time.time() - 50
 157
 158 def calib_set_integ(integ):
 159     globals()["calib_integ_setting"] = integ
 160     globals()["calib_then_tpdat"] = time.time() - 50
 161
 162 def calib_set_bw(bw):
 163     globals()["calib_bw_setting"] = bw
 164     globals()["calib_then_tpdat"] = time.time() - 50
 165
 166 def calib_set_freq(freq):
 167     globals()["calib_freq_setting"] = freq
 168     globals()["calib_then_tpdat"] = time.time() - 50
 169
 170 def calib_set_avg_alpha(alpha):
 171     globals()["calib_avg_alpha"] = alpha
 172
 173 def calib_set_interesting(inter):
 174     globals()["calib_is_interesting"] = inter
 175
 176 def calib_set_decln(dec):
 177     globals()["calib_decln"] = dec
 178     globals()["calib_then_tpdat"] = time.time() - 50
 179
 180 def calib_set_prefix(pfx):
 181     globals()["calib_prefix"] = pfx