def enu2itrf(self, antennas=None, lon=None, lat=None):
        """ Converts a list of ENU positions to ITRF positions. 
            Requires a reference positions (lon,lat) in radians. 
           
            Returns the ITRF reference position and the ITRF list of positions.
        """
        
        dm  = pyrap.measures.measures()
        
        antennas = self.antennas if antennas is None else antennas
        lon = lon or self.lon
        lat = lat or self.lat
        
        # convtert reference position to itrf system
        refpos_wgs84 = dm.position('wgs84', dq.quantity(lon, 'rad'),
                       dq.quantity(lat, 'rad'))
        refpos = dm.measure(refpos_wgs84, 'itrf')

        lon,lat,rad = [ refpos[x]['value'] for x in 'm0 m1 m2'.split() ]

        xyz0 = rad*numpy.array( [math.cos(lat)*math.cos(lon),
                        math.cos(lat)*math.sin(lon), math.sin(lat)])
        
        # 3x3 transformation matrix. Each row is a normal vector, i.e the rows are (dE,dN,dU)
        xform = numpy.array([
            [-math.sin(lon), math.cos(lon), 0],
            [-math.cos(lon)*math.sin(lat), -math.sin(lon)*math.sin(lat), math.cos(lat)],
            [math.cos(lat)*math.cos(lon), math.cos(lat)*math.sin(lon), math.sin(lat)]
])
        antennas = numpy.array(antennas)
        xyz = xyz0[numpy.newaxis,:] + antennas.dot(xform)

        return xyz0,xyz

 def riseset(self, crd, ev="5deg"):
     a = self.rise(crd, ev)
     if isinstance(a, str):
         a = a.split()
         return { "rise": { "last": a[0], "utc": a[0] },
                  "set" : { "last": a[1], "utc": a[1] },
                  "solved": False }
     ofe = self.measure(self._framestack["epoch"], "utc")
     if not is_measure(ofe):
         ofe = self.epoch('utc', 'today')
     x = a.copy()
     for k in x.keys():
         x[k] = self.measure(self.epoch("last", dq.totime(a[k]),
                                        off=self.epoch("r_utc",
                                                       (dq.quantity(ofe["m0"]) + dq.quantity("0.5d")))),
                             "utc")
     return { "rise": { "last": self.epoch("last",
                                           dq.totime(a["rise"])),
                        "utc": x["rise"] },
              
              "set": { "last": self.epoch("last",
                                          dq.totime(a["set"])),
                        "utc": x["set"] },
              "solved" : True
              }

def get_time_range(start_time,end_time,timestep,time_in_sec,TIME_OFFSET=0):
    if HAS_PYRAP:
      try:
        start_time = qu.quantity(start_time).get_value()
        end_time = qu.quantity(end_time).get_value()
        print ('**** specified start and end time ', start_time, end_time)
        reference_time = start_time * 86400.0 - TIME_OFFSET
        st = reference_time - timestep 
        et = end_time * 86400.0 +  timestep + TIME_OFFSET
        #print "getting string",reference_time
        str_start_time =  obtain_observation_year_month_day_hms(reference_time)
        timerange= [st, et]
      except:
        print ('no time range given')
        print ('exiting')
        return -1,-1,-1
    elif HAS_EPHEM:
      if time_in_sec:
        dublin_start = start_time / 86400.0 -15019.5
        dublin_end = end_time / 86400.0 -15019.5
        start_time = ephem.julian_date(ephem.Date(dublin_start)) - 2400000.5
        end_time = ephem.julian_date(ephem.Date(dublin_end)) - 2400000.5
      else:
        start_time = ephem.julian_date(ephem.Date(start_time)) - 2400000.5
        end_time = ephem.julian_date(ephem.Date(end_time)) - 2400000.5
      print ('ephem start and end time ', start_time, end_time)
      reference_time = start_time * 86400.0 - TIME_OFFSET
      st = reference_time - timestep 
      et = end_time * 86400.0 + timestep + TIME_OFFSET
      str_start_time =  obtain_observation_year_month_day_hms(reference_time)
      timerange= [st, et]
    else:
      print ('unable to get time range so exiting!')
      return -1,-1,-1
    return timerange,str_start_time,reference_time

 def direction(self, rf='', v0='0..', v1='90..', off=False):
     loc = { 'type': 'direction' , 'refer':  rf}
     loc['m0'] = dq.quantity(v0)
     loc['m1'] = dq.quantity(v1)
     if is_measure(off):
         if not off['type'] == "direction":
             raise TypeError('Illegal offset type specified.')
         loc["offset"] = off
     return self.measure(loc, rf)

