""" The SNR Calculations implemented is very similar to SIGNAL code
given at following location
https://www.ing.iac.es//astronomy/instruments/signal/signal_code.html """

from __future__ import print_function
import numpy as np
import sys
from math import exp,log,atan2,asin,cos,sin
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.integrate import quad,dblquad

from WEAVEify_routines.step_02_common_routines.python_2_3_compatibility_routines \
import genfromtxt_compatible

def opr3_mock_noise_routines_main(inp_dict_org):
    inp_dict = step_1_define_parameters(inp_dict_org)
    inp_dict = step_2_compute_telescope_area(inp_dict)
    inp_dict = step_3_compute_telescope_throughput_OPR3(inp_dict)
    inp_dict = step_4_compute_telescope_effective_area(inp_dict)
    inp_dict = step_5_compute_dispersion_in_Ang_per_pixels(inp_dict)
    inp_dict = step_6_fraction_of_light_passed_by_fiber_vignette(inp_dict)
    inp_dict = step_7_get_waveband_normalized_simulated_flux(inp_dict)
    inp_dict = step_8_compute_snr(inp_dict)
    inp_dict = step_9_compute_flux_normalization_from_magnitude(inp_dict)
    inp_dict = step_10_add_noise_to_spectra(inp_dict)
    return inp_dict







def step_1_define_parameters(inp_dict_org):
    inp_dict = inp_dict_org.copy()
    fFP      = float(inp_dict["fFP"])
    fcam     = float(inp_dict["fcam"])
    fcol     = float(inp_dict["fcol"])

    dia_fibre_arcsec    = float(inp_dict["dia_fibre_arcsec"])
    ccd_pix_size_micron = float(inp_dict["ccd_pix_size_micron"])

    # Resolving power of the spectrograph
    resolution_mode = inp_dict["resolution_mode"]
    R_val_dict      = {"LR":5750,"HR":21000}
    R_val           = R_val_dict[resolution_mode]


    # Diameter of Fibre in microns
    dia_fibre_micron = dia_fibre_arcsec * fFP * 4.2 * 1e6 / 206265

    # Diameter of Fibre in units of CCD pixels
    dia_fibre_pix       = dia_fibre_micron * (fcam/fcol) / ccd_pix_size_micron

    wave_band_dict      = initilize_waveband()


    out_dict            = {"wave_band_dict":wave_band_dict,\
                           "dia_fibre_micron":dia_fibre_micron,\
                           "dia_fibre_pix":dia_fibre_pix,\
                           "R_val":R_val}

    inp_dict.update(out_dict)
    return inp_dict


##########################################################################
##########################################################################
##########################################################################

def step_2_compute_telescope_area(inp_dict):
    dia_telescope  = float(inp_dict["dia_telescope"])
    r_telescope    = dia_telescope * 0.5
    area_telescope = np.pi * (r_telescope**2.0)
    inp_dict.update({"area_telescope":area_telescope})
    return inp_dict

##########################################################################
##########################################################################
##########################################################################

def step_3_compute_telescope_throughput_OPR3(inp_dict):
    resolution_mode = inp_dict["resolution_mode"]

    if resolution_mode == "LR":
        inp_dict = step_3a_compute_telescope_throughput_LR(inp_dict)

    elif resolution_mode == "HR":
        inp_dict = step_3b_compute_telescope_throughput_HR(inp_dict)

    return inp_dict


