# -*- coding: utf-8 -*-
"""
Created on Thu Jan 28 14:52:44 2021

@author: ppantina
"""

##Import the required libraries
import glob
import sys
import os
import h5py
import datetime
import numpy          as     np
from matplotlib       import pyplot as plt
from scipy            import ndimage  #for convolutions
from silx.io.dictdump import dicttoh5 #for hdf writing       
from datetime         import timezone
import L1B_sub_add_horiz_wind_hrrr_fullflight    #to calculate winds
import L1B_sub_along_cross_winds2     #to calculate along/cross winds
import L1B_sub_plotter
import sub_cmap        #for colorbar
import sub_params      #for parameters
import sub_flightLines #for flightLines
import L1B_sub_EXRAD_SCAN
##Change these variables as needed
radName       = input('Which radar? CRS, EXRAD, HIWRAP_KU, HIWRAP_KA: ')
flightDate    = input('Which date? yyyymmdd (append flight number if needed) ')
versionLetter = input('Which version?: A, B, C, ... ')
makePlots     = input('Run the plots?: yes, no, ... ')


##Next, read the parameter file to fine radar-relevant variables.
fileNameParams = 'radar_params_irma.xlsx'
config = sub_params.sub_params(fileNameParams, flightDate, radName)

##Define the colorbar
cmap_n0q = sub_cmap.make_cmap('n0q', bit = True)

######Define paths and filenames######
##First locate L1A paths/files
L1A_dir  = config['Directory_L1A']
if radName== 'EXRAD_SCAN':##<------------------------------------new, the files have been renamed to just EXRAD
    L1A_file = np.sort(glob.glob(L1A_dir + '*' + 'EXRAD' + '*' + flightDate + '*.h5'))
    L1A      = h5py.File(L1A_file[0], 'r')

elif radName== 'EXRAD_NADIR' :##<------------------------------------new, the files have been renamed to just EXRAD
    L1A_file = glob.glob(L1A_dir + '*' + 'EXRAD' + '*' + flightDate + '*.h5')[0]
    print(L1A_file)
    L1A      = h5py.File(L1A_file, 'r')
else:
    L1A_file = glob.glob(L1A_dir + '*' + radName + '*' + flightDate + '*.h5')[0]
    print(L1A_file)
    L1A      = h5py.File(L1A_file, 'r')

##Next locate L1B paths/files. Create the output dir if necessary.
L1B_dir      = config['Directory_L1B']
L1B_file     = config['experimentName'] + '_' + radName + '_L1B_Rev'+ versionLetter + '_' + flightDate + '.h5'
L1B_filepath = L1B_dir + L1B_file
L1B_plotpath = L1B_dir + 'plots/' + radName + '/' + flightDate + '/'
if not(os.path.isdir(os.path.dirname(L1B_filepath))): os.makedirs(os.path.dirname(L1B_filepath)) ##make the output dir if needed.
if not(os.path.isdir(os.path.dirname(L1B_plotpath))): os.makedirs(os.path.dirname(L1B_plotpath)) ##make the output dir if needed.

##Finally locate the HRRR filepath
m_filepath   = config['Directory_hrrr'] + '/nc3d/*.nc'

##Define some vars to index the ranges.
rangeMin    = config['rangeMin'   ]
rangeMax    = config['rangeMax'   ]
deg_offset  = config['degOffset'  ] ##roll_offset for wind correction.
##calCoeffs defined separately. For KA, there is a second offset.
##this is a roundabout way to add tuples.
extCalCoeff = np.array(eval(config['extCalCoeff'])) + np.array(eval(config['extCalCoeff2'])) ##<---------------------------------------------------------------------------------------------

##Set up range indices and calibration coefficient arrays.
if radName == 'CRS':
    rangeIndex  = np.where((L1A['Products/CoPolChirp/Information/Range'][0]>rangeMin) & (L1A['Products/CoPolChirp/Information/Range'][0]<rangeMax))[0]#<-------------------------------------------------------------------------------------------------------
if radName == 'EXRAD_NADIR':
    rangeIndex  = np.where((L1A['Products/Nadir/Information/Range'][0]>rangeMin) & (L1A['Products/Nadir/Information/Range'][0]<rangeMax))[0] ##<-------------------------------
    
##Read in the appropriate flight lines.
flightLines = sub_flightLines.sub_flightLines (flightDate)


#######################Add /Information to h5#######################
##To write to dict, all strings must be placed in arrays, as far as I can tell.
print ('Adding Information')
L1B                = {}
L1B['Information'] = {}
L1B['Information/ExperimentName'   ] = ['IMPACTS2023']#L1A['ExperimentName' ]<--------------------------
L1B['Information/L1A_ProcessDate'  ] = ['']#L1A['L1A_ProcessDate']<--------------------------
L1B['Information/FlightDate'       ] = [flightDate          ]
L1B['Information/RadarName'        ] = [radName             ]
L1B['Information/L1B_ProcessDate'  ] = [datetime.datetime.now().strftime("%Y%m%d")]
L1B['Information/L1B_Revision'     ] = ['Rev' + versionLetter]
L1B['Information/DataContact'      ] = [config['dataContact' ]]
L1B['Information/MissionPI'        ] = [config['missionPI'   ]]
L1B['Information/Aircraft'         ] = [config['aircraft'    ]]
L1B['Information/InstrumentPI'     ] = [config['instrumentPI']]
L1B['Information/L1B_Revision_Note'] = [config['revisionNote']]