def print_facet_centre_list(facet_centres, num_facet_centres):
    if num_facet_centres != 0:
        print "REQUESTED FACET CENTRES:"
        for i, c in enumerate(facet_centres):
            print "\tFACET %d RA: %s DEC: %s" % (
                i,
                quantity(c[0], "arcsec").get("deg"),
                quantity(c[1], "arcsec").get("deg"),
            )

 def earthmagnetic(self, rf='', v0='0G', v1='0..', v2='90..', off=False):
     loc = { 'type': "earthmagnetic", 'refer': rf }
     loc['m0'] = dq.quantity(v0)
     loc['m1'] = dq.quantity(v1)
     loc['m2'] = dq.quantity(v2)
     if is_measure(off):
         if not off['type'] == "earthmagnetic":
             raise TypeError('Illegal offset type specified.')
         loc["offset"] = off
     return self.measure(loc, rf)

 def todoppler(self, rf, v0, rfq=False):
     if is_measure(rfq) and rfq['type'] == 'frequency':
         rfq = dq.from_dict(rfq['m0'])
     if is_measure(v0):
         if v0['type'] == 'radialvelocity':
             return self.todop(v0, dq.to_dict(dq.quantity(1.,'Hz')))
         elif v0['type'] == 'frequency' and  isinstance(rfq, dq._quanta.Quantity) \
                  and rfq.conforms(dq.quantity('Hz')):
             return self.todop(v0, dq.to_dict(rfq))
         else:
             raise TypeError('Illegal Doppler or rest frequency specified')
     else:
         raise TypeError('Illegal Frequency specified')

def parallactic_angles(times, antenna_positions, field_centre):
    """
    Computes parallactic angles per timestep for the given
    reference antenna position and field centre.

    Arguments:
        times: ndarray
            Array of unique times with shape (ntime,),
            obtained from TIME column of MS table
        antenna_positions: ndarray of shape (na, 3)
            Antenna positions, obtained from POSITION
            column of MS ANTENNA sub-table
        field_centre : ndarray of shape (2,)
            Field centre, should be obtained from MS PHASE_DIR

    Returns:
        An array of parallactic angles per time-step

    """
    import pyrap.quanta as pq

    try:
        # Create direction measure for the zenith
        zenith = pm.direction('AZEL','0deg','90deg')
    except AttributeError as e:
        if pm is None:
            raise ImportError("python-casacore import failed")

        raise

    # Create position measures for each antenna
    reference_positions = [pm.position('itrf',
        *(pq.quantity(x,'m') for x in pos))
        for pos in antenna_positions]

    # Compute field centre in radians
    fc_rad = pm.direction('J2000',
        *(pq.quantity(f,'rad') for f in field_centre))

    return np.asarray([
            # Set current time as the reference frame
            pm.do_frame(pm.epoch("UTC", pq.quantity(t, "s")))
            and
            [   # Set antenna position as the reference frame
                pm.do_frame(rp)
                and
                pm.posangle(fc_rad, zenith).get_value("rad")
                for rp in reference_positions
            ]
        for t in times])

def printJones(stnName,obsTimes,srcDir,freqs,model):
    print "Frequency","Time","J00","J01","J10","J11"
    for freq in freqs:
        Jn=getJonesByAntFld(model,obsTimes,stnName,srcDir,freq)
        #plotJonesGnuplot(quantity(obsTimes.get_value()[0],obsTimes.get_unit()).formatted("YMD"),stepTime.seconds,Jn)
        for ti in range(0,len(obsTimes.get_value() )):
            print freq, quantity(obsTimes.get_value()[ti],obsTimes.get_unit()).formatted("YMD"), Jn[ti,0,0], Jn[ti,0,1],Jn[ti,1,0],Jn[ti,1,1]