def step_3a_compute_telescope_throughput_LR(inp_dict):
    """ The regions of throughput are overlapped for LR mode """
    #lambda_arr_obs  = inp_dict["lambda_arr_obs"]
    inp_dict_Flux       = inp_dict["inp_dict_Flux"]
    flux_arr_cont_added = inp_dict_Flux["Flux_arr_cont_added"]
    lambda_arr_obs      = flux_arr_cont_added[:,0]

    inp_path        = inp_dict["inp_path_throughput"]
    inp_file        = inp_dict["inp_file_throughput"]
    throughput_type = inp_dict["throughput_type"]
    resolution_mode = inp_dict["resolution_mode"]
    data_arr_all    = genfromtxt_compatible(inp_path,inp_file,object)

    eff_tel_arr_all    = data_arr_all[:,:]
    lambda_arr_eff_all = data_arr_all[:,3].astype(np.float)*10.0
    tot_eff_arr_all    = np.ones(lambda_arr_eff_all.shape[0])
    for col_indx in range(4,10):
        eff_arr_all_indx = eff_tel_arr_all[:,col_indx].astype(np.float)
        tot_eff_arr_all  = tot_eff_arr_all * (eff_arr_all_indx / 100.0)


    bool_arr_1        = np.char.upper(data_arr_all[:,0].astype(str)) == resolution_mode
    bool_arr_2        = np.char.lower(data_arr_all[:,2].astype(str)) == throughput_type


    eff_tel_arr_intrp = np.zeros(lambda_arr_obs.shape[0])


    bool_arr_arm_blue = np.char.lower(data_arr_all[:,1].astype(str)) == "blue"
    bool_arr_arm_red  = np.char.lower(data_arr_all[:,1].astype(str)) == "red"

    bool_arr_blue     = bool_arr_1 & bool_arr_2 & bool_arr_arm_blue
    bool_arr_red      = bool_arr_1 & bool_arr_2 & bool_arr_arm_red

    lambda_arr_eff_blue = lambda_arr_eff_all[bool_arr_blue]
    lambda_arr_eff_red  = lambda_arr_eff_all[bool_arr_red]

    tot_eff_arr_blue  = tot_eff_arr_all[bool_arr_blue]
    tot_eff_arr_red   = tot_eff_arr_all[bool_arr_red]


    intrp_funct_blue = interp1d(lambda_arr_eff_blue,tot_eff_arr_blue,\
                       bounds_error=False,fill_value=(tot_eff_arr_blue[0],0.0))
    intrp_funct_red  = interp1d(lambda_arr_eff_red,tot_eff_arr_red,\
                       bounds_error=False,fill_value=(0.0,tot_eff_arr_red[-1]))

    eff_tel_arr_intrp_blue = intrp_funct_blue(lambda_arr_obs)
    eff_tel_arr_intrp_red  = intrp_funct_red(lambda_arr_obs)

    # Addition because there is overlap
    eff_tel_arr_intrp      = eff_tel_arr_intrp_blue + eff_tel_arr_intrp_red

    inp_dict.update({"eff_tel_arr":eff_tel_arr_intrp})
    return inp_dict

def step_3b_compute_telescope_throughput_HR(inp_dict):
    """ The regions of throughput do not have overlap for HR mode """
    #lambda_arr_obs  = inp_dict["lambda_arr_obs"]
    inp_dict_Flux       = inp_dict["inp_dict_Flux"]
    flux_arr_cont_added = inp_dict_Flux["Flux_arr_cont_added"]
    lambda_arr_obs      = flux_arr_cont_added[:,0]
    inp_path        = inp_dict["inp_path_throughput"]
    inp_file        = inp_dict["inp_file_throughput"]
    throughput_type = inp_dict["throughput_type"]
    resolution_mode = inp_dict["resolution_mode"]
    data_arr_all    = genfromtxt_compatible(inp_path,inp_file,object)

    eff_tel_arr_all    = data_arr_all[:,:]
    lambda_arr_eff_all = data_arr_all[:,3].astype(np.float)*10.0
    tot_eff_arr_all    = np.ones(lambda_arr_eff_all.shape[0])
    for col_indx in range(4,10):
        eff_arr_all_indx = eff_tel_arr_all[:,col_indx].astype(np.float)
        tot_eff_arr_all  = tot_eff_arr_all * (eff_arr_all_indx / 100.0)


    bool_arr_1        = np.char.upper(data_arr_all[:,0].astype(str)) == resolution_mode
    bool_arr_2        = np.char.lower(data_arr_all[:,2].astype(str)) == throughput_type


    eff_tel_arr_intrp = np.zeros(lambda_arr_obs.shape[0])

    bool_arr_arm_blue  = np.char.lower(data_arr_all[:,1].astype(str)) == "blue"
    bool_arr_arm_green = np.char.lower(data_arr_all[:,1].astype(str)) == "green"
    bool_arr_arm_red   = np.char.lower(data_arr_all[:,1].astype(str)) == "red"

    bool_arr_blue     = bool_arr_1 & bool_arr_2 & bool_arr_arm_blue
    bool_arr_green    = bool_arr_1 & bool_arr_2 & bool_arr_arm_green
    bool_arr_red      = bool_arr_1 & bool_arr_2 & bool_arr_arm_red

    lambda_arr_eff_blue  = lambda_arr_eff_all[bool_arr_blue]
    lambda_arr_eff_green = lambda_arr_eff_all[bool_arr_green]
    lambda_arr_eff_red   = lambda_arr_eff_all[bool_arr_red]
    lambda_arr_eff_arm   = np.concatenate([lambda_arr_eff_blue,lambda_arr_eff_green,lambda_arr_eff_red])

    tot_eff_arr_blue   = tot_eff_arr_all[bool_arr_blue]
    tot_eff_arr_green  = tot_eff_arr_all[bool_arr_green]
    tot_eff_arr_red    = tot_eff_arr_all[bool_arr_red]

    # Since there is no overlap we can concatenate the fields
    tot_eff_arr_arm    = np.concatenate([tot_eff_arr_blue,tot_eff_arr_green,tot_eff_arr_red])

    intrp_funct_arm    = interp1d(lambda_arr_eff_arm,tot_eff_arr_arm,\
                         bounds_error=False,fill_value=(tot_eff_arr_arm[0],\
                         tot_eff_arr_arm[-1]))

    eff_tel_arr_intrp  = intrp_funct_arm(lambda_arr_obs)

    inp_dict.update({"eff_tel_arr":eff_tel_arr_intrp})
    return inp_dict

