From afe459e0ead458b2cbd45a339f4824aacf0c2bc4 Mon Sep 17 00:00:00 2001
From: Carlos Paniagua <cpaniaguam@gmail.com>
Date: Sun, 14 Jan 2024 22:02:33 -0500
Subject: [PATCH] chore: update class definitions

---
 analyis_publish/PB05_define_angle.py          | 198 ++----------------
 analyis_publish/PB05_define_angle_supl.py     | 170 ++-------------
 analysis/SB04_2d_wavefield_emulator.py        |  61 +-----
 analysis_db/B05_define_angle.py               |   2 +-
 .../ICEsat2_SI_tools/angle_optimizer.py       |   2 +-
 .../ICEsat2_SI_tools/generalized_FT.py        |   4 +-
 src/icesat2_tracks/ICEsat2_SI_tools/io.py     |   5 +-
 .../ICEsat2_SI_tools/spectral_estimates.py    |  10 +-
 .../X03_MCMC_surface_smapling.py              |   2 +-
 .../local_modules/m_colormanager_ph3.py       |   2 +-
 .../local_modules/m_general_ph3.py            |  16 +-
 .../local_modules/m_spectrum_ph3.py           |   8 +-
 12 files changed, 66 insertions(+), 414 deletions(-)

diff --git a/analyis_publish/PB05_define_angle.py b/analyis_publish/PB05_define_angle.py
index c3226674..fc75cebb 100644
--- a/analyis_publish/PB05_define_angle.py
+++ b/analyis_publish/PB05_define_angle.py
@@ -1,17 +1,13 @@
-# %%
-import os, sys
-#execfile(os.environ['PYTHONSTARTUP'])
-
 """
 This file open a ICEsat2 track applied filters and corections and returns smoothed photon heights on a regular grid in an .nc file.
 This is python 3
 """
 
+import os, sys
+
 exec(open(os.environ['PYTHONSTARTUP']).read())
 exec(open(STARTUP_2021_IceSAT2).read())
 
-#%matplotlib inline
-
 import ICEsat2_SI_tools.convert_GPS_time as cGPS
 import h5py
 import ICEsat2_SI_tools.io as io
@@ -36,50 +32,27 @@
 col.colormaps2(21)
 
 col_dict = col.rels
-#import s3fs
-# %%
 track_name, batch_key, test_flag = io.init_from_input(sys.argv) # loads standard experiment
-#track_name, batch_key, test_flag = '20190605061807_10380310_004_01', 'SH_batch01', False
-#track_name, batch_key, test_flag = '20190601094826_09790312_004_01', 'SH_batch01', False
-#track_name, batch_key, test_flag = '20190207111114_06260210_004_01', 'SH_batch02', False
-#track_name, batch_key, test_flag = '20190219073735_08070210_004_01', 'SH_batch02', False
-#track_name, batch_key, test_flag = '20190215184558_07530210_004_01', 'SH_batch02', False
 
-# good track
-#track_name, batch_key, test_flag = '20190502021224_05160312_004_01', 'SH_batch02', False
-#track_name, batch_key, test_flag = '20190502050734_05180310_004_01', 'SH_batch02', False
-#track_name, batch_key, test_flag = '20190216200800_07690212_004_01', 'SH_batch02', False
-
-#track_name, batch_key, test_flag = '20190213133330_07190212_004_01', 'SH_batch02', False
 
 # main text figure
 track_name, batch_key, test_flag = 'SH_20190502_05160312', 'SH_publish', False
 
 #suppl. figures:
-#track_name, batch_key, test_flag = 'SH_20190219_08070210', 'SH_publish', False
-
-#print(track_name, batch_key, test_flag)
 hemis, batch = batch_key.split('_')
-#track_name= '20190605061807_10380310_004_01'
 ATlevel= 'ATL03'
 
 
 
-#plot_path   = mconfig['paths']['plot'] + '/'+hemis+'/'+batch_key+'/' + track_name + '/B05_angle/'
 plot_path   = mconfig['paths']['plot'] + '/'+hemis+'/'+batch_key+'/publish/' + track_name + '/'
 MT.mkdirs_r(plot_path)
 bad_track_path =mconfig['paths']['work'] +'bad_tracks/'+ batch_key+'/'
-# %%
 
 all_beams   = mconfig['beams']['all_beams']
 high_beams  = mconfig['beams']['high_beams']
 low_beams   = mconfig['beams']['low_beams']
 beam_groups = mconfig['beams']['groups']
 group_names = mconfig['beams']['group_names']
-#Gfilt   = io.load_pandas_table_dict(track_name + '_B01_regridded', load_path) # rhis is the rar photon data
-
-# load_path   = mconfig['paths']['work'] +'/B01_regrid_'+hemis+'/'
-# G_binned    = io.load_pandas_table_dict(track_name + '_B01_binned' , load_path)  #
 load_path   = mconfig['paths']['work'] +batch_key +'/B02_spectra/'
 Gk          = xr.load_dataset(load_path+ '/B02_'+track_name + '_gFT_k.nc' )  #
 
@@ -91,22 +64,6 @@
 Prior       = MT.load_pandas_table_dict('/A02_'+track_name, load_path)['priors_hindcast']
 
 
-# font_for_print()
-# F = M.figure_axis_xy(5.5, 3, view_scale= 0.8)
-# plt.suptitle(track_name)
-# ax1 =  plt.subplot(2, 1, 1)
-# plt.title('Data in Beam', loc= 'left')
-#
-# xi =1
-
-#data = Marginals.isel(x=xi).sel(beam_group= 'group1').marginals
-# angle_mask = Marginals.angle[2:-2]
-#
-#data.T.plot(cmap= plt.cm.OrRd)
-
-# %%
-
-
 def derive_weights(weights):
     weights = (weights-weights.mean())/weights.std()
     weights = weights - weights.min()
@@ -127,7 +84,6 @@ def weighted_means(data, weights, x_angle, color='k'):
         k           = wi.k.data
         data_k      = data.sel(k=k).squeeze()
         data_weight = (data_k * wi)
-        #plt.stairs(data_weight.sum('k')/ weight_norm , x_angle, linewidth=1 , color ='k')
         if data_k.k.size > 1:
             for k in data_k.k.data:
                 plt.stairs(data_weight.sel(k=k) / weight_norm, x_angle, color ='gray', alpha =0.5)
@@ -190,10 +146,6 @@ def weighted_means(data, weights, x_angle, color='k'):
     group_weight = Gweights.isel(x =xi)
 
     ax_list= dict()
-    #ax_sum = F.fig.add_subplot(gs[1, 0])
-    # #ax_sum.tick_params(labelbottom=False)
-    #
-    # ax_list['sum'] = ax_sum
 
     data_collect = dict()
     for group, gpos in zip(Marginals.beam_group.data, [ gs[0, 0], gs[1, 0], gs[2, 0]] ):
@@ -215,9 +167,7 @@ def weighted_means(data, weights, x_angle, color='k'):
         else:
             data_wmean = weighted_means(data, weights, x_angle, color= col_dict[group] )
             plt.stairs(data_wmean , x_angle, color =col_dict[group], alpha =1)
-        # test if density is correct
-        # if np.round(np.trapz(data_wmean) * d_angle, 2) < 0.90:
-        #     raise ValueError('weighted mean is not a density anymore')
+
 
         if group == 'group1':
             t_string = group_names[group] +' pair' #group.replace('group',
@@ -225,31 +175,18 @@ def weighted_means(data, weights, x_angle, color='k'):
             t_string = group_names[group]+' pair'  #group.replace('group', +' ')
 
         plt.title(next(fn) + t_string, loc ='left')
-        #plt.sca(ax_sum)
-
-        # if data_collect is None:
-        #     data_collect = data_wmean
-        # else:
         data_collect[group] = data_wmean
-        #ax0.set_yscale('log')
+
 
     data_collect = xr.concat(data_collect.values(), dim='beam_group')
     final_data   = (group_weight * data_collect).sum('beam_group')/group_weight.sum('beam_group').data
 
-    # plt.sca(ax_sum)
-    # plt.stairs( final_data , x_angle, color = 'k', alpha =1, linewidth =0.8)
-    # ax_sum.set_xlabel('Angle (rad)')
-    # plt.title('Weighted mean over group & wavenumber', loc='left')
-
     # get relevant priors
     for axx in ax_list.values():
         axx.set_ylim(0, final_data.max() * 1.5)
-        #figureaxx.set_yscale('log')
         axx.set_xticks(xticks_pi)
         axx.set_xticklabels(xtick_labels_pi)
         axx.set_xlim(-np.pi/2, np.pi/2)
-        #ax_final.set_xticks(xticks_pi)
-        #ax_final.set_xticklabels(xtick_labels_pi)
 
 
     try:
@@ -257,7 +194,6 @@ def weighted_means(data, weights, x_angle, color='k'):
         ax_list['group2'].set_ylabel('PDF')
         ax_list['group3'].set_ylabel('PDF')
         ax_list['group1'].tick_params(labelbottom=True)