def plotJones(stnName,obsTimes,freqs):
    #frequencys=np.linspace(0,100e6,512)
    #freqs=frequencys[150:350]
    timespy=[datetime.fromtimestamp(quantity(t,'d').to_unix_time()) for t in obsTimes.get_value()]
    freq=freqs[0]
    Jn=getJonesByAntFld(model,obsTimes,stnName,srcDir,freq)
    p_ch = np.abs(Jn[:,0,0].squeeze())**2+np.abs(Jn[:,0,1].squeeze())**2
    q_ch = np.abs(Jn[:,1,1].squeeze())**2+np.abs(Jn[:,1,0].squeeze())**2
    #For testing purposes:
    #p_ch = np.real(Jn[:,0,1].squeeze())
    #q_ch = np.real(Jn[:,0,1].squeeze())
    #In dB:
    p_ch = 10*np.log10(p_ch)
    q_ch = 10*np.log10(q_ch)
    
    plt.figure()
    plt.subplot(211)
    plt.plot(timespy, p_ch)
    plt.title('p channel')
    #plt.clim(-9, -3)
    plt.subplot(212)
    plt.plot(timespy, q_ch)
    plt.title('q-channel')
    #plt.clim(-9, -3)
    plt.xlabel('Time')
    plt.show()

	def on_mouse_move(self,widget,event):
		if self._low_res_image != None: 
			sbrMain = self._builder.get_object("sbrMain")
			sbrMain.remove_all(sbrMain.get_context_id("CursorPos"))
			handle = self._builder.get_object("cvsLowRes")
			handle.queue_draw()
			rect = handle.get_allocation()
                        img_height = self._low_res_image.get_height()
                        img_width = self._low_res_image.get_width()
                        
                        self._l_pos = int(event.x/float(rect.width) * img_width) 
			self._m_pos = int(event.y/float(rect.height) * img_height) 
			
			sbrMain.push(sbrMain.get_context_id("CursorPos"),"Pixel position: (x,y) = (%d,%d), (ra,dec) = (%s,%s)" % (self._l_pos, self._m_pos,
							    		 quantity(quantity(self._phase_centres[self._field_id,0,0],"rad").get_value("arcsec") + (-self._l_pos + img_width/2)*self._img_cell_l,"arcsec").get("deg").formatted("[+-]dd.mm.ss.t.."),
							    		 quantity(quantity(self._phase_centres[self._field_id,0,1],"rad").get_value("arcsec") + (-self._m_pos + img_height/2)*self._img_cell_m,"arcsec").get("deg").formatted("[+-]dd.mm.ss.t..")))

def updateMSmetadata(msfile):
  #Update history to show that this script has modified original data
  tms = pt.table(msfile,readonly=False,ack=False)
  th = pt.table(tms.getkeyword('HISTORY'), readonly=False, ack=False)
  nr=th.nrows()
  th.addrows(1)
  tr=th.row()
  tr.put(nr,{'TIME': quantity('today').get('s').get_value(),
             'OBSERVATION_ID':0,
             'MESSAGE': 'Applied polarization modifications',
             'PRIORITY': 'NORMAL',
             'ORIGIN': '%s: version = %s' % (__file__,__version__),
             'OBJECT_ID':0, 
             'APPLICATION':__file__,
             'CLI_COMMAND':sys.argv,
             'APP_PARAMS': ['']})

  if not options.linear:
     #Change metadata information to be circular feeds
     feed = pt.table(tms.getkeyword('FEED'),readonly=False,ack=False)
     for tpart in feed.iter('ANTENNA_ID'):
      tpart.putcell('POLARIZATION_TYPE',0,['R','L'])

     polariz = pt.table(tms.getkeyword('POLARIZATION'),readonly=False,ack=False)
     polariz.putcell('CORR_TYPE',0,[5,6,7,8])
     tms.close()

 def epoch(self, rf='', v0='0.0d', off=False):
     loc = { 'type': 'epoch' , 'refer':  rf}
     loc['m0'] = dq.quantity(v0)
     if is_measure(off):
         if not off['type'] == "epoch":
             raise TypeError('Illegal offset type specified.')        
         loc["offset"] = off
     return self.measure(loc, rf)

def insert_grid_points(survey, filename, calibrator=False):
    with open(filename, 'r') as f:
        lines = f.readlines()

    for line in lines:
        split_string = line.split(None, 3)
        name = split_string[0]
        ra = quantity(split_string[1]).get_value('rad')
        dec = quantity(split_string[2].replace(":", ".", 2)).get_value('rad')
        if len(split_string) == 4:
            description = split_string[3].strip()
        else:
            description = ""
        f = Field.objects.create(
            name=name, ra=ra, dec=dec, description=description,
            survey=survey, calibrator=calibrator
        )

 def tofrequency(self, rf, v0, rfq):
     if is_measure(rfq) and rfq['type'] == 'frequency':
         rfq = dq.from_dict(rfq['m0'])
     if is_measure(v0) and  v0['type'] == 'doppler' \
            and  isinstance(rfq, dq._quanta.Quantity) \
            and  rfq.conforms(dq.quantity('Hz')):
         return self.doptofreq(v0, rf, dq.to_dict(rfq))
     else:
         raise TypeError('Illegal Doppler or rest frequency specified')

 def doppler(self, rf='', v0='0', off=False):
     loc = { 'type': "doppler",
             'refer': rf,
             'm0': dq.quantity(v0) }
     if is_measure(off):
         if not off['type'] == "doppler":
             raise TypeError('Illegal offset type specified.')        
         loc["offset"] = off
     return self.measure(loc, rf)