##########################################################################
##########################################################################
##########################################################################

def step_4_compute_telescope_effective_area(inp_dict):
    area_telescope         = inp_dict["area_telescope"]
    eff_tel_arr            = inp_dict["eff_tel_arr"]
    eff_area_telescope_arr = area_telescope * eff_tel_arr
    inp_dict.update({"eff_area_telescope_arr":eff_area_telescope_arr})
    return inp_dict

##########################################################################
##########################################################################
##########################################################################


def step_5_compute_dispersion_in_Ang_per_pixels(inp_dict):
    #def dispersion(self,cwave,fcam,fcol,beamsize):
    fcam             = float(inp_dict["fcam"])
    fcol             = float(inp_dict["fcol"])
    beamsize         = float(inp_dict["beamsize"])
    R_val            = float(inp_dict["R_val"])
    dia_fibre_micron = float(inp_dict["dia_fibre_micron"])
    ccd_pix_size_micron = float(inp_dict["ccd_pix_size_micron"])
    cwave_arr           = step_5a_compute_central_wavelength_array(inp_dict)

    # compute dispersion in Ang/pixel at central wavelength
    # Focal length of camera
    L_focal          = beamsize * fcam


    # Grating angle calculation in radians
    gamma            = atan2(1.e-3 * dia_fibre_micron * R_val, 2 * beamsize * fcol)

    # l/mm of grating
    l_per_mm         = 2 * sin(gamma) / (1.e-7*cwave_arr)
    l_per_micron     = 1.e-3 * l_per_mm


    # Angle subtended by 1 pixel
    pix_angle        = ccd_pix_size_micron * 1.e-3 / L_focal

    # diffracted angle of one pixel (approximation good to very high precision!)
    diff_angle       = np.arcsin(np.sin(gamma) - 1.e-4 * cwave_arr * l_per_micron)

    # dispersion in ang/pixel
    disp_A_per_pixel_arr = 1e4 * pix_angle * np.cos(diff_angle) / l_per_micron

    inp_dict.update({"disp_A_per_pixel_arr":disp_A_per_pixel_arr})

    return inp_dict

def step_5a_compute_central_wavelength_array(inp_dict):
    lambda_param_dict = inp_dict["lambda_arm_param_dict"]
    inp_dict_Flux       = inp_dict["inp_dict_Flux"]
    flux_arr_cont_added = inp_dict_Flux["Flux_arr_cont_added"]
    lambda_arr          = flux_arr_cont_added[:,0]
    resolution_mode     = inp_dict["resolution_mode"]

    if resolution_mode == "LR":
        arm_type_list = ["blue","red"]

    elif resolution_mode == "HR":
        arm_type_list = ["blue","green","red"]

    cwave_arr      = np.zeros(lambda_arr.shape[0])
    lambda_min     = 0

    for arm_indx in range(len(arm_type_list)):
        cwave_prev = lambda_param_dict[(resolution_mode,arm_type_list[arm_indx])]["cwave"]
        if arm_indx == (len(arm_type_list) - 1):
            lambda_max = 1e5
        else:
            lambda_max = lambda_param_dict[(resolution_mode,arm_type_list[arm_indx+1])]["CRVAL1"]
        bool_arr = (lambda_min <= lambda_arr) & (lambda_arr < lambda_max)
        cwave_arr[bool_arr]   = cwave_prev
        lambda_min = lambda_max
    return cwave_arr