-        #ax_list['group3'].set_xlabel('Angle (rad)')
     except:
         pass
 
@@ -271,9 +207,6 @@ def weighted_means(data, weights, x_angle, color='k'):
     plt.stairs( final_data , x_angle, color = 'k', alpha =0.5, linewidth =0.8, zorder= 12)
 
     final_data_smth = lanczos.lanczos_filter_1d(x_angle,final_data, 0.1)
-    #
-    # for group in Marginals.beam_group.data:
-    #     plt.stairs( data_collect.sel(beam_group= group) * group_weight.sel(beam_group= group) /group_weight.sum('beam_group').data, x_angle, color =col_dict[group], alpha =1)
 
     plt.plot(x_angle[0:-1], final_data_smth, color = 'black', linewidth= 0.8)
 
@@ -299,9 +232,6 @@ def weighted_means(data, weights, x_angle, color='k'):
 Gpdf = xr.merge([M_final,M_final_smth])
 
 Gpdf.weighted_angle_PDF_smth.plot()
-#Gpdf.isel( x=slice(0, 3 )).weighted_angle_PDF_smth.mean('x')
-#Gpdf.angle[Gpdf.mean('x').weighted_angle_PDF_smth.argmax()].data
-
 Gpdf.mean('x').weighted_angle_PDF_smth.plot()
 best_guess_angle = Gpdf.angle[Gpdf.mean('x').weighted_angle_PDF_smth.argmax()].data
 
@@ -309,19 +239,8 @@ def weighted_means(data, weights, x_angle, color='k'):
 
 best_guess_angle/np.pi
 Gpdf.mean('x').weighted_angle_PDF_smth.plot()
-#Gpdf.weighted_angle_PDF.where(~np.isnan(Gpdf.weighted_angle_PDF),0 ).plot()
-
-# if len(Gpdf.x) < 2:
-#     print('not enough x data, exit')
-#     MT.json_save('B05_fail', plot_path+'../',  {'time':time.asctime( time.localtime(time.time()) ) , 'reason': 'not enough x segments'})
-#     print('exit()')
-#     exit()
-
-
 
-
-# %%
-class plot_polarspectra(object):
+class plot_polarspectra:
         def __init__(self,k, thetas, data, data_type='fraction' ,lims=None,  verbose=False):
 
             """
@@ -334,22 +253,21 @@ def __init__(self,k, thetas, data, data_type='fraction' ,lims=None,  verbose=Fal
 
             #self.sample_unit=sample_unit if sample_unit is not None else 'df'
             # decided on freq limit
-            self.lims= lims = [self.k.min(),self.k.max()] if lims is None else lims #1.0 /lims[1], 1.0/ lims[0]
+            self.lims= lims = [self.k.min(),self.k.max()] if lims is None else lims
             freq_sel_bool=M.cut_nparray(self.k, lims[0], lims[1] )
 
-            self.min=np.round(np.nanmin(data[freq_sel_bool,:]), 2)#*0.5e-17
+            self.min=np.round(np.nanmin(data[freq_sel_bool,:]), 2)
             self.max=np.round(np.nanmax(data[freq_sel_bool,:]), 2)
             if verbose:
                 print(str(self.min), str(self.max) )
 
-            self.klabels=np.linspace(self.min, self.max, 5) #np.arange(10, 100, 20)
+            self.klabels=np.linspace(self.min, self.max, 5)
 
             self.data_type=data_type
             if data_type == 'fraction':
                 self.clevs=np.linspace(np.nanpercentile(dir_data.data, 1), np.ceil(self.max* 0.9), 21)
             elif data_type == 'energy':
                 self.ctrs_min=self.min+self.min*.05
-                #self.clevs=np.linspace(self.min, self.max, 21)
                 self.clevs=np.linspace(self.min+self.min*.05, self.max*.60, 21)
 
 
@@ -359,35 +277,27 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
             """
             if ax is None:
                 ax          = plt.subplot(111, polar=True)
-                #self.title  = plt.suptitle('  Polar Spectrum', y=0.95, x=0.5 , horizontalalignment='center')
             else:
                 ax=ax
-            ax.set_theta_direction(-1)  #right turned postive
+            ax.set_theta_direction(-1)
             ax.set_theta_zero_location("W")
 
             grid=ax.grid(color='k', alpha=.5, linestyle='-', linewidth=.5)
 
             if self.data_type == 'fraction':
-                cm=plt.cm.RdYlBu_r #brewer2mpl.get_map( 'RdYlBu','Diverging', 4, reverse=True).mpl_colormap
+                cm=plt.cm.RdYlBu_r
                 colorax = ax.contourf(self.thetas,self.k, self.data, self.clevs, cmap=cm, zorder=1)# ,cmap=cm)#, vmin=self.ctrs_min)
             elif self.data_type == 'energy':
-                cm=plt.cm.Paired#brewer2mpl.get_map( 'Paired','Qualitative', 8).mpl_colormap
+                cm=plt.cm.Paired
                 cm.set_under='w'
                 cm.set_bad='w'
-                colorax = ax.contourf(self.thetas,self.k, self.data, self.clevs, cmap=cm, zorder=1)#, vmin=self.ctrs_min)
-            #divider = make_axes_locatable(ax)
-            #cax = divider.append_axes("right", size="5%", pad=0.05)
+                colorax = ax.contourf(self.thetas,self.k, self.data, self.clevs, cmap=cm, zorder=1)
             self.colorax = colorax
 
             if cbar_flag:
                 cbar = plt.colorbar(colorax, fraction=0.046, pad=0.1, orientation="horizontal")
-                # if self.data_type == 'fraction':
-                #     cbar.set_label('Energy Distribution', rotation=0, fontsize=fontsize)
-                # elif self.data_type == 'energy':
-                #     cbar.set_label('Energy Density ('+self.unit+')', rotation=0, fontsize=fontsize)
                 cbar.ax.get_yaxis().labelpad = 30
                 cbar.outline.set_visible(False)
-                #cbar.ticks.
                 clev_tick_names, clev_ticks =MT.tick_formatter(FP.clevs, expt_flag= False, shift= 0, rounder=4, interval=1)
                 cbar.set_ticks(clev_ticks[::5])
                 cbar.set_ticklabels(clev_tick_names[::5])
@@ -399,7 +309,6 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
                 radial_ticks = np.arange(60, 1000, 20)
             print(radial_ticks)
             xx_tick_names, xx_ticks = MT.tick_formatter( radial_ticks , expt_flag= False, shift= 1, rounder=0, interval=1)
-            #xx_tick_names, xx_ticks = MT.tick_formatter( np.arange( np.floor(self.k.min()),self.k.max(), 20) , expt_flag= False, shift= 1, rounder=0, interval=1)
             xx_tick_names = ['  '+str(d)+'m' for d in xx_tick_names]
 
             ax.set_yticks(xx_ticks[::1])
@@ -415,26 +324,22 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
             lines, labels = plt.thetagrids(degrange, labels=degrange_label)#, frac = 1.07)
 
             for line in lines:
-                #L=line.get_xgridlines
                 line.set_linewidth(5)
-                #line.set_linestyle(':')
 
-            #ax.set_yscale('log')
             ax.set_ylim(self.lims)
             ax.spines['polar'].set_color("none")
             ax.set_rlabel_position(87)
             self.ax=ax
 
-# %%
 font_for_print()
 fn = copy.copy(lstrings)
 
 
 F = M.figure_axis_xy(fig_sizes['two_column_square'][0], fig_sizes['two_column_square'][1], view_scale= 0.7, container = True)
-gs = GridSpec(8,6,  wspace=0.1,  hspace=2.1)#figure=fig,
+gs = GridSpec(8,6,  wspace=0.1,  hspace=2.1)
 col.colormaps2(21)
 
-cmap_spec= col.white_base_blgror #plt.cm.ocean_r
+cmap_spec= col.white_base_blgror
 clev_spec = np.linspace(-8, -1, 21) *10
 
 cmap_angle= col.cascade_r
@@ -452,26 +357,18 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 k_low = (k_low_limits + k_low_limits.diff('k')[0]/2).data
 weighted_spec_sub['k_bins'] = k_low[0:-1]
 weighted_spec_sub = weighted_spec_sub.rename({'k_bins': 'k'})
-#weighted_spec_sub = weighted_spec
 
 lam_p = 2 *np.pi/k_low_limits
 lam = lam_p * np.cos(best_guess_angle)
 k               = 2 * np.pi/lam
-
-#weighted_spec.k/np.cos(best_guess_angle)
-
-#xlims = x_spec[0]-12.5/2, x_spec[-5]
 xlims = x_spec[0], x_spec[-5]
 
-#weighted_spec.plot()
-#clev_spec = np.linspace(-8, -1, 21) *10
 clev_spec = np.linspace(-80, (10* np.log(weighted_spec)).max() * 0.9, 21)
 
 
 
 dd = 10* np.log(weighted_spec_sub)
 clev_log = M.clevels( [dd.quantile(0.01).data * 0.3, dd.quantile(0.98).data * 2.5], 31)* 1