def GiveStrDate(tt):
    time_start = qa.quantity(tt, 's')
    me = pm.measures()
    dict_time_start_MDJ = me.epoch('utc', time_start)
    time_start_MDJ=dict_time_start_MDJ['m0']['value']
    JD=time_start_MDJ+2400000.5-2415020
    d=ephem.Date(JD)

    return d.datetime().isoformat().replace("T","/")

def updatehistory(outms):
  """
  Update history to show that this script has modified original data
  """
  tc = pt.table(outms,readonly=False)
  th = pt.table(tc.getkeyword('HISTORY'), readonly=False, ack=False)
  nr=th.nrows()
  th.addrows(1)
  tr=th.row()
  tr.put(nr,{'TIME': quantity('today').get('s').get_value(), 'OBSERVATION_ID':0,'MESSAGE': ' ', 'PRIORITY': ' ', 'ORIGIN': ' ','OBJECT_ID':0, 'APPLICATION':'mslin2circ','CLI_COMMAND':[''],'APP_PARAMS': ['']})

	def on_click(self,widget,event):
		if self._busy_loading:
		  message = Gtk.MessageDialog(self._builder.get_object("frmMain"),0,Gtk.MessageType.INFO,Gtk.ButtonsType.OK,"Busy synthesizing low res image... please hang on")
		  message.run()
		  message.destroy()
		  return
		if self._busy_faceting:
		  message = Gtk.MessageDialog(self._builder.get_object("frmMain"),0,Gtk.MessageType.INFO,Gtk.ButtonsType.OK,"Already busy faceting... please hang on")
		  message.run()
		  message.destroy()
		  return
		if self._low_res_image != None:
			self._busy_faceting = True
			self._builder.get_object("cvsLowRes").queue_draw()
			self.__change_visibilities()
			
			model = self._builder.get_object("lstFacetingProperties")
			itr = model.get_iter(0)
			facet_size_l = int(model.get_value(itr,1))
			itr = model.iter_next(itr)
			facet_size_m = int(model.get_value(itr,1))
			itr = model.iter_next(itr)
			facet_cell_l = float(model.get_value(itr,1))
			itr = model.iter_next(itr)
			facet_cell_m = float(model.get_value(itr,1))
			itr = model.iter_next(itr)
			conv_support = int(model.get_value(itr,1))
			itr = model.iter_next(itr)
			conv_oversample = int(model.get_value(itr,1))
			
			rect = self._builder.get_object("cvsLowRes").get_allocation()
                        img_height = self._low_res_image.get_height()
                        img_width = self._low_res_image.get_width()
                        facet_ra = quantity(self._phase_centres[self._field_id,0,0],"rad").get_value("arcsec") + (-int(event.x/float(rect.width) * img_width) + img_width/2)*self._img_cell_l
                        facet_dec = quantity(self._phase_centres[self._field_id,0,1],"rad").get_value("arcsec") + (-int(event.y/float(rect.height) * img_height) + img_height/2)*self._img_cell_m	
			facet_centres = np.array([[facet_ra,facet_dec]],dtype=np.float32)
			
			cmd_string = ("python bullseye.py \"%s\" --output_prefix \"%s\" --output_format png --npix_l %d --npix_m %d --cell_l %d"
				      " --cell_m %d --pol %s --conv_sup %d --conv_oversamp %d --facet_centres \(%f,%f\) --field_id %d --average_all 1" % (
				      self._ms_name,self.FACET_TMP_FILE_NAME,facet_size_l,facet_size_m,facet_cell_l,
				      facet_cell_m,self._polarization,conv_support,conv_oversample,
				      int(facet_ra),int(facet_dec),self._field_id))
			threading.Thread(target=self.cmd_call,args=(cmd_string,self.on_finished_faceting)).start()

    def rise(self, crd, ev='5deg'):
        if  not is_measure(crd):
            raise TypeError('No rise/set coordinates specified')
        ps = self._getwhere()
        self._fillnow()
        hd = self.measure(crd, "hadec")
        c = self.measure(crd, "app")
        evq = dq.quantity(ev)
        hdm1 = dq.from_dict(hd["m1"])
        psm1 = dq.from_dict(ps["m1"])
        ct = (dq.quantity(ev) - dq.sin(hdm1) * dq.sin(psm1)) / (dq.cos(hdm1) * dq.cos(psm1))

        if ct.get_value() >= 1:
            return "below below"
        if ct.get_value() <= -1:
            return "above above"
        a = dq.acos(ct)
        return { "rise": dq.norm(dq.quantity(c["m0"]), 0) - a,
                 "set" : dq.norm(dq.quantity(c["m0"]), 0) + a
                 }