##########################################################################

def step_6_fraction_of_light_passed_by_fiber_vignette(inp_dict):
    source_type    = inp_dict["source_type"]

    if source_type == "extended":
        dia_fibre_arcsec   = float(inp_dict["dia_fibre_arcsec"])
        fibre_area_arcsec2 = np.pi * dia_fibre_arcsec**2.0 / 4.0
        f_vign             = fibre_area_arcsec2


    elif source_type == "point":

        fwhm_seeing       = float(inp_dict["fwhm_seeing"])
        PRI_fwhm          = float(inp_dict["PRI_fwhm"])
        r_offset          = float(inp_dict["r_offset"])
        dia_fibre_arcsec  = float(inp_dict["dia_fibre_arcsec"])
        seeing_profile    = inp_dict["seeing_profile"]
        fwhm_seeing_total = (fwhm_seeing**2.0 + PRI_fwhm**2.0)**0.5
        r_fibre_arcsec    = 0.5 * dia_fibre_arcsec
        sigma_seeing      = fwhm_seeing_total / 2.355

        if seeing_profile.lower() =='gaussian':
            f_vign        = dblquad(gaussian_2d, -r_fibre_arcsec, r_fibre_arcsec, \
                            lambda x: -(r_fibre_arcsec**2.0 - x**2.0)**0.5, \
                            lambda x: +(r_fibre_arcsec**2.0 - x**2.0)**0.5, \
                            args=(sigma,r_offset))[0]

        elif seeing_profile.lower() =='moffat':
            beta_moffat   = inp_dict["beta_moffat"]
            alpha_moffat  = fwhm_seeing_total / (2.0 * (2.0**(1./beta_moffat) - 1)**0.5)
            f_vign        = dblquad(moffat_2d, -r_fibre_arcsec, r_fibre_arcsec, \
                            lambda x: -(r_fibre_arcsec**2.0 - x**2.0)**0.5, \
                            lambda x: +(r_fibre_arcsec**2.0 - x**2.0)**0.5, \
                            args=(alpha_moffat,beta_moffat,r_offset))[0]

    inp_dict.update({"f_vign":f_vign})

    return inp_dict

##########################################################################

def gaussian_2d(y,x,sigma,r_offset):
    r     = ((x-r_offset)**2 + y**2.0)**0.5
    gauss = exp(-r**2.0 / (2 * sigma**2.0)) / (2 * np.pi * sigma**2.0)
    return gauss

##########################################################################

def moffat_2d(y,x,alpha_moffat,beta_moffat,r_offset):
    r = ((x-r_offset)**2 + y**2.0)**0.5

    # See Moffat function in Mathematica
    if (alpha_moffat >= 0.0) & (beta_moffat > 1.0):
        norm = (2 * np.pi * alpha_moffat**2) / (2*(beta_moffat-1))
    else:
        raise ValueError('alpha_moffat and/or beta_moffat out of bounds')

    moffat = (1.0 + (r/alpha_moffat)**2.0)**(-beta_moffat) / norm
    return moffat

##########################################################################
##########################################################################
##########################################################################

def step_7_get_waveband_normalized_simulated_flux(inp_dict):
    inp_dict_Flux    = inp_dict["inp_dict_Flux"]
    flux_arr_cont    = inp_dict_Flux["Flux_arr_cont_added"]
    inp_path_filter  = inp_dict["inp_path_trans_filter"]
    inp_file_filter  = inp_dict["inp_file_trans_filter"]
    wave_band_source = inp_dict["wave_band_source"]
    filter_arr_all   = np.loadtxt(inp_path_filter + inp_file_filter)

    band_indx_dict   = {"U":1,"B":2,"V":3,"R":4,"I":5}
    indx             = band_indx_dict[wave_band_source.upper()]

    intrp_funct      = interp1d(filter_arr_all[:,0],filter_arr_all[:,indx],bounds_error=False,fill_value=(0.0,0.0))
    filter_arr       = intrp_funct(flux_arr_cont[:,0])
    filter_arr_norm  = filter_arr / np.sum(filter_arr)
    flux_norm_val    = np.sum(filter_arr_norm * flux_arr_cont[:,1])

    flux_sim_arr_norm      = flux_arr_cont.copy()
    flux_sim_arr_norm[:,1] = flux_sim_arr_norm[:,1] / flux_norm_val
    flux_sim_arr_norm[:,2] = flux_sim_arr_norm[:,2] / flux_norm_val
    inp_dict.update({"flux_sim_arr_norm":flux_sim_arr_norm})

    return inp_dict