-#plt.pcolor(x_spec, k, dd ,vmin= clev_spec[0], vmax= clev_spec[-1],  cmap =cmap_spec )
 plt.pcolormesh(x_spec, lam, dd, cmap=cmap_spec , vmin = clev_log[0], vmax = clev_log[-1])
 
 
@@ -486,19 +383,14 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 cbar = plt.colorbar( fraction=0.018, pad=0.01, orientation="vertical", label ='Power')
 cbar.outline.set_visible(False)
 clev_ticks = np.round(clev_spec[::3], 0)
-#clev_tick_names, clev_ticks =MT.tick_formatter(clev_spec, expt_flag= False, shift= 0, rounder=1, interval=2)
 cbar.set_ticks(clev_ticks)
 cbar.set_ticklabels(clev_ticks)
 
 plt.ylabel('corrected wavelength $(m)$')
-#plt.xlabel('x (km)')
-
-#plt.colorbar()
 ax2 = F.fig.add_subplot(gs[3:5, :])
 ax2.tick_params(labelleft=True)
 
-#Gpdf.weighted_angle_PDF.where(~np.isnan(Gpdf.weighted_angle_PDF),0 ).T.plot()
-dir_data = Gpdf.interp(x= weighted_spec.x).weighted_angle_PDF_smth.T#.rolling(angle=5, min_periods= 1, center=True).mean()
+dir_data = Gpdf.interp(x= weighted_spec.x).weighted_angle_PDF_smth.T
 
 x = Gpdf.x/1e3
 angle = Gpdf.angle[::10]
@@ -508,24 +400,16 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 dir_data_sub['angle_bins']  = angle_low[0:-1]
 dir_data_sub            = dir_data_sub.rename({'angle_bins': 'angle'})
 plt.pcolormesh(dir_data_sub.x/1e3, dir_data_sub.angle, dir_data_sub , vmin= clev_angle[0], vmax= clev_angle[-1], cmap = cmap_spec)
-
-#plt.pcolormesh(dir_data.x/1e3, dir_data.angle, dir_data , vmin= clev_angle[0], vmax= clev_angle[-1], cmap = cmap_spec)
-
-
 cbar = plt.colorbar( fraction=0.02, pad=0.01, orientation="vertical", label ='Density')
 cbar.outline.set_visible(False)
 plt.title(next(fn) + 'Direction PDFs', loc='left')
-
-
 plt.ylabel('Angle')
 plt.xlabel('X (km)')
 
-
 ax2.set_yticks(xticks_pi)
 ax2.set_yticklabels(xtick_labels_pi)
 ax2.set_ylim(angle[0], angle[-1])
 
-
 x_ticks  = np.arange(0, xlims[-1].data, 25)
 x_tick_labels, x_ticks = MT.tick_formatter(x_ticks, expt_flag= False, shift= 0, rounder=1, interval=2)
 
@@ -534,32 +418,19 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 ax1.set_xticklabels(x_tick_labels)
 ax2.set_xticklabels(x_tick_labels)
 
-#ax1.set_yscale('log')
 lam_lim= lam[-1].data, 550
 ax1.set_ylim(lam_lim)
 
 ax1.set_xlim(xlims)
 ax2.set_xlim(xlims)
-#ax2.set_yscale('log')
 ax2.axhline(best_guess_angle, color=col.orange, linewidth=0.8)
 
 
-
-#xx_list = np.insert(weighted_spec.x.data, 0, 0)
-# x_pos_list = spec.create_chunk_boundaries( 1,  xx_list.size,  iter_flag= False )
-# #x_pos_list = spec.create_chunk_boundaries( int(xx_list.size/3),  xx_list.size,  iter_flag= False )
-# x_pos_list = x_pos_list[:, ::2]
-# x_pos_list[-1, -1] = xx_list.size-1re
-#x_pos_list#.shape
-
-x_pos_list =  [0, 1, 2]#np.arange(0,9, 1)#np.vstack([np.arange(1,3), np.arange(0,3)+1])
-#x_pos_list
+x_pos_list =  [0, 1, 2]
 lsrtrings = iter(['c)', 'd)', 'e)'])
 
 dir_ax_list =list()
 for x_pos, gs in zip( x_pos_list , [ gs[-3:, 0:2], gs[-3:, 2:4], gs[-3:, 4:]] ):
-    #print( x_pos)
-    #print( xx_list[x_pos])
     x_range = weighted_spec.x.data[x_pos] + 12.5e3/2 #, x_pos[-1]]]
     print(x_range)
     ax1.axvline(x_range/1e3, linestyle= '-', color= col.green, linewidth=0.9, alpha = 0.8)
@@ -567,37 +438,17 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 
     i_lstring = next(lsrtrings)
     ax1.text(x_range/1e3, np.array(lam_lim).mean()* 1.2 * 3/2, ' '+ i_lstring, fontsize= 8, color =col.green)
-    #ax2.text(x_range/1e3, weighted_spec.k.mean().data, ' a', fontsize= 8)
-
-
-    # ax1.axvline(x_range[-1]/1e3, color = 'gray', alpha = 0.5)
-    #
-    # ax2.axvline(x_range[0]/1e3, linestyle= ':', color= 'white', alpha = 0.5)
-    # ax2.axvline(x_range[-1]/1e3, color = 'gray', alpha = 0.5)
 
 
-    # i_spec  = weighted_spec.sel(x= slice(x_range[0], x_range[-1]) )
-    # i_dir   = corrected_marginals.sel(x= slice(x_range[0], x_range[-1]) )
     i_spec  = weighted_spec.isel(x= x_pos )
     i_dir   = corrected_marginals.interp(x= weighted_spec.x).isel(x= x_pos )
     print(i_spec.x.data, i_spec.x.data)
     dir_data  = (i_dir * i_dir.N_data).sum([ 'beam_group'])/ i_dir.N_data.sum([ 'beam_group'])
     lims = dir_data.k[ (dir_data.sum('angle')!=0) ][0].data, dir_data.k[ (dir_data.sum('angle')!=0)  ][-1].data
 
-    #dir_data.plot()
-    #dir_data.rolling(angle =5,  min_periods= 1, center=True ).mean().plot()
-
     N_angle = i_dir.angle.size
-    dir_data2 =  dir_data#.where( dir_data.sum('angle') !=0, 1/N_angle/d_angle )
-
-    plot_data  = dir_data2  * i_spec#.mean('x')
-    
-    # angle_low = dir_data2.angle[::5]
-    # k_low = dir_data2.k[::5]
-    # plot_data  = dir_data2.groupby_bins('angle' , angle_low).mean().groupby_bins('k', k_low).mean()
-    # plot_data = plot_data.rename({'k_bins':'k', 'angle_bins': 'angle'})
-    # plot_data['k'] = (k_low + k_low.diff('k')[0]/2).data[0:-1]
-    # plot_data['angle'] =(angle_low + angle_low.diff('angle')[0]/2).data[0:-1]
+    dir_data2 =  dir_data
+    plot_data  = dir_data2  * i_spec
     
     plot_data  = dir_data2.rolling(angle =2, k =15, min_periods= 1, center=True ).median()  * i_spec#.mean('x')
     plot_data = plot_data.sel(k=slice(lims[0],lims[-1] ) )
@@ -605,16 +456,12 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
     lam_p = 2 *np.pi/plot_data.k.data
     lam = lam_p * np.cos(best_guess_angle)
 
-    #F = M.figure_axis_xy(5, 4)
-    #ax = plt.subplot(1, 1, 1, polar=True)
-    #
     if np.nanmax(plot_data.data) != np.nanmin(plot_data.data):
 
         ax3 = F.fig.add_subplot(gs, polar=True)
         FP= plot_polarspectra(lam, plot_data.angle, plot_data, lims=[lam[-1], 138 ] , verbose= False, data_type= 'fraction')
         FP.clevs=np.linspace(np.nanpercentile(plot_data.data, 1), np.round(plot_data.max(), 4), 21)
         FP.linear(ax = ax3, cbar_flag=False)
-        #plt.show()
         plt.title('\n\n'+i_lstring,y=1.0, pad=-6, color=col.green)
 
         dir_ax_list.append(ax3)
@@ -630,11 +477,7 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 cbar.set_label('Energy Density \n(10$^3$ (m/m)$^2$ k$^{-1}$ deg$^{-1}$ )', rotation=90)#, fontsize=10)
 
 
-# F.save_pup(path = plot_path, name = 'B05_dir_ov_'+track_name)
-# F.save_light(path = plot_path, name = 'B05_dir_ov_'+track_name)
-
-
-# %% shift simple
+# shift simple
 
 font_for_print()
 fn = copy.copy(lstrings)
@@ -658,6 +501,3 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 plt.ylabel('$m^2/\lambda$')
 F.save_pup(path = plot_path, name = 'B05_dir_ov_'+track_name+'_1d')
 
