text: import numpy
  import datetime
  import pandas
  
  class RT_SCT(object):
      """
      Class to convert ray-tracing results 
      into power in beam/range domain.
      An explanation for these calculations
      is given in Sebastien's dissertation.
      """
      
      def __init__(self, rt_dict, max_gates = 75):
          """
          Initialize the ray tracing params
          for a radar!
          rt_dict is obtained by running the RayTrace Class!
          """
          import rt_utils
          self.rt_dict = rt_dict
          self.max_gates = max_gates
          self.ranges = 180 + 45*numpy.arange(max_gates+1,dtype=numpy.int)
          self.cnst_dict = rt_utils.get_rt_consts()
          
      def get_gnd_sct_df(self):
          """
          Get estimated power from ground-scatter
          """
          lag_power_arr = []
          range_arr = []
          gates_arr = []
          date_arr = []
          azim_arr = []
          elev_arr = []
          df_arr = []
          print(self.rt_dict.keys())
          for _dt in self.rt_dict.keys():
              for _az in self.rt_dict[_dt].keys():
                  grp_ran = []
                  weights = []
                  sct_el_arr = []
                  
                  for _el in self.rt_dict[_dt][_az].keys():
                      gscat = self.rt_dict[_dt][_az][_el]['ground_scatter']
                      if gscat is None:
                          continue
                      sct_el_arr.append(_el)
                      _wght = 1/(gscat["grpran"]**3)
                      grp_ran.append(gscat["grpran"]*1e-3)
                      weights.append(_wght)
                  if len(grp_ran) == 0:
                      continue
                  lag_power, bins = numpy.histogram(\
                              grp_ran,\
                              bins=self.ranges,\
                              weights=weights) 
  #                 if _dt == datetime.datetime(2012,12,5,21):
  #                     print (lag_power,grp_ran)
  #                 lag_power = lag_power#10.*numpy.log10(lag_power/numpy.max(lag_power))
  #                 lag_power_arr += list(lag_power)
  #                 range_arr += list(bins)
                  _gates = list(numpy.arange(self.max_gates+1,dtype=numpy.int))
                  _dates = [_dt for x in range(len(list(bins)))]
                  _azims = [_az for x in range(len(list(bins)))]
                  _elevs = [numpy.median(sct_el_arr) for x in range(len(bins))]
                  _df = pandas.DataFrame(list(zip(lag_power, bins, _gates,\
                                 _dates, _azims, _elevs)),
                              columns =['lag_power', 'range', \
                               'gate', 'date', 'azim', 'median_elev'])
                  print(_df.head())
                  df_arr.append(_df)
                  
  #         gnd_sct_df = pandas.DataFrame(list(zip(\
  #                 date_arr, lag_power_arr, range_arr,\
  #                 gates_arr,azim_arr,elev_arr\
  #                     )), 
  #                    columns =['date', 'lag_power', 'range', \
  #                              'gate', 'azim', 'median_elev'])
          gnd_sct_df = pandas.concat(df_arr)
  #         print(gnd_sct_df[gnd_sct_df["date"] == datetime.datetime(2012,12,5,21)])
          # Now we need to normalize the power!
          # for comparing it with actual data (in some sense)
          # Refer to Sebastien's dissertation for 
          # more details!
          # IDL code: http://davit1.ece.vt.edu/doc/rt/rt_run-code.html
          if gnd_sct_df.shape[0] >= 1:
              gnd_sct_df['day'] = gnd_sct_df['date'].dt.date
              max_power_df = gnd_sct_df[['day', 'lag_power']].groupby(\
                                              ['day']\
                                              ).max().reset_index()
              max_power_df.columns = ['day', 'max_lag_power']
              gnd_sct_df = pandas.merge(gnd_sct_df, max_power_df, on=['day'])
              gnd_sct_df['lag_power'] = 10.*numpy.log10(\
                                              gnd_sct_df['lag_power']/gnd_sct_df['max_lag_power']\
                                              )
          
          return gnd_sct_df
          
      def get_iono_sct_df(self):
          """
          Get estimated power from ionospheric-scatter
          """
          weights_arr = []
          grp_ran_arr = []
          range_arr = []
          gates_arr = []
          aspect_arr = []
          alt_arr = []
          date_arr = []
          azim_arr = []
          elev_arr = []
          df_arr = []
  
          for _dt in self.rt_dict.keys():
              for _az in self.rt_dict[_dt].keys():
                  grp_ran = []
                  weights = []
                  aspect = []
                  altitude = []
                  sct_el_arr = []
                  for _el in self.rt_dict[_dt][_az].keys():
                      ionos = self.rt_dict[_dt][_az][_el]['ionospheric_scatter']
                      if ionos['nstep'] > 0.:
                          weights += list(ionos['weights'])
                          grp_ran += list(ionos['gran']*1e-3)
                          aspect += list(ionos['aspect'])
                          altitude += list((ionos['altitude'] - self.cnst_dict["Rav"])/1000.)
                          sct_el_arr.append(_el)
                  if len(weights) == 0.:
                      continue
  #                 adj_weights = weights/numpy.max(weights)
  #                 lag_power, bins = numpy.histogram(\
  #                                 grp_ran,\
  #                                 bins=self.ranges,\
  #                                 weights=adj_weights) 
  #                 lag_power = 10.*numpy.log10(lag_power/numpy.max(lag_power))
  #                 weights_arr += weights
  #                 grp_ran_arr += grp_ran
  # #                 range_arr += list(self.ranges)
  # #                 gates_arr += list(numpy.arange(self.max_gates+1,dtype=numpy.int))
  #                 date_arr += [_dt for x in range(len(grp_ran))]
  #                 azim_arr += [_az for x in range(len(grp_ran))]
  #                 elev_arr += [numpy.median(sct_el_arr) for x in range(len(grp_ran))]
  #                 aspect_arr += [numpy.median(aspect) for x in range(len(grp_ran))]
                  
                  _df = pandas.DataFrame(list(zip(weights, grp_ran,\
                                 [_dt for x in range(len(grp_ran))],\
                                 [_az for x in range(len(grp_ran))],\
                                 [numpy.median(sct_el_arr) for x in range(len(grp_ran))],\
                                 [numpy.median(altitude) for x in range(len(grp_ran))],\
                                 [numpy.median(aspect) for x in range(len(grp_ran))])),
                              columns = ['weights', 'grp_ran','date', \
                               'azim', 'median_elev', 'median_alt','median_aspect'])
                  df_arr.append(_df)
  
  #         iono_sct_df = pandas.DataFrame(list(zip(\
  #                 date_arr, weights_arr, grp_ran_arr,\
  #                 azim_arr, elev_arr, aspect_arr\
  #                     )), 
  #                    columns =['date', 'weights', 'grp_ran',\
  #                              'azim', 'median_elev', 'median_aspect'])
          iono_sct_df = pandas.concat(df_arr)
          
          if iono_sct_df.shape[0] >= 1:
              iono_sct_df['day'] = iono_sct_df['date'].dt.date        
              # Now we need to normalize the weights by day
              # and calculate the lag_power, and then recreate 
              # the dataframe as we take a histogram 
              # (of group ranges and weights) to calc
              # lag_power and it has a different length
              # compared to the group ranges.
              max_wght_df = iono_sct_df[['day', 'weights']].groupby(\
                                      ['day']\
                                      ).max().reset_index()
              max_wght_df.columns = ['day', 'max_weights']
              iono_sct_df = pandas.merge(iono_sct_df, max_wght_df, on=['day'])
              iono_sct_df["norm_weights"] = iono_sct_df["weights"]/iono_sct_df["max_weights"]
  
              lag_power_arr = []
              gates_arr = []
              range_arr = []
              date_arr = []
              azim_arr = []
              elev_arr = []
              aspect_arr = []
              alt_arr = []
              
              fin_df_arr = []
  
              max_gates = 75
              grps = iono_sct_df.groupby(["date", "azim"])
              ranges = 180 + 45*numpy.arange(75+1,dtype=numpy.int)
              addtnl_range = numpy.max(ranges) + 45
              for _d, _az in list(grps.indices.keys()):
                  sel_df = iono_sct_df[ (iono_sct_df["date"] == _d) & (iono_sct_df["azim"] == _az) ]
  
                  lag_power, bins = numpy.histogram( sel_df["grp_ran"], bins=self.ranges,\
                                                    weights=sel_df["norm_weights"])
                  _curr_el = sel_df["median_elev"].unique()[0]
                  _curr_asp = sel_df["median_aspect"].unique()[0]
                  _curr_alt = sel_df["median_alt"].unique()[0]
                  # Note the numbers of bins is one less
                  # than the lag_power, since we are getting
                  # a histogram. So When onverting to a DF we'll
                  # need an additional element in the rage!
  #                 lag_power_arr += list(lag_power)
  #                 gates_arr += list(numpy.arange(1,self.max_gates+1,dtype=numpy.int))# 
  #                 range_arr += list(self.ranges[:-1])# 
  #                 azim_arr += [_az for x in range(lag_power.shape[0])]
  #                 date_arr += [_d for x in range(lag_power.shape[0])]
  #                 elev_arr += [_curr_el for x in range(lag_power.shape[0])]
  #                 aspect_arr += [_curr_asp for x in range(lag_power.shape[0])]
                  
                  
                  _gates = list(numpy.arange(1,self.max_gates+1,dtype=numpy.int))# 
                  _dates = [_d for x in range(lag_power.shape[0])]
                  _ranges = list(self.ranges[:-1])
                  _azims = [_az for x in range(len(list(bins)))]
                  _elevs = [_curr_el for x in range(lag_power.shape[0])]
                  _aspects = [_curr_asp for x in range(lag_power.shape[0])]
                  _alts = [_curr_alt for x in range(lag_power.shape[0])]
                  
                  _df = pandas.DataFrame(list(zip(_dates, lag_power, _gates, _ranges,\
                                 _azims, _elevs, _aspects, _alts)),
                              columns =['date', 'lag_power', 'gates', 'range',\
                                       'azim', 'median_elev', 'median_aspect', 'median_alt'])
                  fin_df_arr.append(_df)
                  
                  
              iono_df = pandas.concat(fin_df_arr)
  #             iono_df = pandas.DataFrame(list(zip(\
  #                         date_arr, lag_power_arr, gates_arr, range_arr,\
  #                         azim_arr, elev_arr, aspect_arr\
  #                             )), 
  #                            columns =['date', 'lag_power', 'gates', 'range',\
  #                                      'azim', 'median_elev', 'median_aspect'])
  
              iono_df["day"] = iono_df["date"].dt.date
              max_power_df = iono_df[['day', 'lag_power']].groupby(\
                                              ['day']\
                                              ).max().reset_index()
              max_power_df.columns = ['day', 'max_lag_power']
              iono_df = pandas.merge(iono_df, max_power_df, on=['day'])
  
              iono_df['lag_power'] = 10.*numpy.log10(\
                                      iono_df['lag_power']/iono_df['max_lag_power']\
                                      )
  
              return iono_df
          else:
              return iono_sct_df
          
      def get_ionos_scatter(self,tol=1.):
          """
          This is a redundant function and may not be needed.
          Basically, here we identify the rays which with 
          good aspect conditions!
          """
          tht0 = 2000./6370./2.
          Rav = 6370.
          iono_rays = {}
          for _dt in self.rt_dict.keys():
              iono_rays[_dt] = {}
              for _az in self.rt_dict[_dt].keys():
                  iono_rays[_dt][_az] = {}
                  for _el in self.rt_dict[_dt][_az].keys():
                      iono_rays[_dt][_az][_el] = {}
                      dip = self.rt_dict[_dt][_az][_el]['dip']
                      ndip = dip.shape[0]
                      thtdip = numpy.arange(ndip)*2500./Rav/ndip
                      rays = self.rt_dict[_dt][_az][_el]['ray_parameters']
                      # single param keys
                      iono_rays[_dt][_az][_el]['nsteps'] = rays['nsteps']
                      iono_rays[_dt][_az][_el]['rayhour'] = rays['rayhour']
                      iono_rays[_dt][_az][_el]['rayazim'] = rays['rayazim']
                      iono_rays[_dt][_az][_el]['rayelev'] = rays['rayelev']
                      # list keys
                      iono_rays[_dt][_az][_el]['r'] = []
                      iono_rays[_dt][_az][_el]['th'] = []
                      iono_rays[_dt][_az][_el]['gran'] = []
                      iono_rays[_dt][_az][_el]['nr'] = []
                      for _n in range(1,len(rays['r'])):
                          if (rays['gran'][_n] > 0.) & (rays['th'][_n] < 2000./6370.) & (rays['r'][_n]*1e-3 > 6370.+80.):
                              # Calculate k vector
                              kx = rays['r'][_n]*numpy.sin(rays['th'][_n]-rays['th'][_n-1])
                              kz = rays['r'][_n]*numpy.cos(rays['th'][_n]-rays['th'][_n-1]) - rays['r'][_n-1]
                              kvect = numpy.sqrt( kx**2. + kz**2. )
                              # Middle of the step: position and index in B grid
                              midtht = (rays['th'][_n]-rays['th'][_n-1])/2. + rays['th'][_n-1]
                              thtind = numpy.argmin(numpy.abs(midtht-thtdip))
                      #         print(midtht,thtind,thtdip[thtind])
                              middip = dip[thtind,0]
                              middec = dip[thtind,1]
                              #calculate vector magnetic field
                              Bx = numpy.cos(numpy.deg2rad(-1*middip)) * numpy.cos(\
                                              numpy.deg2rad(rays['rayazim']) - numpy.deg2rad(middec)\
                                                          )
                              Bz = numpy.sin(numpy.deg2rad(-1*middip))
                              # calculate cosine of aspect angle
                              cos_aspect = (Bx*kx + Bz*kz)/kvect
                              if (abs(cos_aspect) < numpy.cos(numpy.pi/2. - numpy.deg2rad(tol))):
                                  'nsteps', 'rayhour', 'rayazim', 'rayelev', 'r', 'th', 'gran', 'nr'
                                  iono_rays[_dt][_az][_el]['r'].append(rays['r'][_n])
                                  iono_rays[_dt][_az][_el]['th'].append(rays['th'][_n])
                                  iono_rays[_dt][_az][_el]['gran'].append(rays['gran'][_n])
                                  iono_rays[_dt][_az][_el]['nr'].append(rays['nr'][_n])
          return iono_rays
  