##########################################################################
##########################################################################
##########################################################################

def step_8_compute_snr(inp_dict):
    # number of pixels along the slit subtended by the fiber
    npix_spatial       = inp_dict["dia_fibre_pix"]
    # number of pixels in 1 Angstrom
    npix_spectral_arr  = 1.0 / inp_dict["disp_A_per_pixel_arr"]
    pixArea_arr        = npix_spatial * npix_spectral_arr
    dia_fibre_pix      = inp_dict["dia_fibre_pix"]
    t_exp              = inp_dict["t_exp"]
    read_noise         = inp_dict["read_noise"]
    dark_current       = inp_dict["dark_current"]
    flux_sim_arr_norm  = inp_dict["flux_sim_arr_norm"]

    s_source_arr       = step_7a_compute_N_photon_for_given_magnitude(inp_dict,flag_source_sky="source")
    s_sky_arr          = step_7a_compute_N_photon_for_given_magnitude(inp_dict,flag_source_sky="sky")
    #s_source_arr       = s_source_arr * flux_sim_arr_norm[:,2]
    s_source_arr       = s_source_arr * flux_sim_arr_norm[:,1]


    # SNR per pixel
    snr_per_pixel_arr  = s_source_arr / (s_source_arr + s_sky_arr + npix_spatial * \
                         (read_noise**2.0 + dark_current * t_exp / 3600.0))**0.5

    # SNR per angstroms
    snr_per_A_arr      = npix_spectral_arr * s_source_arr / (npix_spectral_arr * s_source_arr + \
                         npix_spectral_arr * s_sky_arr + pixArea_arr * (read_noise**2.0 + \
                         dark_current * t_exp / 3600.0))**0.5

    # SNR per resolution element
    snr_per_res_arr    = dia_fibre_pix * s_source_arr / (dia_fibre_pix * s_source_arr + \
                         dia_fibre_pix * s_sky_arr + dia_fibre_pix**2.0 * \
                         (read_noise**2.0 + dark_current * t_exp / 3600.0))**0.5


    out_dict = {"snr_per_pixel_arr":snr_per_pixel_arr,\
                "snr_per_A_arr":snr_per_A_arr,\
                "snr_per_res_arr":snr_per_res_arr}

    inp_dict.update(out_dict)

    return inp_dict

##########################################################################

def step_7a_compute_N_photon_for_given_magnitude(inp_dict,flag_source_sky="source"):
    eff_area_arr         = inp_dict["eff_area_telescope_arr"]
    disp_A_per_pixel_arr = inp_dict["disp_A_per_pixel_arr"]
    wave_band_dict       = inp_dict["wave_band_dict"]
    wave_band_source     = inp_dict["wave_band_source"]
    t_exp                = inp_dict["t_exp"]
    dia_fibre_arcsec     = inp_dict["dia_fibre_arcsec"]
    fibre_area_arcsec2   = np.pi * dia_fibre_arcsec**2.0 / 4.0


    if flag_source_sky.lower() == "source":
        f_vign         = inp_dict["f_vign"]
        m_source_sky   = float(inp_dict["m_source"])
        extinct        = wave_band_dict[wave_band_source]["extinct"]
        airmass        = inp_dict["airmass"]
        extinct_factor = 10**(-extinct * airmass / 2.5)
        f_source_sky   = extinct_factor * f_vign

    elif flag_source_sky.lower() == "sky":
        m_source_sky   = float(inp_dict["m_sky"])
        f_source_sky   = fibre_area_arcsec2


    N_photon_mag_0   = sub_step_7b_compute_N_photon_for_0_magnitude(inp_dict,flag_source_sky)

    term_area        = eff_area_arr * t_exp * f_source_sky * disp_A_per_pixel_arr
    term_mag         = N_photon_mag_0 * 10**(-m_source_sky/2.5)
    s_source_sky_arr = term_mag * term_area
    return s_source_sky_arr