-# %%
-#F.save_pup(path = plot_path, name = 'B05_dir_ov_'+track_name)
-# MT.json_save('B05_success', plot_path + '../', {'time':time.asctime( time.localtime(time.time()) )})
diff --git a/analyis_publish/PB05_define_angle_supl.py b/analyis_publish/PB05_define_angle_supl.py
index 315248dc..7f0880ae 100644
--- a/analyis_publish/PB05_define_angle_supl.py
+++ b/analyis_publish/PB05_define_angle_supl.py
@@ -36,51 +36,23 @@
 col.colormaps2(21)
 
 col_dict = col.rels
-#import s3fs
-# %%
 track_name, batch_key, test_flag = io.init_from_input(sys.argv) # loads standard experiment
-#track_name, batch_key, test_flag = '20190605061807_10380310_004_01', 'SH_batch01', False
-#track_name, batch_key, test_flag = '20190601094826_09790312_004_01', 'SH_batch01', False
-#track_name, batch_key, test_flag = '20190207111114_06260210_004_01', 'SH_batch02', False
-#track_name, batch_key, test_flag = '20190219073735_08070210_004_01', 'SH_batch02', False
-#track_name, batch_key, test_flag = '20190215184558_07530210_004_01', 'SH_batch02', False
-
-# good track
-#track_name, batch_key, test_flag = '20190502021224_05160312_004_01', 'SH_batch02', False
-#track_name, batch_key, test_flag = '20190502050734_05180310_004_01', 'SH_batch02', False
-#track_name, batch_key, test_flag = '20190216200800_07690212_004_01', 'SH_batch02', False
-
-#track_name, batch_key, test_flag = '20190213133330_07190212_004_01', 'SH_batch02', False
-
-# main text figure
-#track_name, batch_key, test_flag = 'SH_20190502_05160312', 'SH_publish', False
-
-#suppl. figures:
-#track_name, batch_key, test_flag = 'SH_20190219_08070210', 'SH_publish', False
 track_name, batch_key, test_flag = 'SH_20190224_08800210', 'SH_publish', False
 
-#print(track_name, batch_key, test_flag)
 hemis, batch = batch_key.split('_')
-#track_name= '20190605061807_10380310_004_01'
 ATlevel= 'ATL03'
 
 
 
-#plot_path   = mconfig['paths']['plot'] + '/'+hemis+'/'+batch_key+'/' + track_name + '/B05_angle/'
 plot_path   = mconfig['paths']['plot'] + '/'+hemis+'/'+batch_key+'/publish/' + track_name + '/'
 MT.mkdirs_r(plot_path)
 bad_track_path =mconfig['paths']['work'] +'bad_tracks/'+ batch_key+'/'
-# %%
 
 all_beams   = mconfig['beams']['all_beams']
 high_beams  = mconfig['beams']['high_beams']
 low_beams   = mconfig['beams']['low_beams']
 beam_groups = mconfig['beams']['groups']
 group_names = mconfig['beams']['group_names']
-#Gfilt   = io.load_pandas_table_dict(track_name + '_B01_regridded', load_path) # rhis is the rar photon data
-
-# load_path   = mconfig['paths']['work'] +'/B01_regrid_'+hemis+'/'
-# G_binned    = io.load_pandas_table_dict(track_name + '_B01_binned' , load_path)  #
 load_path   = mconfig['paths']['work'] +batch_key +'/B02_spectra/'
 Gk          = xr.load_dataset(load_path+ '/B02_'+track_name + '_gFT_k.nc' )  #
 
@@ -92,22 +64,6 @@
 Prior       = MT.load_pandas_table_dict('/A02_'+track_name, load_path)['priors_hindcast']
 
 
-# font_for_print()
-# F = M.figure_axis_xy(5.5, 3, view_scale= 0.8)
-# plt.suptitle(track_name)
-# ax1 =  plt.subplot(2, 1, 1)
-# plt.title('Data in Beam', loc= 'left')
-#
-# xi =1
-
-#data = Marginals.isel(x=xi).sel(beam_group= 'group1').marginals
-# angle_mask = Marginals.angle[2:-2]
-#
-#data.T.plot(cmap= plt.cm.OrRd)
-
-# %%
-
-
 def derive_weights(weights):
     weights = (weights-weights.mean())/weights.std()
     weights = weights - weights.min()
@@ -128,7 +84,6 @@ def weighted_means(data, weights, x_angle, color='k'):
         k           = wi.k.data
         data_k      = data.sel(k=k).squeeze()
         data_weight = (data_k * wi)
-        #plt.stairs(data_weight.sum('k')/ weight_norm , x_angle, linewidth=1 , color ='k')
         if data_k.k.size > 1:
             for k in data_k.k.data:
                 plt.stairs(data_weight.sel(k=k) / weight_norm, x_angle, color ='gray', alpha =0.5)
@@ -187,14 +142,9 @@ def weighted_means(data, weights, x_angle, color='k'):
     tname = track_name.split('_')[1]+'\non '+ track_name.split('_')[0][0:8]
     plt.suptitle('Weighted marginal PDFs for \n$X_i$='+ x_str +' km for track '+tname, y= 1.03, x = 0.125, horizontalalignment= 'left')
 
-    #plt.suptitle('Weighted marginal PDFs\nx='+ x_str +'\n'+track_name, y= 1.05, x = 0.125, horizontalalignment= 'left')
     group_weight = Gweights.isel(x =xi)
 
     ax_list= dict()
-    #ax_sum = F.fig.add_subplot(gs[1, 0])
-    # #ax_sum.tick_params(labelbottom=False)
-    #
-    # ax_list['sum'] = ax_sum
 
     data_collect = dict()
     for group, gpos in zip(Marginals.beam_group.data, [ gs[0, 0], gs[1, 0], gs[2, 0]] ):
@@ -216,41 +166,24 @@ def weighted_means(data, weights, x_angle, color='k'):
         else:
             data_wmean = weighted_means(data, weights, x_angle, color= col_dict[group] )
             plt.stairs(data_wmean , x_angle, color =col_dict[group], alpha =1)
-        # test if density is correct
-        # if np.round(np.trapz(data_wmean) * d_angle, 2) < 0.90:
-        #     raise ValueError('weighted mean is not a density anymore')
 
         if group == 'group1':
-            t_string = group_names[group] +' pair' #group.replace('group',
+            t_string = group_names[group] +' pair'
         else:
-            t_string = group_names[group]+' pair'  #group.replace('group', +' ')
+            t_string = group_names[group]+' pair'
 
         plt.title(next(fn) + t_string, loc ='left')
-        #plt.sca(ax_sum)
-
-        # if data_collect is None:
-        #     data_collect = data_wmean
-        # else:
         data_collect[group] = data_wmean
-        #ax0.set_yscale('log')
 
     data_collect = xr.concat(data_collect.values(), dim='beam_group')
     final_data   = (group_weight * data_collect).sum('beam_group')/group_weight.sum('beam_group').data
 
-    # plt.sca(ax_sum)
-    # plt.stairs( final_data , x_angle, color = 'k', alpha =1, linewidth =0.8)
-    # ax_sum.set_xlabel('Angle (rad)')
-    # plt.title('Weighted mean over group & wavenumber', loc='left')
-
     # get relevant priors
     for axx in ax_list.values():
         axx.set_ylim(0, final_data.max() * 1.5)
-        #figureaxx.set_yscale('log')
         axx.set_xticks(xticks_pi)
         axx.set_xticklabels(xtick_labels_pi)
         axx.set_xlim(-np.pi/2, np.pi/2)
-        #ax_final.set_xticks(xticks_pi)
-        #ax_final.set_xticklabels(xtick_labels_pi)
 
 
     try:
@@ -258,7 +191,6 @@ def weighted_means(data, weights, x_angle, color='k'):
         ax_list['group2'].set_ylabel('PDF')
         ax_list['group3'].set_ylabel('PDF')
         ax_list['group1'].tick_params(labelbottom=True)
-        #ax_list['group3'].set_xlabel('Angle (rad)')
     except:
         pass
 
@@ -272,9 +204,6 @@ def weighted_means(data, weights, x_angle, color='k'):
     plt.stairs( final_data , x_angle, color = 'k', alpha =0.5, linewidth =0.8, zorder= 12)
 
     final_data_smth = lanczos.lanczos_filter_1d(x_angle,final_data, 0.1)
-    #
-    # for group in Marginals.beam_group.data:
-    #     plt.stairs( data_collect.sel(beam_group= group) * group_weight.sel(beam_group= group) /group_weight.sum('beam_group').data, x_angle, color =col_dict[group], alpha =1)
 
     plt.plot(x_angle[0:-1], final_data_smth, color = 'black', linewidth= 0.8)
 
@@ -300,8 +229,6 @@ def weighted_means(data, weights, x_angle, color='k'):
 Gpdf = xr.merge([M_final,M_final_smth])
 
 Gpdf.weighted_angle_PDF_smth.plot()