if radName == 'EXRAD_SCAN':
    L1B_ex_last = L1B_sub_EXRAD_SCAN.L1B_sub_EXRAD_SCAN (flightDate, L1A_file, rangeMin, rangeMax, flightLines, extCalCoeff, versionLetter, L1B, L1B_dir, cmap_n0q, makePlots)
    print ('SCAN finishes by exiting the script early')
    sys.exit()
##For HIWRAP, save data every 0.5 seconds instead of every 0.25. It is too large otherwise. For other radars, save all times.
if 'HIWRAP' in radName: times2save = np.where((L1A['Time/Data/TimeUTC'][:]%1 == 0) | (L1A['Time/Data/TimeUTC'][:]%1 == 0.5))[0]
else:                   times2save = np.where( L1A['Time/Data/TimeUTC'][:] > 0)[0] 

#######################Add Time/Data and Time/Information#######################
print ('Adding Time')
flightYear = flightDate[0:4] #extract the year from flightDate. Use it for printing below.
L1B['Time']             = {}
L1B['Time/Data']        = {}
L1B['Time/Information'] = {}
L1B['Time/Data/TimeUTC'                          ] = L1A['Time/Data/TimeUTC'][:][times2save].T[0]#<---------------------------------------- ##Select only the relevant times
L1B['Time/Information/TimeUTC_units'             ] = ['seconds']
L1B['Time/Information/TimeUTC_description'       ] = L1A['Time/Information/Description']
L1B['Time/Information/TimeUTC_01Jan' + flightYear] = [datetime.datetime(int(flightYear),1,1).replace(tzinfo=timezone.utc).timestamp()]
L1B['Time/Information/TimeUTC_01Jan' + flightYear + '_description'] = ['Time of 0 UTC, Jan 01, ' + flightYear + ', for reference if the user does not have an easy linux time converter']

#######################Add Products/Information#######################
print ('Adding Product Information')
L1B['Products'            ] = {}
L1B['Products/Data'       ] = {}
L1B['Products/Information'] = {}

##Redefine the ranges (channels are different between radars), then calculate GateSpacing
if radName == 'CRS'          :L1B['Products/Information/Range'] = L1A['Products/CoPolChirp/Information/Range'][:][0][rangeIndex]#<---------------------------------
if radName == 'EXRAD_NADIR'  :L1B['Products/Information/Range'] = L1A['Products/Nadir/Information/Range'][:][0][rangeIndex] ##<-------------------------------
ranges = np.copy(L1B['Products/Information/Range'])            ##Use this as a shortcut later
L1B['Products/Information/GateSpacing'] = np.diff(ranges[0:2]) ##Just the difference between first two gates

##Populate aircraft motion. Use the smoothed motion where available.
if radName == 'CRS': ##no smoothing of AircraftMotion
    L1B['Products/Information/AircraftMotion'            ] = L1A['Products/Information/AircraftMotion'][:].T[0]#<------------------
    L1B['Products/Information/AircraftMotion_units'      ] = ['m/s']
    L1B['Products/Information/AircraftMotion_description'] = ['Estimated aircraft motion subtracted from Doppler estimate']
if radName == 'EXRAD_NADIR':
    L1B['Products/Information/AircraftMotion'            ] = L1A['Products/Information/AircraftMotion'][:]#<------------------
    L1B['Products/Information/AircraftMotion_units'      ] = ['m/s']
    #L1B['Products/Information/AircraftMotion_description'] = ['Estimated aircraft motion normal to the beam subtracted from Doppler estimate. Smoothed to a 2-second average motion.']
    L1B['Products/Information/AircraftMotion_description'] = ['Estimated aircraft motion normal to the beam subtracted from Doppler estimate.']

##Added turn flag Dec 2022
turnFlag = np.zeros_like(L1B['Time/Data/TimeUTC'][:], dtype = np.int8)
for ff in flightLines:
    idx = np.where((L1B['Time/Data/TimeUTC'][:]>ff[0])&(L1B['Time/Data/TimeUTC'][:]<ff[1]))
    turnFlag[idx] = 1
L1B['Products/Information/AircraftTurnFlag'            ] = turnFlag
L1B['Products/Information/AircraftTurnFlag_units'      ] = ['0 or 1']
L1B['Products/Information/AircraftTurnFlag_description'] = ['Flag is 1 when plane is flying level without turns. Flag is 0 elsewhere.']