##########################################################################

def sub_step_7b_compute_N_photon_for_0_magnitude(inp_dict,flag_source_sky="source"):

    if flag_source_sky.lower() == "source":
        wave_band  = inp_dict["wave_band_source"]

    elif flag_source_sky.lower() == "sky":
        wave_band  = inp_dict["wave_band_sky"]

    #lambda_arr_obs = inp_dict["lambda_arr_obs"]
    wave_band_dict = inp_dict["wave_band_dict"]
    lambda_c       = wave_band_dict[wave_band]["lambda_c"]
    flux_jy        = wave_band_dict[wave_band]["flux_jy"]
    # Factor 10000.0 is because N_photons are define in cm^-2 and not m^-2
    #unit_energy    = 1e-23 / h_planck_cgs
    #N_photon_mag_0 = flux_jy * unit_energy/ (lambda_c)
    N_photon_mag_0 = flux_jy * 1e4 * 1e4 / (6.6*lambda_c)

    return N_photon_mag_0

def sub_step_7b_compute_N_photon_for_0_magnitude_backup(inp_dict,flag_source_sky="source"):

    if flag_source_sky.lower() == "source":
        wave_band  = inp_dict["wave_band_source"]

    elif flag_source_sky.lower() == "sky":
        wave_band  = inp_dict["wave_band_sky"]

    wave_band_dict = inp_dict["wave_band_dict"]
    lambda_c       = wave_band_dict[wave_band]["lambda_c"]
    flux_jy        = wave_band_dict[wave_band]["flux_jy"]
    # Factor 10000.0 is because N_photons are define in cm^-2 and not m^-2
    #unit_energy    = 1e-23 / h_planck_cgs
    #N_photon_mag_0 = flux_jy * unit_energy/ (lambda_c)
    N_photon_mag_0 = flux_jy * 1e4 * 1e4 / (6.6*lambda_c)

    return N_photon_mag_0

##########################################################################
##########################################################################
##########################################################################

def step_9_compute_flux_normalization_from_magnitude(inp_dict):
    flux_m0_jy_source,lambda_c_source = step_9a_compute_flux_in_jy_for_given_source(inp_dict,"source")
    flux_m0_jy_sky,lambda_c_sky       = step_9a_compute_flux_in_jy_for_given_source(inp_dict,"sky")

    m_source          = inp_dict["m_source"]
    m_sky             = inp_dict["m_sky"]
    # Flux units are erg / s / cm^2 / Hz
    flux_source_jy    = 10**(-0.4*m_source) * flux_m0_jy_source * 1e-23
    flux_sky_jy       = 10**(-0.4*m_sky) * flux_m0_jy_sky * 1e-23

    # Flux units are erg / s / cm^2 / Angstroms
    factor            = 3e8 / 1e-10 # light speed by Angstrom to meter conversion
    flux_source_Ang   = flux_source_jy * factor / lambda_c_source**2.0
    flux_sky_Ang      = flux_sky_jy * factor / lambda_c_sky**2.0

    flux_norm         = flux_source_Ang + flux_sky_Ang


    inp_dict.update({"flux_norm_(erg/s/cm^2/A)":flux_norm})

    return inp_dict

def step_9a_compute_flux_in_jy_for_given_source(inp_dict,flag_source_sky):
    if flag_source_sky.lower() == "source":
        wave_band  = inp_dict["wave_band_source"]

    elif flag_source_sky.lower() == "sky":
        wave_band  = inp_dict["wave_band_sky"]

    wave_band_dict = inp_dict["wave_band_dict"]
    lambda_c       = wave_band_dict[wave_band]["lambda_c"]
    flux_jy        = wave_band_dict[wave_band]["flux_jy"]
    return flux_jy,lambda_c



##########################################################################
##########################################################################
##########################################################################