-#Gpdf.isel( x=slice(0, 3 )).weighted_angle_PDF_smth.mean('x')
-#Gpdf.angle[Gpdf.mean('x').weighted_angle_PDF_smth.argmax()].data
 
 Gpdf.mean('x').weighted_angle_PDF_smth.plot()
 best_guess_angle = Gpdf.angle[Gpdf.mean('x').weighted_angle_PDF_smth.argmax()].data
@@ -310,17 +237,8 @@ def weighted_means(data, weights, x_angle, color='k'):
 
 best_guess_angle/np.pi
 Gpdf.mean('x').weighted_angle_PDF_smth.plot()
-#Gpdf.weighted_angle_PDF.where(~np.isnan(Gpdf.weighted_angle_PDF),0 ).plot()
-
-# if len(Gpdf.x) < 2:
-#     print('not enough x data, exit')
-#     MT.json_save('B05_fail', plot_path+'../',  {'time':time.asctime( time.localtime(time.time()) ) , 'reason': 'not enough x segments'})
-#     print('exit()')
-#     exit()
-
 
-# %%
-class plot_polarspectra(object):
+class plot_polarspectra:
         def __init__(self,k, thetas, data, data_type='fraction' ,lims=None,  verbose=False):
 
             """
@@ -330,9 +248,6 @@ def __init__(self,k, thetas, data, data_type='fraction' ,lims=None,  verbose=Fal
             self.k      =k
             self.data   =data
             self.thetas =thetas
-
-            #self.sample_unit=sample_unit if sample_unit is not None else 'df'
-            # decided on freq limit
             self.lims= lims = [self.k.min(),self.k.max()] if lims is None else lims #1.0 /lims[1], 1.0/ lims[0]
             freq_sel_bool=M.cut_nparray(self.k, lims[0], lims[1] )
 
@@ -348,7 +263,6 @@ def __init__(self,k, thetas, data, data_type='fraction' ,lims=None,  verbose=Fal
                 self.clevs=np.linspace(np.nanpercentile(dir_data.data, 1), np.ceil(self.max* 0.9), 21)
             elif data_type == 'energy':
                 self.ctrs_min=self.min+self.min*.05
-                #self.clevs=np.linspace(self.min, self.max, 21)
                 self.clevs=np.linspace(self.min+self.min*.05, self.max*.60, 21)
 
 
@@ -358,7 +272,6 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
             """
             if ax is None:
                 ax          = plt.subplot(111, polar=True)
-                #self.title  = plt.suptitle('  Polar Spectrum', y=0.95, x=0.5 , horizontalalignment='center')
             else:
                 ax=ax
             ax.set_theta_direction(-1)  #right turned postive
@@ -367,25 +280,18 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
             grid=ax.grid(color='k', alpha=.5, linestyle='-', linewidth=.5)
 
             if self.data_type == 'fraction':
-                cm=plt.cm.RdYlBu_r #brewer2mpl.get_map( 'RdYlBu','Diverging', 4, reverse=True).mpl_colormap
-                colorax = ax.contourf(self.thetas,self.k, self.data, self.clevs, cmap=cm, zorder=1)# ,cmap=cm)#, vmin=self.ctrs_min)
+                cm=plt.cm.RdYlBu_r
+                colorax = ax.contourf(self.thetas,self.k, self.data, self.clevs, cmap=cm, zorder=1)
             elif self.data_type == 'energy':
-                cm=plt.cm.Paired#brewer2mpl.get_map( 'Paired','Qualitative', 8).mpl_colormap
+                cm=plt.cm.Paired
                 cm.set_under='w'
                 cm.set_bad='w'
-                colorax = ax.contourf(self.thetas,self.k, self.data, self.clevs, cmap=cm, zorder=1)#, vmin=self.ctrs_min)
-            #divider = make_axes_locatable(ax)
-            #cax = divider.append_axes("right", size="5%", pad=0.05)
+                colorax = ax.contourf(self.thetas,self.k, self.data, self.clevs, cmap=cm, zorder=1)
 
             if cbar_flag:
                 cbar = plt.colorbar(colorax, fraction=0.046, pad=0.1, orientation="horizontal")
-                # if self.data_type == 'fraction':
-                #     cbar.set_label('Energy Distribution', rotation=0, fontsize=fontsize)
-                # elif self.data_type == 'energy':
-                #     cbar.set_label('Energy Density ('+self.unit+')', rotation=0, fontsize=fontsize)
                 cbar.ax.get_yaxis().labelpad = 30
                 cbar.outline.set_visible(False)
-                #cbar.ticks.
                 clev_tick_names, clev_ticks =MT.tick_formatter(FP.clevs, expt_flag= False, shift= 0, rounder=4, interval=1)
                 cbar.set_ticks(clev_ticks[::5])
                 cbar.set_ticklabels(clev_tick_names[::5])
@@ -397,7 +303,6 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
                 radial_ticks = np.arange(100, 1000, 50)
             print(radial_ticks)
             xx_tick_names, xx_ticks = MT.tick_formatter( radial_ticks , expt_flag= False, shift= 1, rounder=0, interval=1)
-            #xx_tick_names, xx_ticks = MT.tick_formatter( np.arange( np.floor(self.k.min()),self.k.max(), 20) , expt_flag= False, shift= 1, rounder=0, interval=1)
             xx_tick_names = ['  '+str(d)+'m' for d in xx_tick_names]
 
             ax.set_yticks(xx_ticks[::1])
@@ -413,26 +318,22 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
             lines, labels = plt.thetagrids(degrange, labels=degrange_label)#, frac = 1.07)
 
             for line in lines:
-                #L=line.get_xgridlines
                 line.set_linewidth(5)
-                #line.set_linestyle(':')
 
-            #ax.set_yscale('log')
             ax.set_ylim(self.lims)
             ax.spines['polar'].set_color("none")
             ax.set_rlabel_position(87)
             self.ax=ax
 
-# %%
 font_for_print()
 fn = copy.copy(lstrings)
 
 
 F = M.figure_axis_xy(fig_sizes['two_column_square'][0], fig_sizes['two_column_square'][1], view_scale= 0.7, container = True)
-gs = GridSpec(8,6,  wspace=0.1,  hspace=2.1)#figure=fig,
+gs = GridSpec(8,6,  wspace=0.1,  hspace=2.1)
 col.colormaps2(21)
 
-cmap_spec= col.white_base_blgror #plt.cm.ocean_r
+cmap_spec= col.white_base_blgror
 clev_spec = np.linspace(-8, -1, 21) *10
 
 cmap_angle= col.cascade_r
@@ -448,16 +349,12 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 lam_p = 2 *np.pi/weighted_spec.k
 lam = lam_p * np.cos(best_guess_angle)
 k               = 2 * np.pi/lam
-#weighted_spec.k/np.cos(best_guess_angle)
 
 xlims = x_spec[0]-12.5/2, x_spec[-5]
-#weighted_spec.plot()
-#clev_spec = np.linspace(-8, -1, 21) *10
 clev_spec = np.linspace(-80, (10* np.log(weighted_spec)).max() * 0.9, 21)
 
 dd = 10* np.log(weighted_spec.rolling(k=10, min_periods= 1, center=True).mean())
 clev_log = M.clevels( [dd.quantile(0.01).data * 0.3, dd.quantile(0.98).data * 1.5], 31)* 1
-#plt.pcolor(x_spec, k, dd ,vmin= clev_spec[0], vmax= clev_spec[-1],  cmap =cmap_spec )
 plt.pcolormesh(x_spec, lam, dd, cmap=cmap_spec , vmin = clev_log[0], vmax = clev_log[-1])
 
 
@@ -473,19 +370,14 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 cbar = plt.colorbar( fraction=0.018, pad=0.01, orientation="vertical", label ='Power')
 cbar.outline.set_visible(False)
 clev_ticks = np.round(clev_spec[::3], 0)
-#clev_tick_names, clev_ticks =MT.tick_formatter(clev_spec, expt_flag= False, shift= 0, rounder=1, interval=2)
 cbar.set_ticks(clev_ticks)
 cbar.set_ticklabels(clev_ticks)
 
 plt.ylabel('corrected wavelength $(m)$')
-#plt.xlabel('x (km)')
-
-#plt.colorbar()
 ax2 = F.fig.add_subplot(gs[3:5, :])
 ax2.tick_params(labelleft=True)
 
-#Gpdf.weighted_angle_PDF.where(~np.isnan(Gpdf.weighted_angle_PDF),0 ).T.plot()
-dir_data = Gpdf.interp(x= weighted_spec.x).weighted_angle_PDF_smth.T#.rolling(angle=5, min_periods= 1, center=True).mean()
+dir_data = Gpdf.interp(x= weighted_spec.x).weighted_angle_PDF_smth.T
 
 x = Gpdf.x/1e3
 angle = Gpdf.angle
@@ -513,81 +405,49 @@ def linear(self, radial_axis='period', ax=None, cbar_flag=True):
 ax1.set_xticklabels(x_tick_labels)
 ax2.set_xticklabels(x_tick_labels)
 