##These are common to all radars
L1B['Products/Information/PRI'               ] = [str(int(config['L1B_HighPRF']*1e6)) + ' us / '+str(int(config['L1B_LowPRF']*1e6))+' us staggered']
L1B['Products/Information/AveragedPulses'    ] = np.int32  ([config['L1B_averagedPulses'       ]])
L1B['Products/Information/AntennaBeamwidth'  ] = np.float32([config['L1B_beamWidth'            ]])
L1B['Products/Information/AntennaSize'       ] = np.float32([config['AntennaSize'              ]])
if 'EXRAD' in radName: #<-----------------
    L1B['Products/Information/Frequency'         ] = L1A['Products/Nadir/Information/Frequency'][:] #<-----------------
else:
    L1B['Products/Information/Frequency'         ] = L1A['Products/CoPolChirp/Information/Frequency'][:] #<-----------------
L1B['Products/Information/Wavelength'        ] = 3e8/L1B['Products/Information/Frequency'     ][:]

##Populate pulses, frequency, wavelength units/descriptions
L1B['Products/Information/AntennaBeamwidth_units'      ] = ['degrees']
L1B['Products/Information/AntennaBeamwidth_description'] = ['Antenna one-way 3dB beamwidth']
L1B['Products/Information/AntennaSize_units'           ] = ['meters']
L1B['Products/Information/AntennaSize_description'     ] = ['none'  ]
L1B['Products/Information/AveragedPulses_units'        ] = ['unitless']
L1B['Products/Information/AveragedPulses_description'   ] = ['Number of pulses averaged per profile, centered around profile time']
L1B['Products/Information/Frequency_units'             ] = ['Hz']
L1B['Products/Information/Frequency_description'       ] = ['Radar transmit frequency']
L1B['Products/Information/GateSpacing_units'           ] = ['meters']
L1B['Products/Information/GateSpacing_description'     ] = ['Range gate spacing']
L1B['Products/Information/PRI_units'                   ] = ['microseconds']
L1B['Products/Information/PRI_description'             ] = ['Pulse repetition interval.']
L1B['Products/Information/Range_units'                 ] = ['meters']
L1B['Products/Information/Range_description'           ] = ['Range in meters from the aircraft']
L1B['Products/Information/Wavelength_description'      ] = ['Radar transmit wavelength']
L1B['Products/Information/Wavelength_units'            ] = ['meters']
L1B['Products/Information/NominalAntennaPointing'      ] = ['Nadir' ]

#######################Populate the Navigation/Data #######################
print ('Adding Navigation')
L1B['Navigation'            ] = {}
L1B['Navigation/Data'       ] = {}
L1B['Navigation/Information'] = {}

##Add navigation data
L1B['Navigation/Data/Drift'        ] = L1A['Navigation/Data/Drift'        ][:].T[0][times2save]##<---
L1B['Navigation/Data/EastVelocity' ] = L1A['Navigation/Data/EastVelocity' ][:].T[0][times2save]##<---
L1B['Navigation/Data/Heading'      ] = L1A['Navigation/Data/Heading'      ][:].T[0][times2save]##<---
L1B['Navigation/Data/Height'       ] = L1A['Navigation/Data/Height'       ][:].T[0][times2save]##<---
L1B['Navigation/Data/Longitude'    ] = L1A['Navigation/Data/Longitude'    ][:].T[0][times2save]##<---##<---NONE IN FILE
L1B['Navigation/Data/NorthVelocity'] = L1A['Navigation/Data/NorthVelocity'][:].T[0][times2save]##<---
L1B['Navigation/Data/Latitude'     ] = L1A['Navigation/Data/Latitude'     ][:].T[0][times2save]##<---##<---NONE IN FILE
L1B['Navigation/Data/Pitch'        ] = L1A['Navigation/Data/Pitch'        ][:].T[0][times2save]##<---
L1B['Navigation/Data/Roll'         ] = L1A['Navigation/Data/Roll'         ][:].T[0][times2save]##<---
L1B['Navigation/Data/Track'        ] = L1A['Navigation/Data/Track'        ][:].T[0][times2save]##<---
L1B['Navigation/Data/UpVelocity'   ] = L1A['Navigation/Data/UpVelocity'   ][:].T[0][times2save]##<---
L1B['Navigation/Data/dxdr'         ] = L1A['Navigation/Data/dx'           ][:].T[0][times2save]##<---
L1B['Navigation/Data/dydr'         ] = L1A['Navigation/Data/dy'           ][:].T[0][times2save]##<---
L1B['Navigation/Data/dzdr'         ] = L1A['Navigation/Data/dz'           ][:].T[0][times2save]##<---

##Calculate speed and nominal distance
spd = np.abs(L1B['Navigation/Data/EastVelocity'][:] + 1j*L1B['Navigation/Data/NorthVelocity'][:])
spd[np.isnan(spd)] = np.nanmedian(spd) #replace missing values with median speed
L1B['Navigation/Data/NominalDistance'] = np.float32(np.diff(L1B['Time/Data/TimeUTC'][0:2])*np.cumsum(spd)) #distance = rate * time