def step_10_add_noise_to_spectra(inp_dict):
    inp_dict_Flux       = inp_dict["inp_dict_Flux"]
    Flux_arr_cont_added = inp_dict_Flux["Flux_arr_cont_added"]
    snr_per_pixel_arr   = inp_dict["snr_per_pixel_arr"]
    bool_arr_0          = snr_per_pixel_arr == 0.0
    bool_arr_neg        = snr_per_pixel_arr < 0.0
    bool_arr            = bool_arr_0 | bool_arr_neg
    snr_per_pixel_arr[bool_arr] = np.median(snr_per_pixel_arr)
    noise_arr           = np.random.normal(1.0,1.0 / snr_per_pixel_arr)

    flux_arr_no_noise   = Flux_arr_cont_added[:,1]
    wavelength_arr      = Flux_arr_cont_added[:,0]
    flux_arr            = flux_arr_no_noise * noise_arr # Flux array with noise
    cont_arr            = Flux_arr_cont_added[:,2]

    out_dict_final      = {"wavelength_arr":wavelength_arr,\
                           "flux_arr":flux_arr,\
                           "flux_arr_no_noise":flux_arr_no_noise,\
                           "cont_arr":cont_arr,\
                           "snr_per_pixel_arr":snr_per_pixel_arr}

    inp_dict.update({"out_dict_final":out_dict_final})
    return inp_dict



##########################################################################
############## INSTRUMENT AND DETECTOR PROPERTIES ROUTINES  ##############
##########################################################################

##########################################################################

def initilize_waveband():
    """
    lambda_c -->>  Central wavelength of waveband
    flux_jy  -->>  Flux in jansky in the waveband (1 Jy = 10^-26 W/Hz/m^2)
    exinct   -->>  Extinsion in waveband for airmass = 1.0

    Following table taken from Chris Benn's code
    Central wavelengths of bands used (from output of ~crb/ing/filter.f).
    Typical wavelength ranges are given below.  Actual values of effective
    bandwidth (integral over T(l)dl, where T is transmission, l is wavelength)
    are given in array BAND
                                          BAND values are for:
       U  3600 A  (3300 - 3800 A)         50-mm square and 139-mm liq-CuSO4 U
       B  4300 A  (4000 - 5000 A)         Harris
       V  5500 A  (5000 - 6000 A)         Harris
       R  6500 A  (6000 - 7200 A)         Harris
       I  8200 A  (7200 - 9200 A)         Harris
       Z  9500 A  (8500 - 9700 A)         RGO glass
       J 12500 A  (1.10 - 1.40 mic)       WHIRCAM
       H 16700 A  (1.50 - 1.82 mic)       WHIRCAM
       K 21600 A  (2.00 - 2.30 mic)       WHIRCAM Kshort

    Note that the central wavelengths and wavelength range are different
    for different types of filter.  E.g. the NOT ALFOSC filters have:
       Centre   FWHM Thruput
     U  3620 A   600  60%
     B  4400 A  1000  68%
     V  5300 A   800  85%
     R  6500 A  1300  81%
     i  7970 A  1570  91%

    """
    wave_band_dict = {\
    'U':{"lambda_c":3600.0,  "delta_lambda":500.0,  "flux_jy":1810.0, "sky_mag":22.0, "extinct":0.4554},\
    'B':{"lambda_c":4300.0,  "delta_lambda":720.0,  "flux_jy":4260.0, "sky_mag":22.7, "extinct":0.2175},\
    'V':{"lambda_c":5500.0,  "delta_lambda":860.0,  "flux_jy":3640.0, "sky_mag":21.9, "extinct":0.1020},\
    'R':{"lambda_c":6500.0,  "delta_lambda":1330.0, "flux_jy":3080.0, "sky_mag":21.0, "extinct":0.0575},\
    'I':{"lambda_c":8200.0,  "delta_lambda":1400.0, "flux_jy":2550.0, "sky_mag":20.0, "extinct":0.0157},\
    'Z':{"lambda_c":9500.0,  "delta_lambda":700.0,  "flux_jy":2200.0, "sky_mag":18.8, "extinct":0.0087},\
    #'J':{"lambda_c":12500.0, "delta_lambda":1600.0, "flux_jy":1570.0, "sky_mag":16.1, "extinct":0.10},\
    #'H':{"lambda_c":16700.0, "delta_lambda":2900.0, "flux_jy":1020.0, "sky_mag":14.7, "extinct":0.11},\
    #'K':{"lambda_c":21600.0, "delta_lambda":3200.0, "flux_jy":640.0,  "sky_mag":12.5, "extinct":0.07},\
    #'?':{"lambda_c":0.0,     "delta_lambda":0.0,    "flux_jy":0.0,    "sky_mag":0.0,  "extinct":0.0},\
    }
    return wave_band_dict


##########################################################################
##########################################################################
##########################################################################