-#ax1.set_yscale('log')
 lam_lim= lam[-1].data, lam[10].data
 ax1.set_ylim(lam_lim)
 
 ax1.set_xlim(xlims)
 ax2.set_xlim(xlims)
-#ax2.set_yscale('log')
 ax2.axhline(best_guess_angle, color=col.orange, linewidth=0.8)
 
 
-#xx_list = np.insert(weighted_spec.x.data, 0, 0)
-# x_pos_list = spec.create_chunk_boundaries( 1,  xx_list.size,  iter_flag= False )
-# #x_pos_list = spec.create_chunk_boundaries( int(xx_list.size/3),  xx_list.size,  iter_flag= False )
-# x_pos_list = x_pos_list[:, ::2]
-# x_pos_list[-1, -1] = xx_list.size-1re
-#x_pos_list#.shape
-
-x_pos_list =  [0, 1, 2]#np.arange(0,9, 1)#np.vstack([np.arange(1,3), np.arange(0,3)+1])
+x_pos_list =  [0, 1, 2]
 #x_pos_list
 lsrtrings = iter(['c)', 'd)', 'e)'])
 
 for x_pos, gs in zip( x_pos_list , [ gs[-3:, 0:2], gs[-3:, 2:4], gs[-3:, 4:]] ):
-    #print( x_pos)
-    #print( xx_list[x_pos])
-    x_range = weighted_spec.x.data[x_pos]#, x_pos[-1]]]
+    x_range = weighted_spec.x.data[x_pos]
     print(x_range)
     ax1.axvline(x_range/1e3, linestyle= '-', color= col.green, linewidth=0.9, alpha = 0.8)
     ax2.axvline(x_range/1e3, linestyle= '-', color= col.green, linewidth=0.9, alpha = 0.8)
 
     i_lstring = next(lsrtrings)
     ax1.text(x_range/1e3, np.array(lam_lim).mean()*3/2, ' '+ i_lstring, fontsize= 8, color =col.green)
-    #ax2.text(x_range/1e3, weighted_spec.k.mean().data, ' a', fontsize= 8)
-
-
-    # ax1.axvline(x_range[-1]/1e3, color = 'gray', alpha = 0.5)
-    #
-    # ax2.axvline(x_range[0]/1e3, linestyle= ':', color= 'white', alpha = 0.5)
-    # ax2.axvline(x_range[-1]/1e3, color = 'gray', alpha = 0.5)
 
-
-    # i_spec  = weighted_spec.sel(x= slice(x_range[0], x_range[-1]) )
-    # i_dir   = corrected_marginals.sel(x= slice(x_range[0], x_range[-1]) )
     i_spec  = weighted_spec.isel(x= x_pos )
     i_dir   = corrected_marginals.interp(x= weighted_spec.x).isel(x= x_pos )
     print(i_spec.x.data, i_spec.x.data)
     dir_data  = (i_dir * i_dir.N_data).sum([ 'beam_group'])/ i_dir.N_data.sum([ 'beam_group'])
     lims = dir_data.k[ (dir_data.sum('angle')!=0) ][0].data, dir_data.k[ (dir_data.sum('angle')!=0)  ][-1].data
 
-    #dir_data.plot()
-    #dir_data.rolling(angle =5,  min_periods= 1, center=True ).mean().plot()
-
     N_angle = i_dir.angle.size
-    dir_data2 =  dir_data#.where( dir_data.sum('angle') !=0, 1/N_angle/d_angle )
+    dir_data2 =  dir_data
 
-    plot_data  = dir_data2  * i_spec#.mean('x')
-    plot_data  = dir_data2.rolling(angle =2, k =15, min_periods= 1, center=True ).mean()  * i_spec#.mean('x')
+    plot_data  = dir_data2  * i_spec
+    plot_data  = dir_data2.rolling(angle =2, k =15, min_periods= 1, center=True ).mean()  * i_spec
 
     plot_data = plot_data.sel(k=slice(lims[0],lims[-1] ) )
     xx = 2 * np.pi/plot_data.k
 
-    #F = M.figure_axis_xy(5, 4)
-    #ax = plt.subplot(1, 1, 1, polar=True)
-    #
     if np.nanmax(plot_data.data) != np.nanmin(plot_data.data):
 
         ax3 = F.fig.add_subplot(gs, polar=True)
         FP= plot_polarspectra(xx, plot_data.angle, plot_data, lims=[xx[-1], 340 ] , verbose= False, data_type= 'fraction')
         FP.clevs=np.linspace(np.nanpercentile(plot_data.data, 1), np.round(plot_data.max(), 4), 21)
         FP.linear(ax = ax3, cbar_flag=False)
-        #FP.cbar.set_label('Energy Density ( (m/m)$^2$ k$^{-1}$ deg$^{-1}$ )', rotation=0, fontsize=10)
-        #plt.show()
         plt.title('\n\n'+i_lstring,y=1.0, pad=-6, color=col.green)
 
 F.save_pup(path = plot_path, name = 'B05_dir_ov_'+track_name)
 F.save_light(path = plot_path, name = 'B05_dir_ov_'+track_name)
-# %%
-#F.save_pup(path = plot_path, name = 'B05_dir_ov_'+track_name)
-# MT.json_save('B05_success', plot_path + '../', {'time':time.asctime( time.localtime(time.time()) )})
diff --git a/analysis/SB04_2d_wavefield_emulator.py b/analysis/SB04_2d_wavefield_emulator.py
index ccb030c3..7f2a2037 100644
--- a/analysis/SB04_2d_wavefield_emulator.py
+++ b/analysis/SB04_2d_wavefield_emulator.py
@@ -1,6 +1,5 @@
 
 import os, sys
-#execfile(os.environ['PYTHONSTARTUP'])
 
 """
 This file open a ICEsat2 track applied filters and corections and returns smoothed photon heights on a regular grid in an .nc file.
@@ -10,8 +9,6 @@
     exec(open(os.environ['PYTHONSTARTUP']).read())
     exec(open(STARTUP_2021_IceSAT2).read())
 
-    #%matplotlib inline
-
     import ICEsat2_SI_tools.convert_GPS_time as cGPS
     import h5py
     import ICEsat2_SI_tools.io as io
@@ -22,12 +19,9 @@
     import spicke_remover
     import datetime
     import concurrent.futures as futures
-    #import s3fs
-    # %%
+
     track_name, batch_key, test_flag = io.init_from_input(sys.argv) # loads standard experiment
-    #track_name, batch_key, test_flag = '20190605061807_10380310_004_01', 'SH_batch01', False
-    #track_name, batch_key, test_flag = '20190601094826_09790312_004_01', 'SH_batch01', False
-    #track_name, batch_key, test_flag = '20190207111114_06260210_004_01', 'SH_batch02', False
+
     track_name, batch_key, test_flag = '20190219073735_08070210_004_01', 'SH_batch02', False
 
 
@@ -35,23 +29,21 @@
 
     #print(track_name, batch_key, test_flag)
     hemis, batch = batch_key.split('_')
-    #track_name= '20190605061807_10380310_004_01'
     ATlevel= 'ATL03'
 
     save_path   = mconfig['paths']['work'] + '/B03_spectra_'+hemis+'/'
     save_name   = 'B03_'+track_name
 
-    #plot_path   = mconfig['paths']['plot'] + '/'+hemis+'/'+batch_key+'/' + track_name + '/B_spectra/'
     plot_path   = mconfig['paths']['plot'] + '/phase_fitting_fake/2D_fake/'
     MT.mkdirs_r(plot_path)
     MT.mkdirs_r(save_path)
     bad_track_path =mconfig['paths']['work'] +'bad_tracks/'+ batch_key+'/'
-# %%
+
 if __name__ == '__main__':
     all_beams   = mconfig['beams']['all_beams']
     high_beams  = mconfig['beams']['high_beams']
     low_beams   = mconfig['beams']['low_beams']
-    #Gfilt   = io.load_pandas_table_dict(track_name + '_B01_regridded', load_path) # rhis is the rar photon data
+
 
     load_path   = mconfig['paths']['work'] +'/B01_regrid_'+hemis+'/'
     Gd      = io.load_pandas_table_dict(track_name + '_B01_binned' , load_path)  #
@@ -70,30 +62,21 @@
     Gspec.coords['T'] = 1/Gspec.coords['f']
     Gspec=Gspec.swap_dims({'k': 'f'})
 
-    #Gspec.spectral_power_optm.sel(beam='weighted_mean').plot()
-    # %%
-    #k_lim= 0.02
+
     A, B = Gspec.sel(beam= 'gt2r').Y_model_hat , Gspec.sel(beam= 'gt2l').Y_model_hat
 
     r_ave_kargs={'x':2, 'f':10, 'center':True, 'min_periods':2}
     r_ave_kargs2={'f':10, 'center':True, 'min_periods':2}
-    #(abs(B) - abs(A)).plot()
+
     S_aa = (A*A.conj()).real
     S_bb = (B*B.conj()).real
-    # co_spec = (A.conj() *B) /S_aa/S_bb
-    # abs(co_spec).plot()
+
     co_spec = (A.conj() *B).rolling(**r_ave_kargs).mean()
     np.log(abs(co_spec)).plot(levels=np.arange(-2, 3, 0.1))
 
     (abs(co_spec)).plot(levels=np.exp(np.arange(-3, 2, 0.1)))
 
     abs(co_spec).mean('x').plot()
-
-#(abs(co_spec)/(S_aa *S_bb).rolling(**r_ave_kargs).mean()).plot()
-#
-# (abs(A.conj() *B)/(S_aa *S_bb)).rolling(**r_ave_kargs).mean()[:,:].plot()
-
-# %%
 if __name__ == '__main__':
     L1 = 50
     k1 = 2* np.pi /L1
@@ -111,11 +94,9 @@
     alpha = 35
     kk, ll = np.cos(alpha * np.pi/180) * np.array([0.9*k1, k1, 1.1* k1]), np.sin(alpha * np.pi/180) * np.array([0.9* k1, 1*k1, 1.1* k1])
     M_k, M_l = kk.size, ll.size
-    #y =np.sin(k1* x) + np.sin(k2* x)
     kk_mesh, ll_mesh = np.meshgrid(kk, ll)
     kk_mesh, ll_mesh = kk_mesh.reshape(kk_mesh.size), ll_mesh.reshape(ll_mesh.size)
-    G = np.cos(np.outer(XX, kk_mesh) + np.outer(YY, ll_mesh)).T# + np.sin(np.outer(XX, kk_mesh) + np.outer(YY, ll_mesh)).T
-    #G = np.vstack([ np.cos(np.outer(x, k) + np.outer(y, l)).T ,  np.sin(np.outer(x, k)  + np.outer(y, l) ).T ] ).T
+    G = np.cos(np.outer(XX, kk_mesh) + np.outer(YY, ll_mesh)).T
     G.shape
 
     plt.contourf(x, y, G.sum(0).reshape(Ny, Nx) )
@@ -124,7 +105,6 @@
  # %% radial coordincates
 
 def gaus_2d(x, y, pos_tuple, sigma_g ):
-    #grid = ( (XX - pos_tuple[0])  * (YY - pos_tuple[1]) )
     gx = np.exp(-0.5 * (x - pos_tuple[0])**2 /sigma_g**2  )
     gy = np.exp(-0.5 * (y - pos_tuple[1])**2 /sigma_g**2  )
     return  np.outer(gx ,  gy).T
@@ -140,8 +120,6 @@ def gaus_2d(x, y, pos_tuple, sigma_g ):
     plt.contourf(k_range, l_range, gaus_lk  )
     plt.axis('equal')
 
-# %%
-
 k_0 = 0.03
 l_0 = 0
 dk = 0.01
@@ -203,7 +181,6 @@ def gaus_2d_mesh(XX,YY, pos_tuple, sigma_g ):
     plt.axis('equal')
 
 
-# %% radial coodinates
 def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_flag = True, amp_std= None, random=True):
     """
     inputs:
@@ -263,10 +240,6 @@ def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_fl
 
 if __name__ == '__main__':
     font_for_pres()
-    #k_mesh.shape
-    #for angle in np.arange(-80, 80+20, 40):
-    #for phase in np.arange(0, 2*np.pi, np.pi/3):
-    #for k_abs in np.arange(0.01, 0.09, 0.01):
 
     angle =45
     amp =1
@@ -276,24 +249,16 @@ def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_fl
     for dk in np.arange(0.005, 0.03, 0.002):
 
         for size in [2, 5, 10, 30, 50, 100, 200]:
-            #for angle in np.arange(-80, 80+20, 40):
-            #for k_abs in np.arange(0.01, 0.09, 0.01):
-            #for phase in np.arange(0, 2*np.pi, np.pi/3):
-
-
             F = M.figure_axis_xy(8, 3, view_scale = 0.5, container =False)
             plt.suptitle('k_abs=' + str(k_abs) +' angle=' + str(angle) + ' size=' + str(size) +' dk=' + str(dk) )
             ax = plt.subplot(1, 2, 1)
             k_list, l_list, amp_weights, stancil_shape = get_stancils_polar(k_abs, angle * np.pi/180, size=size, dk = dk, mesh = True , plot_flag= True, random = True)
-
             circle1 = plt.Circle((k_list[0], l_list[0]), dk, color='b', fill=False)
             ax.add_patch(circle1)
 
-
             k_noise, l_noise, amp_noise, stancil_shape = get_stancils_polar(0.8, 0 * np.pi/180, size=20, dk = 0.3, mesh = True , plot_flag= True, random = True)
             amp_noise = (amp_noise *0+1) * 0
 
-
             plt.xlim(0, 0.1)
             #plt.axis('equal')
             plt.ylim(-0.1, 0.1)
@@ -308,10 +273,6 @@ def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_fl
             amp_all = np.concatenate([amp_weights, amp_noise])
             amp_all.shape
             G = np.vstack([ np.cos(np.outer(XX, k_all) + np.outer(YY, l_all)).T ,  np.sin(np.outer(XX, k_all) + np.outer(YY, l_all)).T ] ).T
-            #G = np.vstack([ np.cos(np.outer(x, k) + np.outer(y, l)).T             ,  np.sin(np.outer(x, k)  + np.outer(y, l) ).T ] ).T
-
-            #phase1 = np.random.rand(1, amp_list.size) *  np.pi*2
-            #phase = np.arange(0, amp_list.size) *  np.pi/2
 
             b = np.hstack([ np.cos(phase)*amp_all, np.sin(phase) *amp_all]).squeeze() * amp
             z_model = (G @ b).reshape(Ny, Nx)
@@ -323,9 +284,3 @@ def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_fl
             F.save_light(path = plot_path, name = 'fake_2d_dk' +str(dk) +'_s' + str(int(size)))
             plt.show()
 
-
-    # class twoD_wave_packets(object):
-    #     def __init__(self):
-
-
-    # %%
diff --git a/analysis_db/B05_define_angle.py b/analysis_db/B05_define_angle.py
index 7cf9f5ff..05aa1c37 100644
--- a/analysis_db/B05_define_angle.py
+++ b/analysis_db/B05_define_angle.py
@@ -280,7 +280,7 @@ def weighted_means(data, weights, x_angle, color='k'):
     exit()
 
 # %%
-class plot_polarspectra(object):
+class plot_polarspectra:
         def __init__(self,k, thetas, data, data_type='fraction' ,lims=None,  verbose=False):
 
             """
diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/angle_optimizer.py b/src/icesat2_tracks/ICEsat2_SI_tools/angle_optimizer.py
index ccb858fe..fa167fcc 100644
--- a/src/icesat2_tracks/ICEsat2_SI_tools/angle_optimizer.py
+++ b/src/icesat2_tracks/ICEsat2_SI_tools/angle_optimizer.py
@@ -122,7 +122,7 @@ def simple_log_panelty(x, x0, sigma):
         return  - 0.5 * (cost_sqrt/tot_var + np.log(tot_var) ).sum() + prior_weight * penalties
 
 
-class sample_with_mcmc(object):
+class sample_with_mcmc:
     """
     sample a 2nd surface using mcmc and other methods. its made for getting a quick estimate!
 
diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py b/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py
index 2ca2ca4d..72eda331 100644
--- a/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py
+++ b/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py
@@ -156,7 +156,7 @@ def define_weights(stancil, prior, x, y, dx, k, max_nfev, plot_flag=False):
     return weight, prior_pars
 
 
-class wavenumber_spectrogram_gFT(object):
+class wavenumber_spectrogram_gFT:
     def __init__(self, x, data, L, dx, wavenumber, data_error=None, ov=None):
         """
         returns a wavenumber spectrogram with the resolution L-ov
@@ -720,7 +720,7 @@ def power_from_model(p_hat, dk, M, N_x, N_x_full):
     return spec
 
 
-class generalized_Fourier(object):
+class generalized_Fourier:
     def __init__(self, x, ydata, k):
         """
         non_dimensionalize (bool, default=True) if True, then the data and R_data_uncertainty is non-dimensionalized by the std of the data
diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/io.py b/src/icesat2_tracks/ICEsat2_SI_tools/io.py
index 332ce7f6..521e305c 100644
--- a/src/icesat2_tracks/ICEsat2_SI_tools/io.py
+++ b/src/icesat2_tracks/ICEsat2_SI_tools/io.py
@@ -70,13 +70,10 @@ def ID_to_str(ID_name):
     date
     return IDs[0] +' ' +date +' granule: ' + IDs[2]
 
-class case_ID(object):
+class case_ID:
     """docstring for case_ID"""
     def __init__(self, track_name):
         import re
-        super(case_ID, self).__init__()
-
-        #track_name_pattern = r'(\d{4})(\d{2})(\d{2})(\d{2})(\d{2})(\d{2})_(\d{4})(\d{2})(\d{2})_(\d{3})_(\d{2})'
         track_name_pattern = r'(\D{2}|\d{2})_?(\d{4})(\d{2})(\d{2})(\d{2})?(\d{2})?(\d{2})?_(\d{4})(\d{2})(\d{2})_?(\d{3})?_?(\d{2})?'
         case_ID_pattern = r'(\d{4})(\d{2})(\d{2})_(\d{4})(\d{2})(\d{2})'
 
diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/spectral_estimates.py b/src/icesat2_tracks/ICEsat2_SI_tools/spectral_estimates.py
index 52960e44..d46d9ea2 100644
--- a/src/icesat2_tracks/ICEsat2_SI_tools/spectral_estimates.py
+++ b/src/icesat2_tracks/ICEsat2_SI_tools/spectral_estimates.py
@@ -376,7 +376,7 @@ def get_lon_lat_coords(stancil):
     return coord_positions
 
 
-class wavenumber_spectrogram(object):
+class wavenumber_spectrogram:
     def __init__(self, x_grid, data, Lpoints, ov=None, window=None):
         """
         returns a wavenumber spectrogram with the resolution L-ov
@@ -512,7 +512,7 @@ def get_stancil_var_apply(stancil):
             self.G.attrs['mean_variance_detrended_chunks']  = np.array(stancil_vars).mean()
             self.G.attrs['mean_variance_pwelch_spectrum']   = self.calc_var()
 
-class wavenumber_spectrogram_LS_even(object):
+class wavenumber_spectrogram_LS_even:
     def __init__(self, x, data, L, waven_method = 'fftX2' , dy=None ,  ov=None, window=None, kjumps=1):
         """
         returns a wavenumber spectrogram with the resolution L-ov
@@ -618,7 +618,7 @@ def parceval(self, add_attrs=True ):
     def mean_spectral_error(self, confidence = 0.95 ):
         return wavenumber_spectrogram.mean_spectral_error(self, confidence= confidence )
 
-class wavenumber_spectrogram_LS(object):
+class wavenumber_spectrogram_LS:
     def __init__(self, x, data, L, dx, dy = None, waven_method = 'fftX2', ov=None, window=None):
         """
         returns a wavenumber spectrogram with the resolution L-ov
@@ -925,7 +925,7 @@ def mean_spectral_error(self, mask=None, confidence = 0.95 ):
 
 
 # class for getting standard Pwelch spectrum. old version, deprechiate
-class wavenumber_pwelch(object):
+class wavenumber_pwelch:
     def __init__(self,data, x, L, ov=None, window=None, save_chunks=False, plot_chunks=False):
         """
         returns a wavenumber spectrum using the pwelch method
@@ -1098,7 +1098,7 @@ def calc_var(self):
 
 # %% optimze spectral variance
 
-class conserve_variance(object):
+class conserve_variance:
     def __init__(self,Z, freq, data, nan_mask= None):
 
         """
diff --git a/src/icesat2_tracks/analysis_fake_data/X03_MCMC_surface_smapling.py b/src/icesat2_tracks/analysis_fake_data/X03_MCMC_surface_smapling.py
index 312c157d..4e3f65a9 100644
--- a/src/icesat2_tracks/analysis_fake_data/X03_MCMC_surface_smapling.py
+++ b/src/icesat2_tracks/analysis_fake_data/X03_MCMC_surface_smapling.py
@@ -41,7 +41,7 @@ def cost(x, y):
     z = 4+ np.sin(4*  2 * np.pi *x/Lx) +  np.sin( 3 * np.pi *x/Lx - np.pi/5) + np.cos(1*  2 * np.pi *y/Ly) +  np.sin( 3 * np.pi *y/Ly - np.pi/3)
     return z**2
 
-class sample_with_mcmc(object):
+class sample_with_mcmc:
     """
     sample a 2nd surface using mcmc. its make for getting a quick estimate!
     """
diff --git a/src/icesat2_tracks/local_modules/m_colormanager_ph3.py b/src/icesat2_tracks/local_modules/m_colormanager_ph3.py
index 074e07a7..17163538 100644
--- a/src/icesat2_tracks/local_modules/m_colormanager_ph3.py
+++ b/src/icesat2_tracks/local_modules/m_colormanager_ph3.py
@@ -50,7 +50,7 @@ def json_load(name, path, verbose=False):
 
 
 
-class color(object):
+class color:
         def __init__(self, path=None, name=None):
             self.white=(1,1,1)
             if (path is not None) & (name is not None):
diff --git a/src/icesat2_tracks/local_modules/m_general_ph3.py b/src/icesat2_tracks/local_modules/m_general_ph3.py
index e268658b..6e24a8d9 100644
--- a/src/icesat2_tracks/local_modules/m_general_ph3.py
+++ b/src/icesat2_tracks/local_modules/m_general_ph3.py
@@ -27,7 +27,7 @@
 import matplotlib.pyplot as plt
 
 
-class color(object):
+class color:
         def __init__(self):
 
             self.red=(203/255, 32/255, 39/255)
@@ -61,7 +61,7 @@ def show(self):
 
 # funny massage
 
-class figure_axis_xy(object):
+class figure_axis_xy:
         """define standart  XY Plot with reduced grafics"""
 
         def __init__(self,x_size=None,y_size=None,view_scale=None, size_tuple=None , fig_scale=None, container=False, dpi=180):
@@ -171,7 +171,7 @@ class subplot_routines(figure_axis_xy):
         def __init__(self, ax):
             self.ax=ax
 
-class plot_sprecta(object):
+class plot_sprecta:
             def __init__(self,fs, Xdata,sample_unit=None,data_unit=None):
 
                 self.fs=fs
@@ -261,7 +261,7 @@ def power(self, Color='b', fig_scale=2, fax='f'):
                 self.F.make_clear()
                 plt.grid()
 
-class plot_periodogram(object):
+class plot_periodogram:
             def __init__(self,time,fs, data,clevs=None, sample_unit=None, data_unit=None,
                         ylim=None, time_unit=None, cmap=None):
                 self.fs=fs[1:]
@@ -696,7 +696,7 @@ def set_xaxis_to_days(self, **kwargs):
                 set_timeaxis_days(self.ax, **kwargs)
 
 
-class plot_polarspectra(object):
+class plot_polarspectra:
         def __init__(self,f, thetas, data,unit=None, data_type='fraction' ,lims=None,  verbose=False):
 
             self.f=f
@@ -1256,7 +1256,7 @@ def spickes_to_mean(ts, nloop=None, spreed=1, gaussian=True):
 
 ## Composites
 
-class composite_data(object):
+class composite_data:
         def __init__(self,var, index_weight=None):
             #print(var.shape)
             self.composites=var
@@ -1285,7 +1285,7 @@ def bootstrap(self, ci=[2.5, 50, 97.5], reps=1000,):
             yb = 1/np.arange(1, n+1)[:, None] * np.cumsum(xb, axis=0)
             upper, lower = np.percentile(yb, [2.5, 97.5], axis=1)
 
-class comp_iter(object):
+class comp_iter:
         def __init__(self, span, dt=None, unit=None):
             self.span=list(span)
             for s in self.span:
@@ -1306,7 +1306,7 @@ def __init__(self, span, dt=None, unit=None):
 
                 self.time_iter_string=time_str
 
-class composite(object):
+class composite:
         def __init__(self,index, time=None, weigthing=False, span=None):
             """ Initial Class for bulding composite based on:
                 index       position in the time vector 'time'
diff --git a/src/icesat2_tracks/local_modules/m_spectrum_ph3.py b/src/icesat2_tracks/local_modules/m_spectrum_ph3.py
index 19b49700..f904b63c 100644
--- a/src/icesat2_tracks/local_modules/m_spectrum_ph3.py
+++ b/src/icesat2_tracks/local_modules/m_spectrum_ph3.py
@@ -91,7 +91,7 @@ def spicke_remover(data, nstd=20.0, spreed=500.0, max_loops=10.0 , verbose=False
 
     return data2 , act_flag
 
-class Spectrum(object):
+class Spectrum:
     """ A class that represents a single realization of
             the one-dimensional spectrum  of a given field phi """
 
@@ -149,7 +149,7 @@ def parceval(self):
 
 
 
-class moments(object):
+class moments:
     def __init__(self,data_org,dt, L=None, ov=None,window=None, save_chunks=False, plot_chunks=False, prewhite=None):
         """
         This function calculates the spectral moments from a station (buoy, GPS, or seismic station) that measures
@@ -396,7 +396,7 @@ def cal_MEM(self, theta=None, flim=(0.01, .5)):
 
 
 
-class pwelch(object):
+class pwelch:
     def __init__(self,data,dt,L=None, ov=None,window=None, save_chunks=False, plot_chunks=False, periodogram=False, prewhite=None):
         """
         prewhite    None(default)
@@ -800,7 +800,7 @@ def save_data(self, path=None, S=None):
         save_file(P, path)
 
 
-class save_data_periodogram(object):
+class save_data_periodogram:
     def __init__(self,P, S=None):
         self.meta=S.meta if S is not None else ''
         self.data_unit=S.unit if S is not None else ''