##Add units and information
L1B['Navigation/Information/PitchCorrection'              ] = L1A['Navigation/Information/PitchCorrection'                ]
L1B['Navigation/Information/RollCorrection'               ] = L1A['Navigation/Information/RollCorrection'                 ] 
L1B['Navigation/Information/HeadingCorrection'            ] = L1A['Navigation/Information/HeadingCorrection'              ]
L1B['Navigation/Information/PitchCorrection_units'        ] = ['degrees'                                                  ]
L1B['Navigation/Information/PitchCorrection_description'  ] = ['Pitch offset applied relative to navigation data source'  ]
L1B['Navigation/Information/RollCorrection_units'         ] = ['degrees'                                                  ]
L1B['Navigation/Information/RollCorrection_description'   ] = ['Roll offset applied relative to navigation data source'   ]
L1B['Navigation/Information/HeadingCorrection_units'      ] = ['degrees'                                                  ]
L1B['Navigation/Information/HeadingCorrection_description'] = ['Heading offset applied relative to navigation data source']
L1B['Navigation/Information/NavigationSource'             ] = [config['L1A_navSource']+' Navigation data corrected to nominal radar antenna pointing' ]
L1B['Navigation/Information/Drift_units'                  ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Drift_description'            ] = ['Difference between track and heading'                                                 ]
L1B['Navigation/Information/EastVelocity_units'           ] = ['m/s'                                                                                  ]
L1B['Navigation/Information/EastVelocity_description'     ] = ['Eastward portion of velocity'                                                         ]
L1B['Navigation/Information/Heading_units'                ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Heading_description'          ] = ['Aircraft heading in degrees from north. 90 degrees is eastward pointing'              ]
L1B['Navigation/Information/Height_units'                 ] = ['meters'                                                                               ]
L1B['Navigation/Information/Height_description'           ] = ['Aircraft height above sea level'                                                      ]
L1B['Navigation/Information/NorthVelocity_units'          ] = ['m/s'                                                                                  ]
L1B['Navigation/Information/NorthVelocity_description'    ] = ['Northward portion of velocity'                                                        ]
L1B['Navigation/Information/Latitude_units'               ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Latitude_description'         ] = ['Latitude'                                                                             ]
L1B['Navigation/Information/Longitude_units'              ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Longitude_description'        ] = ['Longitude'                                                                            ]
L1B['Navigation/Information/Pitch_units'                  ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Pitch_description'            ] = ['Pitch'                                                                                ]
L1B['Navigation/Information/Roll_units'                   ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Roll_description'             ] = ['Roll'                                                                                 ]
L1B['Navigation/Information/Track_units'                  ] = ['degrees'                                                                              ]
L1B['Navigation/Information/Track_description'            ] = ['Direction from motion in degrees from north. 90 degrees is eastward motion'           ]
L1B['Navigation/Information/UpVelocity_units'             ] = ['m/s'                                                                                  ]
L1B['Navigation/Information/UpVelocity_description'       ] = ['Upward Velocity'                                                                      ]
L1B['Navigation/Information/dxdr_units'                   ] = ['m/m'                                                                                  ]
L1B['Navigation/Information/dxdr_description'             ] = ['Data cross-track distance from aircraft per radar range. Positive is in the starboard direction.']
L1B['Navigation/Information/dydr_units'                   ] = ['m/m'                                                                                             ]
L1B['Navigation/Information/dydr_description'             ] = ['Data along-track distance from aircraft per radar range. Positive is in the forward direction.'  ]
L1B['Navigation/Information/dzdr_units'                   ] = ['m/m'                                                                                             ]
L1B['Navigation/Information/dzdr_description'             ] = ['Data vertical distance from aircraft per radar range. Positive is in the upward direction.'      ]    
L1B['Navigation/Information/NominalDistance_units'        ] = ['meters'                                                                                          ]
L1B['Navigation/Information/NominalDistance_description'   ] = ['Nominal aircraft travel distance estimated by cumulative summing instantaneous velocities.'      ]

##Define subroutine to calculate horizontal resolution based on beamwidth and range
def res (ranges, config):
    winHann = lambda N: 0.5 * (1-np.cos(2*np.pi*np.arange(1,N+1)/(N+1)))
    import numpy as np
    res = np.zeros(len(ranges), dtype = np.float32)  
    for r in range (len(ranges)): ##for all ranges
        #%x = conv(nominalmotion,antennabeam^2)
        #%res = x>0.25
        x    = np.convolve(np.ones(int(np.floor(200*(config['L1B_LowPRF']+config['L1B_HighPRF'])/2.*config['L1B_averagedPulses']))), winHann(np.floor(2*ranges[r]*np.sin(np.deg2rad(config['L1B_beamWidth']))))**2)
        x    = x/np.max(x) ##normalize to the max
        res[r] = np.sum(x>0.25)
        if r > 1:
            if  res[r] < res[r-1]: res[r] = res[r-1]
    return(res)

##Populate the horiz/vertical resolutions
horiz_res = res(ranges, config)
L1B['Products/Information/ResolutionHorizontal6dB'               ] = horiz_res
L1B['Products/Information/ResolutionHorizontal6dB_units'         ] = ['meters']
L1B['Products/Information/ResolutionHorizontal6dB_description'   ] = ['Approximate horizontal resolution defined as width of spatial weighting after averaging as a function of radar range']
if ((radName == 'CRS') | (radName == 'EXRAD_NADIR')):   
    L1B['Products/Information/ResolutionVertical6dB'             ] = np.float32([config['L1B_verticalResolution']])
    L1B['Products/Information/ResolutionVertical6dB_units'       ] = ['meters']
    L1B['Products/Information/ResolutionVertical6dB_description' ] = ['Approximate vertical resolution defined as width of -6dB points of range weighting function']

##Subroutine to create a single channel mask variable from three individual L1A variables
def mask_sigma(L1A, channelNumber, rangeIndex):
    import numpy as np    
    array = np.zeros(L1A['Products/' + channelNumber + '/Information/Mask1Sigma'].shape, dtype = 'int8') #<--------------------
    array[L1A['Products/' + channelNumber + '/Information/Mask1Sigma'][:]==1] = 1 #<--------------------
    array[L1A['Products/' + channelNumber + '/Information/Mask2Sigma'][:]==1] = 2 #<--------------------
    array[L1A['Products/' + channelNumber + '/Information/Mask3Sigma'][:]==1] = 3 #<--------------------
    array = array[:,rangeIndex] ##cut by rangeIndex
    return(array)
    
##Populate CoPol/CrPol Masks
if radName == 'CRS': ##both CoPol (1)/CrPol (2)
    L1B['Products/Information/MaskCoPol'            ] = mask_sigma(L1A, 'CoPolChirp', rangeIndex) #<--------------------
    L1B['Products/Information/MaskCoPol_description'] = ['Co-polarization SNR Mask: (Mask >= #) corresponds with (SNR > # sigma noise)']
    L1B['Products/Information/MaskCrPol'            ] = mask_sigma(L1A, 'CrossPolChirp', rangeIndex) #<--------------------
    L1B['Products/Information/MaskCrPol_description'] = ['Cr-polarization SNR Mask: (Mask >= #) corresponds with (SNR > # sigma noise)']
if radName == 'EXRAD_NADIR': ##CoPol only, (0)
    L1B['Products/Information/MaskCoPol'            ] = mask_sigma(L1A, 'Nadir', rangeIndex) ##<-------------------------------
    L1B['Products/Information/MaskCoPol_description'] = ['Co-polarization SNR Mask: (Mask >= #) corresponds with (SNR > # sigma noise)']

#Function to calculate dBZe from Zuncal
def dbze  (zuncal, extcalCoeff, wavelength, k_factor): 
    import numpy as np
    db = (10*np.log10(zuncal) + extcalCoeff + 10*np.log10(wavelength**4*10**18/np.pi**5/k_factor)).astype(np.float32)
    return (db)

#Subroutine to process sigma0 from NRCSuncal
def sigma (NRCSuncal, extCalCoeff):
    import numpy as np
    NRCSuncal[NRCSuncal <= 0] = np.nan               #Filter negative numbers
    NRCSuncal = 10*np.log10(NRCSuncal) + extCalCoeff #convert to dB, add external cal coefficient
    return (NRCSuncal)

#######################ADD Products/Data#######################
##For CRS/EXRAD/HIWRAP(uncombined), add Product/Data, using above subroutines when needed
##For EXRAD and HIWRAP Velocities, remove the L1A AircraftMotion and add the smoothed L1B AircraftMotion
#<---------------------------------------------------------------------------------------------------------------------------changing velocity_uncorrected to "velocity" for EX ALOFT, no HRRR correction
#<---------------------------------------------------------------------------------------------------------------------------changing velocity_uncorrected to "velocity" for EX ALOFT, no HRRR correction
#<---------------------------------------------------------------------------------------------------------------------------changing velocity_uncorrected to "velocity" for EX ALOFT, no HRRR correction
#<---------------------------------------------------------------------------------------------------------------------------changing velocity_uncorrected to "velocity" for EX ALOFT, no HRRR correction
#<---------------------------------------------------------------------------------------------------------------------------changing velocity_uncorrected to "velocity" for EX ALOFT, no HRRR correction
#<---------------------------------------------------------------------------------------------------------------------------changing velocity_uncorrected to "velocity" for EX ALOFT, no HRRR correction
print ('Adding Product Data')
if radName == 'CRS': #<----------------- CHANGED CHANNEL NUMBERS
    L1B['Products/Data/dBZe'                 ] = (dbze(L1A['Products/CoPolChirp/Data/Zuncal'        ][:,rangeIndex], extCalCoeff[1], L1B['Products/Information/Wavelength'], config['L1B_Kfactor'])).astype('float32')
    L1B['Products/Data/Velocity_uncorrected' ] = (L1A['Products/CoPolChirp/Data/V'                  ][:,rangeIndex]).astype('float32')
    #L1B['Products/Data/SpectrumWidth'        ] = L1A['Products/CoPolChirp/Data/SpecWid'            ][:,rangeIndex]#<--------------------
    L1B['Products/Information/SNR'           ] = (L1A['Products/CoPolChirp/Information/SNR'         ][:,rangeIndex]).astype('float32')
    L1B['Products/Data/sigma0'               ] = sigma(L1A['Products/CoPolChirp/Data/NRCSuncal'    ][:], extCalCoeff[1])
    L1B['Products/Information/noiseFloor'] = L1A['Products/CoPolChirp/Information/EstimatedNoise'][:].T[0]
    #L1B['Products/Data/LDR'                  ] = (10*np.log10(L1A['Products/CrossPolChirp/Data/Zuncal'][:,rangeIndex].T/L1A['Products/CoPolChirp/Information/CalEnergy'][:].T)).T- 14.75- 10*np.log10(L1A['Products/CoPolChirp/Data/Zuncal'][:,rangeIndex])#<--------------------  
    ##New LDR calculation from Matt M on Aug 16 2023
    #L1B['Products/Data/LDR'                  ] = 10*np.log10(((L1A['Products/CrossPolChirp/Data/Zuncal'][:,rangeIndex].T/L1A['Products/CrossPolChirp/Information/EstimatedNoise'][:,0]).T)/((L1A['Products/CoPolChirp/Data/Zuncal'][:,rangeIndex].T/L1A['Products/CoPolChirp/Information/EstimatedNoise'][:,0]).T))#<--------------------    
    ##New LDR calculation from Matt M on Aug 16 2023    
    L1B['Products/Data/LDR'                  ] = (10*np.log10(L1A['Products/CrossPolChirp/Information/SNR'][:,rangeIndex]/L1A['Products/CoPolChirp/Information/SNR'][:,rangeIndex])).astype('float32')#<--------------------
    
if radName == 'EXRAD_NADIR':#<----------------- CHANGED CHANNEL NUMBERS
    L1B['Products/Data/dBZe'                ] = (dbze(L1A['Products/Nadir/Data/Zuncal'    ][:,rangeIndex], extCalCoeff[0], L1B['Products/Information/Wavelength'], config['L1B_Kfactor'])).astype('float32')
    L1B['Products/Data/Velocity'] = (L1A['Products/Nadir/Data/V'             ][:,rangeIndex]).astype('float32')
    #L1B['Products/Data/SpectrumWidth'       ] = np.real(L1A['Products/Nadir/Data/SpecWid'        ][:,rangeIndex])#<--------------------
    L1B['Products/Information/SNR'          ] = (L1A['Products/Nadir/Information/SNR'     ][:,rangeIndex]).astype('float32')
    L1B['Products/Data/sigma0'              ] = sigma(L1A['Products/Nadir/Data/NRCSuncal'][:], extCalCoeff[0])
    L1B['Products/Information/noiseFloor'] = L1A['Products/Nadir/Information/EstimatedNoise'][:].T[0]
        
    if flightDate == '20220213':
        L1B['Products/Data/Velocity_uncorrected'] = (L1B['Products/Data/Velocity_uncorrected']*np.nan).astype('float32')
        L1B['Products/Data/SpectrumWidth']        = (L1B['Products/Data/SpectrumWidth']*np.nan).astype('float32')
        print ('Special EXRAD case, deleting V and SpecWid arrays')

if radName == 'CRS': ##add LDR for CRS. 2sigma filter
    L1B['Products/Data/LDR'][L1B['Products/Information/MaskCrPol'] <2 ] = np.nan #cross
    L1B['Products/Data/LDR'][L1B['Products/Information/MaskCoPol'] <2 ] = np.nan #co
    L1B['Products/Data/LDR'][np.imag(L1B['Products/Data/LDR'])     !=0] = np.nan

##Filter by 1 sigma for dBZe, Velocity, SpectrumWidth
if ((radName == 'CRS') | (radName == 'EXRAD_NADIR')): 
    L1B['Products/Data/dBZe'                ][L1B['Products/Information/MaskCoPol'] <1 ] = np.nan
    L1B['Products/Data/Velocity'][L1B['Products/Information/MaskCoPol'] <1 ] = np.nan
    #L1B['Products/Data/SpectrumWidth'][L1B['Products/Information/MaskCoPol'] <1 ] = np.nan#<--------------------

##Add units and descriptions
if ((radName == 'CRS') | (radName == 'EXRAD_NADIR')):   
    if radName == 'CRS': ##Only CRS has LDR 
        L1B['Products/Information/dBZe_description'     ] = ['Equivalent reflectivity factor in dB with 1-sigma noise threshold applied. |K|^2 = '+str(config['L1B_Kfactor'])+' rather than 0.93 for consistency with CloudSat.']
        L1B['Products/Information/LDR_units'            ] = ['dB']
        L1B['Products/Information/LDR_description'      ] = ['Linear Depolarization Ratio with 2-sigma co- and cross-polarization noise thresholding applied']
    if radName == 'EXRAD_NADIR': 
        L1B['Products/Information/dBZe_description'            ] = ['Equivalent reflectivity factor in dB with 1-sigma noise threshold applied. |K|^2 = '+str(config['L1B_Kfactor'])]
    L1B['Products/Information/dBZe_units'                      ] = ['10*log10(mm^6/m^3)']
    L1B['Products/Information/Velocityunits'      ] = ['m/s']
    L1B['Products/Information/Velocity_description'] = ['Doppler velocity with aircraft motion correction, NUBF correction, and 1-sigma noise threshold applied. Positive velocity is upward'    ]
    L1B['Products/Information/SpectrumWidth_units'             ] = ['m/s']
    L1B['Products/Information/SpectrumWidth_description'       ] = ['Doppler velocity spectrum width estimate including aircraft motion. 1-sigma noise threshold applied.'    ]
    L1B['Products/Information/SNR_units'                       ] = ['Watt / Watt']
    L1B['Products/Information/SNR_description'                 ] = ['Estimated Signal to Noise Ratio'    ]
    L1B['Products/Information/sigma0_units'                    ] = ['10*log10(m^2/m^2)']
    L1B['Products/Information/sigma0_description'              ] = ['Ocean Normalized Radar Cross Section, only valid over ocean'    ]
    L1B['Products/Information/noiseFloor_units'                ] = ['Relative power']
    L1B['Products/Information/noiseFloor_description'          ] = ['Uncalibrated noise floor estimate']

##Subroutine to perform NUBF correction
def nubf (dist, spd, ranges, beam, tempdbze, gradkern):
    from scipy import ndimage
    import numpy  as np
    
    ##Constant defined in NUBF paper by Loftus and McLinden
    alpha = lambda spd, rng, beamwidth: spd*0.0207*rng*np.deg2rad(beamwidth)**2
    
    ##Make a 2d array (tile) of speed so that it is a function of range.
    ##Use this for the alpha calculation.
    spd2d = np.tile(spd, (len(ranges), 1)).T
    
    ##Calculate the constant based on alpha equation and the spd2d
    constant = alpha(spd2d, ranges/1000., beam)
    
    ##Calculate the dBZ gradient in km
    ##Start with DY (distance)
    DY = np.convolve(dist, -gradkern, 'same')/1000 #m -> km
    
    ##Plane flies at 200m/s for 2 secs. Expect around 0.4km distance. Filter out outliers.
    DY[DY<.3] = np.nan
    DY[DY>.6] = np.nan
    
    ##Now calculate dBZe/DY. ndimage allows convolution on a 2d array.   
    grad = (ndimage.filters.convolve1d(tempdbze.T, -gradkern, axis=1, mode='constant')/DY).T

    ##The offset is the product of the alpha constant * gradient
    offset = grad * constant
    return (offset)
    
##Find radar-specific NUBF variables, then do the calculation
if radName == 'EXRAD_NADIR' : 
    gradkern = np.array([-1, 0, 0, 0, 1]) #8 sample, 2 second kernel
    beam     = L1B['Products/Information/AntennaBeamwidth']
    tempdbze = np.copy(L1B['Products/Data/dBZe'])
    tempdbze[L1B['Products/Information/MaskCoPol'] < 2] = np.nan #do some filtering. #####################################################WHY
    offset = nubf(L1B['Navigation/Data/NominalDistance'], spd, ranges, beam, tempdbze, gradkern)
    
    L1B['Products/Information/Velocity_nubf_offset'                        ] = offset.astype(np.float32) ##RevB was .1438*2, corrected Feb 26 2021
    L1B['Products/Information/Velocity_nubf_offset_units'                  ] = ['m/s']
    L1B['Products/Information/Velocity_nubf_offset_description'            ] = ['The nonuniform beam filling (NUBF) offset used to correct Doppler velocity [Vel_corr = Vel_uncorr - nubf_offset - horizwind_offset]']

    L1B['Products/Data/Velocity'] = L1B['Products/Data/Velocity']  - L1B['Products/Information/Velocity_nubf_offset' ]  #<----------------------------------new for RevB, add NUBF correction to Velocity
    
##Save the L1B dict to HDF file
print ('Writing dict')
dicttoh5(L1B, L1B_filepath)
h = h5py.File(L1B_filepath, 'a')

##Calculate the HRRR Winds along the flight path
print ('Processing winds')
L1B_uwind, L1B_vwind = L1B_sub_add_horiz_wind_hrrr_fullflight.L1B_sub_add_horiz_wind ('HRRR', radName, L1B_filepath, m_filepath)

##Calculate the along/crosstrack winds and return the winds.
##The L1B[Vel] variable is for plotting only.
print ('Reprocessing doppler')
if 'HIWRAP' in radName:
    cross, along, vel_corr_cross, vel_corr_along = \
        L1B_sub_along_cross_winds2.L1B_sub_along_cross_winds (radName, flightDate, L1B, L1B_uwind, L1B_vwind, L1B['Products/Combined/Data/Velocity_uncorrected'][:], cmap_n0q, deg_offset)
else:
    cross, along, vel_corr_cross, vel_corr_along = \
        L1B_sub_along_cross_winds2.L1B_sub_along_cross_winds (radName, flightDate, L1B, L1B_uwind, L1B_vwind, L1B['Products/Data/Velocity_uncorrected'][:], cmap_n0q, deg_offset)
##Reopen the newly created file and add these new variables

print ('Adding cross/along wind information')
h.create_dataset('/Products/Information/HRRR_CrossWind',             data=cross)
h.create_dataset('/Products/Information/HRRR_AlongWind',             data=along)
h.create_dataset('/Products/Information/HRRR_CrossWind_units',       data=np.string_(['m/s']))
h.create_dataset('/Products/Information/HRRR_AlongWind_units',       data=np.string_(['m/s']))
h.create_dataset('/Products/Information/HRRR_CrossWind_description', data=np.string_(['Cross wind from HRRR reanalysis, interpolated to this grid']))
h.create_dataset('/Products/Information/HRRR_AlongWind_description', data=np.string_(['Along wind from HRRR reanalysis, interpolated to this grid']))


#<-------------------------------------------------------------------------------------------------------------------------------------------------------------------- this doesnt run for ALOFT2023
if (radName == 'CRS'):###<-----------Manual DV Correction from Matt, Oct 6 2023
    h.create_dataset('/Products/Data/Velocity_corrected',                                      data=(L1B['Products/Data/Velocity_uncorrected'] - vel_corr_cross -vel_corr_along).astype('float32'))
    #h.create_dataset('/Products/Data/Velocity_corrected',                                      data=L1B['Products/Data/Velocity_uncorrected'] - (np.tile(L1B['Navigation/Data/dydr'], (len(ranges), 1)).T * along) + (np.tile(L1B['Navigation/Data/dxdr'], (len(ranges), 1)).T * cross))
    h.create_dataset('/Products/Information/Velocity_corrected_units',                         data=np.string_(['m/s']))
    h.create_dataset('/Products/Information/Velocity_corrected_description',                   data=np.string_(['Doppler velocity corrected to account for intrusion of horizontal reanalysis winds']))
    h.create_dataset('/Products/Information/Velocity_horizwind_offset',                        data=(vel_corr_cross +vel_corr_along).astype('float32'))
    h.create_dataset('/Products/Information/Velocity_horizwind_offset_units',                  data=np.string_(['m/s']))
    h.create_dataset('/Products/Information/Velocity_horizwind_offset_description',            data=np.string_(['The horizontal wind offset used to correct Doppler velocity [Vel_corr = Vel_uncorr - horizwind_offset]']))
if (radName == 'EXRAD_NADIR'):
    h.create_dataset('/Products/Data/Velocity_corrected',                                      data=(L1B['Products/Data/Velocity_uncorrected'] -L1B['Products/Information/Velocity_nubf_offset']- vel_corr_cross -vel_corr_along).astype('float32')) #NUBF and winds
    h.create_dataset('/Products/Information/Velocity_corrected_units',                         data=np.string_(['m/s']))
    h.create_dataset('/Products/Information/Velocity_corrected_description',                   data=np.string_(['Doppler velocity corrected to account for intrusion of horizontal reanalysis winds and nonuniform beam filling (NUBF)']))
    h.create_dataset('/Products/Information/Velocity_horizwind_offset',                        data=(vel_corr_cross +vel_corr_along).astype('float32'))
    h.create_dataset('/Products/Information/Velocity_horizwind_offset_units',                  data=np.string_(['m/s']))
    h.create_dataset('/Products/Information/Velocity_horizwind_offset_description',            data=np.string_(['The horizontal wind offset used to correct Doppler velocity [Vel_corr = Vel_uncorr - nubf_offset - horizwind_offset]']))

if ((makePlots == 'yes') | (makePlots == 'y')):
    print ('Plotting')
    if radName == 'CRS':
        L1B_sub_plotter.L1B_sub_plotter (radName,       h, flightLines, L1B_file, L1B_plotpath, cmap_n0q, flightDate, versionLetter)
    if radName == 'EXRAD_NADIR':
        L1B_sub_plotter.L1B_sub_plotter ('EXRAD_NADIR', h, flightLines, L1B_file, L1B_plotpath, cmap_n0q, flightDate, versionLetter)
h.close()