# -*- coding: utf-8 -*-
from typing import Any

import math
import numpy as np

try:
    from .splinefit import Splinefit
except:
    from splinefit import Splinefit

try:
    from .guitastrotools import GuitastroTools, GuitastroException
except:
    from guitastrotools import GuitastroTools, GuitastroException

# #####################################################################
# #####################################################################
# #####################################################################
# Class ExposureTimeCalculator
# #####################################################################
# #####################################################################
# Dictionaries:
# self._optics : D F/D Topt Fwhm_psf_opt proce
# self._cameras :
# self._etc["param"]["optic"]["D"] = [0.3 ,"Optic diameter (m)"]
#
#
# #####################################################################

class ExposureTimeCalculatorException(GuitastroException):
    """Exception raised for errors in the ExposureTimeCalculator class.
    """

    KEY_NOT_FOUND = 0
    CAMERA_NOT_FOUND = 1
    OPTIC_NOT_FOUND = 2

    errors = [""]*3
    errors[KEY_NOT_FOUND] = "Key not found"
    errors[CAMERA_NOT_FOUND] = "Camera not found"
    errors[OPTIC_NOT_FOUND] = "Optic not found"


class ExposureTimeCalculator(ExposureTimeCalculatorException, GuitastroTools):
    """Exposure Time Calculator.

    The Exposure Time Calculator (ETC) take input parameters defining an observatory (optics, camera, sky)
    and compute limiting magnitude or other useful values to simulate images.

    First, init the ETC. For example, we choose a camera "ProLine 16803" and an optical tube "Takahashi_180ED":

    ::

        etc = ExposureTimeCalculator()
        camera = "ProLine 16803"
        etc.camera(camera)
        optic = "Takahashi_180ED"
        etc.optics("Takahashi_180ED")

    Second, the Full Width at Half Maximum (FWHM) of the Point Spread Function (PSF) is 15 micro-meters.
    We choice no filter and a local seeing of 2.5 arcsec.

    ::

        etc.params("Fwhm_psf_opt", 15e-6)
        etc.params("band","C")
        etc.params("seeing", 2.5)

    For a given exposure time of 100 seconds,
    a Signal to Noise Ratio (SNR) of 5, we compute the magnitude:

    ::

        etc.inputs("t", 100)
        etc.inputs("snr", 5)
        m = etc.snr2m_computations()

    The m variable contains the computed magnitude.
    You can print all the detailed computation:

    ::

        etc

    ETC have three important methods according known inputs:

        * **snr2m_computations**: To compute apparent magnitude to reach a given signal to noise ratio and an exposure time.
        * **snr2t_computations**: To compute exposure time to reach a given signal to noise ratio and an apparent magnitude.
        * **t2snr_computations**: To compute signal to noise ration according an exposure time and an apparent magnitude.

    Before using these methods, you have to set the input parameters using method params as showed in the next examples:

    ::

        etc.inputs("t", 100)
        etc.inputs("snr", 5)
        m = etc.snr2m_computations()

    ::

        etc.inputs("m", 18)
        etc.inputs("snr", 5)
        m = etc.snr2t_computations()

    ::

        etc.inputs("m", 18)
        etc.inputs("t", 100)
        m = etc.t2snr_computations()

    Note that t is the time exposure expressed in seconds

    The parameters to simulate an image of this star by a Gaussiuan is:

    ::

        print(etc.simu_star_params())

    To get the list of supported cameras:

    ::

        etc = ExposureTimeCalculator()
        cameras = etc.camera()

    Of course one can define its own camera setup using the method set_params.

    To get the list of supported optics:

    ::

        etc = ExposureTimeCalculator()
        optics = etc.optics()

    Of course one can define its own optics setup using the method set_params.

    """
    def __init__(self):
        self.init()

    def __repr__(self):
        return str(self)

    def _set_array_optics(self):
        #
        # D : Optic diameter (m)
        # F/D : Focal diameter ratio
        # Topt : Transmission of the optics in the photometric band (Reflec=0.8, Refrac=0.95)
        # Fwhm_psf_opt : Fwhm of the point spread function in the image plane (m)
        #
        self.OPT_DIAM = 0 # m
        self.OPT_FOND = 1
        self.OPT_TRANSMISSION = 2
        self.OPT_FWHM_PSF = 3  # m
        self.OPT_FULL_DIAM = 4 # (m) Full light diameter at focus plane
        self.OPT_PRICE = 5 # (euros)
        self._optics = {}
        # ---
        # expr pow(0.8,1)*pow(0.95,2)*(1-pow(0.58,2))
        self._optics["TAROT"] = [0.250, 3.40, 0.48, 10e-6, 40e-3, 15000]
        # expr pow(0.8,1)*pow(0.95,3)*(1-pow(0.45,2))
        self._optics["Takahashi_180ED"] = [0.180, 2.80, 0.55,  5e-6, 44e-3,  5900]
        # expr pow(0.8,1)*pow(0.95,2)*(1-pow(0.5,2))
        self._optics["ASA_200mm_Newton"] = [0.200, 2.80, 0.54,  5e-6, 52e-3,  9000]
        # expr pow(0.8,1)*pow(0.95,2)*(1-pow(0.4,2))
        self._optics["ASA_250mm_Newton"] = [0.250, 3.60, 0.61,  5e-6, 50e-3,  9000]
        self._optics["ASA_H400mm"] = [0.400, 2.40, 0.61,  5e-6, 70e-3, 44000]
        # expr pow(0.8,1)*pow(0.95,2)*(1-pow(0.57,2))
        self._optics["OS_RiFast_300"] = [0.300, 3.80, 0.48,  5e-6, 60e-3, 17500]
        # expr pow(0.8,1)*pow(0.95,2)*(1-pow(0.55,2))
        self._optics["OS_RiFast_400"] = [0.400, 3.80, 0.50,  5e-6, 70e-3, 36000]
        self._optics["OS_RiFast_500"] = [0.500, 3.80, 0.50,  5e-6, 80e-3, 62000]
        # expr pow(0.8,1)*pow(0.95,3)*(1-pow(0.55,2))
        self._optics["OS_RH320_AT"] = [0.320, 2.20, 0.47,  5e-6, 52e-3, 40000]
        self._optics["OS_RH350_AT"] = [0.350, 2.80, 0.47,  5e-6, 60e-3, 40000]
        self._optics["OS_RH350_AT_x0.8"] = [0.350, 2.24, 0.47,  5e-6, 40e-3, 43000]
        # expr pow(0.8,1)*pow(0.95,2)*(1-pow(0.32,2))
        self._optics["HTP-400"] = [0.400, 3.00, 0.65,  5e-6, 60e-3, 20250]
        self._optics["HTP-500"] = [0.500, 3.00, 0.65,  5e-6, 60e-3, 30450]
        # expr pow(0.8,1)*pow(0.95,4)*(1-pow(0.42,2))
        self._optics["Celestron_RASA_11"] = [0.280, 2.20, 0.54,  5e-6, 70e-3,  4400]
        #
        self._optics["GB-400A"] = [0.400, 3.00, 0.53, 10e-6, 52e-3, 80000]
        self._optics["GB-400B"] = [0.400, 2.90, 0.52, 20e-6, 90e-3, 80000]

    def optics(self, typeopt=""):
        """Get/Set the list of supported optics

        ::

            etc = ExposureTimeCalculator()
            optics = etc.optics()

        To print the initialized parameters, see the section "param optic" of:

        ::

            etc

        """
        self._set_array_optics()
        if typeopt == "":
            return self._optics.keys()
        if typeopt in self._optics.keys():
            opt = self._optics[typeopt]
            D, FonD, Topt, Fwhm_psf_opt, diamFull, optPrice = opt
            self._etc["param"]["optic"]["D"][0] = D
            self._etc["param"]["optic"]["FonD"][0] = FonD
            self._etc["param"]["optic"]["Topt"][0] = Topt
            self._etc["param"]["optic"]["Fwhm_psf_opt"][0] = Fwhm_psf_opt
            self._etc["param"]["optic"]["diamFull"][0] = diamFull
            self._etc["param"]["optic"]["optPrice"][0] = optPrice
        else:
            msg = f"The asked {typeopt} is not found amongst {self._optics.keys()}"
            raise ExposureTimeCalculatorException(ExposureTimeCalculatorException.OPTIC_NOT_FOUND, msg)

    def _set_array_cameras(self):
        #---------------------------------------------------------
        #  [list camName  {nbcell1 nbcell2 photocell1 photocell2 C_th G N_ro eta Em}]
        #  - nbcell1     == Number of pixels on axis1
        #  - nbcell2     == Number of pixels on axis2
        #  - photocell1 == Pixel size (m)
        #  - photocell2 == Pixel size (m)
        #  - C_th       == Thermic coefficient (electrons/sec/photocell)
        #  - G          == CCD gain (electrons/ADU)
        #  - N_ro       == Readout noise (electrons)
        #  - eta        == CCD Quantum efficiency in the photometric band (electron/photon)
        #  - Em         == Electron multiplier (>1 if EMCCD, else =1)

        self.CAM_NAXIS1 = 0 # Number of pixels on axis1
        self.CAM_NAXIS2 = 1 # Number of pixels on axis2
        self.CAM_PHOTOCELL1 = 2  #  Pixel size (m)
        self.CAM_PHOTOCELL2 = 3  #  Pixel size (m)
        self.CAM_THERM_COEF = 4  # Thermic coefficient (electrons/sec/photocell)
        self.CAM_GAIN = 5  # CCD gain (electrons/ADU)
        self.CAM_READOUT_NOISE = 6  # Readout noise (electrons)
        self.CAM_QUANTUM_EFF = 7  # Quantum efficiency in the photometric band (electron/photon)
        self.CAM_ELECTRON_AMP = 8  #  Electron multiplier (>1 if EMCCD, else =1)
        self.CAM_PRICE = 9  # Camera price (euros)

        self._cameras = {}
        # --- QHYCCD
        #                                               nax1 nax2   cell1  cell2  Th    G   Noise RQE Em
        self._cameras["QHY 163M"] = [4656, 3522,  3.8e-6,  3.8e-6, 0.01,   1.,   1.8, 0.6,  1,   1600]
        self._cameras["QHY 411"] = [14192, 10640, 3.76e-6, 3.76e-6, 0.005,  0.5,  3.7, 0.95, 1, 200000]
        self._cameras["QHY 42"] = [2048, 2048,   11e-6,   11e-6,  0.05,  1,    1.7, 0.95, 1,  11000]
        self._cameras["QHY 600U3G20M"] = [9576, 6388, 3.76e-6, 3.76e-6, 0.005,  0.5,  3.7, 0.95, 1,   5000]
        self._cameras["QHY 45GX"] = [2084, 2085,   24e-6,   24e-6,   0.3,  1,   25,   0.73, 1,  40000]
        self._cameras["QHY 50GX"] = [8176, 6132,    6e-6,    6e-6, 0.003,  1.5, 11.5, 0.2,  1,  20000]

        # --- ZWO ASI 6200MM
        #                                               nax1 nax2   cell1  cell2  Th    G   Noise RQE Em
        self._cameras["ZWO 6200MM"] = [9576, 6388, 3.76e-6, 3.76e-6, 0.005,  0.5,  3.7, 0.95, 1,  5000]

        # --- Andor
        self._cameras["Andor DW436"] = [2048, 2048, 13.5e-6, 13.5e-6, 0.045,   2.8,   9.2, 0.85, 1,  50000] ; # measured on TAROT at 1MHz
        self._cameras["Andor Neo sCMOS"] = [2560, 2160,  6.5e-6,  6.5e-6, 0.07,    4.,    1.4, 0.57, 1,  30000]
        self._cameras["Andor Lucas R DL-604"] = [1004, 1002,  8.0e-6,  8.0e-6,   0.07,    3.,   18, 0.65, 1,  30000]
        self._cameras["Andor iKon-L 936 Z-BV"] = [2048, 2048, 13.5e-6, 13.5e-6, 0.00013, 2.8,   2.9, 0.9,  1,  65000] ; # 20s redout time
        self._cameras["Andor DW936 BV"] = [2048, 2048, 13.5e-6, 13.5e-6, 0.045,   1.8,   6.8, 0.85, 1,  65000] ; # measured on Zadko at 1 MHz
        self._cameras["Andor Marana sCMOS"] = [2048, 2048, 11.0e-6, 11.0e-6, 0.2,     1.5,   1.6, 0.95, 1,  50000]
        self._cameras["Andor Balor sCMOS"] = [4129, 4104, 12.0e-6, 12.0e-6, 0.08,    1.5,   3.0, 0.61, 1, 200000]

        # --- Audine
        self._cameras["Audine Kaf401ME"] = [ 768,  512,  9e-6,    9e-6,   0.2,   2.1,  12,   0.5,  1,  1000]
        self._cameras["Audine Kaf1600"] = [1536, 1024,  9e-6,    9e-6,   0.2,   2.1,  12,   0.5,  1,  2000]

        # --- SBig
        self._cameras["ST-402ME"] = [765,  510,  9e-6,    9e-6,   0.1,   1.5,  13.8, 0.75, 1,  1000]
        self._cameras["ST-1603ME"] = [1530, 1020,  9e-6,    9e-6,   0.1,   1.7,  18.0, 0.75, 1,  2000]
        self._cameras["ST-3200ME"] = [2184, 1472,  6.8e-6,  6.8e-6, 0.1,   1.0,  10.0, 0.75, 1,  3000]
        #
        self._cameras["STF-8300M"] = [3326, 2504,  5.4e-6,  5.4e-6, 0.05,  1.0,  15.0, 0.45, 1,  2000]
        self._cameras["STT-8300M"] = [3326, 2504,  5.4e-6,  5.4e-6, 0.02,  0.37,  9.3, 0.45, 1,  2000]
        self._cameras["STT-1603ME"] = [1536, 1024,  9e-6,    9e-6,   0.1,   2.3,  15.0, 0.75, 1,  3000]
        self._cameras["STT-3200ME"] = [2184, 1472,  6.8e-6,  6.8e-6, 0.06,  1.0,  10.0, 0.75, 1,  4000]
        #
        self._cameras["STXL-11002"] = [4008, 2672,  9e-6,    9e-6,   0.07,  0.80, 12.0, 0.45, 1, 11000] ; # measured Les Makes bin1x1 -20°C
        self._cameras["STXL-6303E"] = [3072, 2048,  9e-6,    9e-6,   0.3,   1.47, 11,   0.65, 1,  9000]
        #
        self._cameras["STX-16803"] = [4096, 4096,  9e-6,    9e-6,   0.02,  1.27, 10,   0.6,  1, 12000]
        #
        self._cameras["ST-2000XM"] = [1600, 1200,  7.4e-6,  7.4e-6, 0.1,   0.6,   7.9, 0.35, 1,  2000]

        # --- Princeton
        self._cameras["Peregrine 486"] = [4096, 4096, 15e-6,   15e-6,   0.01,  1.8,  10,   0.9,  1, 130000]

        # --- Apogee
        self._cameras["Alta F230"] = [2048, 2048, 15e-6,   15e-6,   0.4,   1.5,  12,   0.85, 1, 60000]
        self._cameras["Alta F42"] = [2048, 2048, 13.5e-6, 13.5e-6, 1,     1.5,   9,   0.85, 1, 50000]
        #
        self._cameras["Alta F16M"] = [4096, 4096,  9e-6,    9e-6,   0.2,   1.5,   9,   0.6,  1,  13000]
        self._cameras["Alta F9000"] = [3056, 3056, 12e-6,   12e-6,   0.6,   1.29, 11,   0.6,  1,  12000] ; # measured. readout time 1MHz -35°C
        self._cameras["Alta F4320"] = [2048, 2048, 24e-6,   24e-6,   2,     1.5,  12,   0.6,  1,  30000]
        self._cameras["Alta F6"] = [1024, 1024, 24e-6,   24e-6,   0.5,   1.5,   8,   0.6,  1,  20000]
        #
        self._cameras["Alta F16000"] = [4872, 3248,  7.4e-6,  7.4e-6, 0.01,  1.5,   9,   0.4,  1,  20000]
        self._cameras["Alta F29050"] = [6576, 4384,  5.5e-6,  5.5e-6, 0.15,  1.5,   6,   0.4,  1,  20000]

        # --- FLI
        self._cameras["ProLine PL09000"] = [3056, 3056, 12e-6,   12e-6,   0.1,   1.5,  10,   0.65, 1,  12000]
        self._cameras["ProLine 16801"] = [4096, 4096,  9e-6,    9e-6,   0.08,  1.5,   9,   0.6,  1,  13000] ; # readout time 8 MHz -35°C
        self._cameras["ProLine 16803"] = [4096, 4096,  9e-6,    9e-6,   0.005, 1.5,  10,   0.6,  1,  13000]
        self._cameras["ProLine 4301"] = [2084, 2084, 24e-6,   24e-6,   0.4,   1.5,   8,   0.6,  1,  20000]
        #
        self._cameras["ProLine 230 Midband"] = [2048, 2048, 15e-6,   15e-6,   0.4,   1.5,   9.5, 0.85, 1,  50000]
        self._cameras["ProLine 3041 Broadband"] = [2048, 2048, 15e-6, 15e-6,   0.3,   1.5,   8,   0.85, 1,  50000]
        self._cameras["ProLine 4240 Midband"] = [2048, 2048, 13.5e-6, 13.5e-6, 0.2,   1.5,   8,   0.85, 1,  45000]
        self._cameras["ProLine 4720"] = [1024, 1024, 13e-6,   13e-6,   0.02,  1.5,  10,   0.85, 1,  30000]
        self._cameras["ProLine 4710 Deep D."] = [1024, 1024, 13e-6,   13e-6,   7,     1.5,  10,   0.9,  1,  30000]
        #
        self._cameras["MicroLine 50100"] = [8176, 6132,  6e-6,    6e-6,   0.003, 1.5,  11.5, 0.2,  1,  24000]
        #
        self._cameras["Kepler 400"] = [2048, 2048, 11e-6,   11e-6,   0.6,   1.5,   1.6, 0.95, 1,  26000] ; # + GPS
        self._cameras["Kepler 4040"] = [4096, 4096,  9e-6,    9e-6,   0.15,  1.0,   3.7, 0.74, 1,  20500] ; # + GPS
        self._cameras["Kepler 6060 BSI"] = [6144, 6200, 10e-6,   10e-6,   0.6,   1.5,   3.0, 0.95, 1, 134000] ; # grade 2 (189000 € for grade 1)
        self._cameras["Kepler 6060 FSI"] = [6144, 6200, 10e-6,   10e-6,   0.6,   1.5,   4.6, 0.72, 1,  56000] ; # grade 1 (107000 € for grade 0)

        # --- Atik
        self._cameras["Atik Titan"] = [659, 494,   7.4e-6,  7.4e-6,   0.2,   1.5,   5,   0.6,  1,  635]
        self._cameras["Atik 314EX"] = [1392, 1040, 4.65e-6, 4.65e-6,  0.2,   1.5,   3,   0.6,  1, 1429]
        self._cameras["Atik 320E"] = [1620, 1220, 4.40e-6, 4.40e-6,  0.2,   1.5,   3,   0.6,  1, 1040]
        self._cameras["Atik 314L+"] = [1392, 1040, 6.45e-6, 6.45e-6,  5e-4, 0.267,  3.7, 0.6,  1, 1347] ; # T=0°C http://www.dangl.at/ausruest/atik_314/atik_314_e.htm
        self._cameras["Atik 383L+"] = [3362, 2504, 5.40e-6, 5.40e-6,  0.2,   1.5,   7,   0.6,  1, 1999]
        self._cameras["Atik 420"] = [1620, 1220, 4.40e-6, 4.40e-6,  0.1,   1.5,   4,   0.6,  1, 1045]
        self._cameras["Atik 450"] = [2448, 2050, 3.45e-6, 3.45e-6,  0.1,   1.5,   5,   0.6,  1, 2356]
        self._cameras["Atik 428EX"] = [1932, 1452, 4.54e-6, 4.54e-6,  0.1,   1.5,   5,   0.6,  1, 1771]
        self._cameras["Atik 460EX"] = [2750, 2200, 4.54e-6, 4.54e-6,  0.1,   1.5,   5,   0.6,  1, 2356]
        self._cameras["Atik 490EX"] = [3380, 2704, 3.69e-6, 3.69e-6,  0.1,   1.5,   5,   0.6,  1, 2630]
        self._cameras["Atik 4000LE"] = [2048, 2048, 7.4e-6,  7.4e-6,   0.01,  1.5,  11,   0.6,  1, 4026]
        self._cameras["Atik 11000"] = [4008, 2672, 9e-6,    9e-6,     0.03,  1.5,  13,   0.6,  1, 5483]

        # --- Raptor Photonics
        self._cameras["OWL VIS-SWIR 320 0deg"] = [320,  256, 30e-6,   30e-6,    264000,   6.55,  131,   0.8, 1, 15000] ; # T=0°C
        self._cameras["OWL VIS-SWIR 320 -40deg"] = [320,  256, 30e-6,   30e-6,    131,      6.55,  131,   0.8, 1, 15000] ; # T=-40°C
        self._cameras["OWL SWIR 640 0deg"] = [640,  512, 15e-6,   15e-6,    264000/4.,   6.55,  131/4.,   0.75, 1, 60000] ; # T=0°C
        self._cameras["OWL SWIR 640 -40deg"] = [640,  512, 15e-6,   15e-6,    131/4.,    6.55,  131/4.,   0.75, 1, 60000] ; # T=-40°C

        #--   Point Grey U3->USB3 Capteurs CCD Sony
        self._cameras["GS-U3-28S4M ICX687"] = [1928,  1448, 3.69e-6,   3.69e-6,   0.79,   0.16,  11.01,   0.71, 1, 4000] ; # WellDepth=9387 e-
        self._cameras["GS-U3-28S5M ICX674"] = [1920,  1440, 4.54e-6,   4.54e-6,   1.27,   0.24,   9.39,   0.67, 1, 4000] ; # WellDepth=14693 e-
        self._cameras["GS-U3-60S6M ICX694"] = [2736,  2192, 4.54e-6,   4.54e-6,   0.82,   0.23,  10.54,   0.73, 1, 4000] ; # WellDepth=14446 e-
        self._cameras["GS-U3-91S6M ICX814"] = [3376,  2704, 3.69e-6,   3.69e-6,   0.42,   0.16,   9.43,   0.75, 1, 4000] ; # WellDepth=9996 e-

        # --- Spectral instruments
        self._cameras["Spectral Instruments 230-84 500KHz"] = [4112, 4096, 15e-6, 15e-6, 0.01,  1.5,  7.0, 0.92, 1, 130000]
        self._cameras["Spectral Instruments 230-84 2MHz"] = [4112, 4096, 15e-6, 15e-6, 0.01,  1.5, 14.7, 0.92, 1, 130000]
        self._cameras["Spectral Instruments 231-84 100KHz"] = [4112, 4096, 15e-6, 15e-6, 0.0003,  1.5,  2.1, 0.92, 1, 130000]
        self._cameras["Spectral Instruments 231-84 1.5MHz"] = [4112, 4096, 15e-6, 15e-6, 0.0003,  1.5,  9.7, 0.92, 1, 130000]

        # --- Watec
        self._cameras["Watec 120N+"] = [720, 576, 8.6e-6, 8.3e-6, 10,  5,  30, 0.7, 1, 250]

        # --- https://lytid.com/case-study/astrophysical-observations-with-siris/ Lytid SIRIS
        self._cameras["Lytid SIRIS"] = [640, 480, 15e-6, 15e-6, 1.0,  6.0,  10.0, 0.9, 1, 110000]

    def camera(self, typecam=""):
        """Get/Set the list of supported cameras

        ::

            etc = ExposureTimeCalculator()
            cameras = etc.camera()

        To print the initialized parameters, see the section "param ccd" of:

        ::

            etc

        """
        self._set_array_cameras()
        if typecam == "":
            return self._cameras.keys()
        if typecam in self._cameras.keys():
            cam = self._cameras[typecam]
            nbcell1, nbcell2, photocell1, photocell2, C_th, G, N_ro, eta, Em, camPrice = cam
            self._etc["param"]["ccd"]["nbcell1"][0] = nbcell1
            self._etc["param"]["ccd"]["nbcell2"][0] = nbcell2
            self._etc["param"]["ccd"]["photocell1"][0] = photocell1
            self._etc["param"]["ccd"]["photocell2"][0] = photocell2
            self._etc["param"]["ccd"]["eta"][0] = eta
            self._etc["param"]["ccd"]["N_ro"][0] = N_ro
            self._etc["param"]["ccd"]["C_th"][0] = C_th
            self._etc["param"]["ccd"]["G"][0] = G
            self._etc["param"]["ccd"]["Em"][0] = Em
            self._etc["param"]["ccd"]["camPrice"][0] = camPrice
        else:
            msg = f"The asked {typecam} is not found amongst {self._cameras.keys()}"
            raise ExposureTimeCalculatorException(ExposureTimeCalculatorException.CAMERA_NOT_FOUND, msg)

    def _params_defaults(self, band="V", moon_age=0.0):
        #---------------------------------------------------------
        #  brief set default input parameters
        #  remark called by init
        #  param band (default) V
        #  param moon_age (default) 0
        #  private
        #

       self._etc["param"] = {}
       self._etc["param"]["local"] = {}
       self._etc["param"]["object"] = {}
       self._etc["param"]["optic"] = {}
       self._etc["param"]["ccd"] = {}
       self._etc["param"]["filter"] = {}

       self._etc["param"]["local"]["moon_age"] = [moon_age , "Age of the Moon (day)", "MOON_AGE"]
       self._etc["param"]["object"]["band"] = [band , "Photometric system symbol", "PHOTBAND"]
       self._modify_band()
       Tatm0 = self._params_Tatm0(300, 0.10)
       self._etc["param"]["local"]["Tatm0"] = [Tatm0 , "Zenith transmission of the atmosphere in the photometric band", "ZENTRANS"]
       self._etc["param"]["local"]["Elev"] = [65 , "Elevation above horizon (deg)", "ELEV"]
       self._etc["param"]["local"]["seeing"] = [3.0 , "Fwhm of the seeing (arcsec)", "SEEING"]

       self._etc["param"]["optic"]["D"] = [0.3 ,"Optic diameter (m)", "TELDIAM"]
       self._etc["param"]["optic"]["FonD"] = [4.0 ,"Focal diameter ratio", "FOND"]
       self._etc["param"]["optic"]["Topt"] = [0.8*0.8*0.95*0.95 ,"Transmission of the optics in the photometric band (Reflec=0.8, Refrac=0.95)", "OPTTRANS"]
       self._etc["param"]["optic"]["Fwhm_psf_opt"] = [15e-6 ,"Fwhm of the point spread function in the image plane (m)", "OPTFWHM"]
       self._etc["param"]["optic"]["diamFull"] = [40e-3 ,"Full light diameter at image plane (m)", "FULLDIAM"]
       self._etc["param"]["optic"]["optPrice"] = [10000 ,"Optical price (euro)", "OPTPRICE"]

       self._etc["param"]["ccd"]["nbcell1"] = [2048 ,"Number of photocells on an axis1", "NBCELL1"]
       self._etc["param"]["ccd"]["nbcell2"] = [2048 ,"Number of photocells on an axis2", "NBCELL2"]
       self._etc["param"]["ccd"]["photocell1"] = [13.5e-6 ,"Photocell size on axis1 (m)", "CELLDIM1"]
       self._etc["param"]["ccd"]["photocell2"] = [13.5e-6 ,"Photocell size on axis2 (m)", "CELLDIM2"]
       self._etc["param"]["ccd"]["bin1"] = [1 ,"Binning on axis1 (photocells/pixel)", "BIN1"]
       self._etc["param"]["ccd"]["bin2"] = [1 ,"Binning on axis2 (photocells/pixel)", "BIN2"]
       self._etc["param"]["ccd"]["eta"] = [0.9 ,"CCD Quantum efficiency in the photometric band (electron/photon)", "CCD_ETA"]
       self._etc["param"]["ccd"]["N_ro"] = [8.5, "Readout noise (electrons/pixel)", "CCD_RON"]
       self._etc["param"]["ccd"]["C_th"] = [0.002, "Thermic coefficient (electrons/sec/photocell)", "CCD_TRM"]
       self._etc["param"]["ccd"]["G"] = [1.8, "CCD gain (electrons/ADU)", "CCD_GAIN"]
       self._etc["param"]["ccd"]["Em"] = [1.0, "Electron multiplier (>1 if EMCCD, else =1)", "CCD_EM"]
       self._etc["param"]["ccd"]["camPrice"] = [3000, "Camera price (euro)", "CAMPRICE"]


    def _params_Tatm0(self, altitude_m=300, Aerosol_Optical_Depth=0.10):
        #---------------------------------------------------------
        #  brief set zenith transmission of the atmosphere in the photometric band
        #  param altitude altitude (m) 300 par défaut
        #  param Aerosol_Optical_Depth 0.10 par défaut
        #  return Tatm0 value
        #
        z = 1.
        h = altitude_m*1e-3
        AOD = Aerosol_Optical_Depth ; # 0.07=hiver 0.21=ete : Aerosol Optical Depth (AOD)
        wvl = self._etc["param"]["filter"]["l"][0]
        convert_unit = 1.
        lmu = 1.*wvl*convert_unit
        n1n1 = 0.23465+(1.076e2/(146-1/lmu/lmu))+(0.93161/(41-1/lmu/lmu))
        Ar = 9.4977e-3*math.pow(lmu,-4)*n1n1*n1n1*math.exp(-h/7.996)
        Tz = math.exp(-2*0.0168*math.exp(-15*abs(lmu-0.59)))
        Ao = -2.5*math.log10(Tz)
        Ao0 = 1.5*math.exp(-math.pow((lmu-0.300)/0.012,2))
        Ao1 = 0.03/(1+math.pow((lmu-0.59)/0.07,2))
        Ao2 = 0.011/(1+math.pow((lmu-0.576)/0.006,2))
        Ao3 = 0.009/(1+math.pow((lmu-0.604)/0.01,2))
        Ao4 = 0.01/(1+math.pow((lmu-0.630)/0.013,2))
        Ao5 = 0.0048/(1+math.pow((lmu-0.531)/0.006,2))
        Ao6 = 0.003/(1+math.pow((lmu-0.545)/0.009,2))
        Ao7 = 0.003/(1+math.pow((lmu-0.564)/0.01,2))
        Ao8 = 0.004/(1+math.pow((lmu-0.572)/0.01,2))
        Ao9 = 0.003/(1+math.pow((lmu-0.506)/0.003,2))
        Ao10 = 0.004/(1+math.pow((lmu-0.477)/0.0025,2))
        Ao = Ao0+Ao1+Ao2+Ao3+Ao4+Ao5+Ao6+Ao7+Ao8+Ao9+Ao10
        Aa = 2.5*math.log10( math.exp(AOD*math.pow(lmu/0.55,-1.3)))
        AA = Ar+Ao+Aa
        TT = pow(10,-0.4*AA*z)
        return TT

    def _params_msky(self, band="V", moon_age=0):
        #---------------------------------------------------------
        #  brief set msky parameter from band and moon_age
        #  code _params_msky 20
        #  endcode
        #  param band (default) V
        #  param moon_age (default) 0
        #
        value = 20.0
        if moon_age<=1.5:
           if band=="U":
               value = 22.0
           if band=="B":
               value = 22.7
           if (band=="V") or (band=="g"):
               value = 21.8
           if band=="C":
               value = 21.4
           if (band=="R") or (band=="r"):
               value = 20.9
           if (band=="I") or (band=="i"):
                value = 19.9
        elif moon_age<=5:
           if band=="U":
               value = 21.5
           if band=="B":
               value = 22.4
           if (band=="V") or (band=="g"):
               value = 21.7
           if band=="C":
               value = 21.3
           if (band=="R") or (band=="r"):
               value = 20.8
           if (band=="I") or (band=="i"):
               value = 19.9
        elif moon_age<=8.5:
           if band=="U":
               value = 19.9
           if band=="B":
               value = 21.6
           if (band=="V") or (band=="g"):
               value = 21.4
           if band=="C":
               value = 21.0
           if (band=="R") or (band=="r"):
               value = 20.6
           if (band=="I") or (band=="i"):
               value = 19.7
        elif moon_age<=12:
           if band=="U":
               value = 18.5
           if band=="B":
               value = 20.7
           if (band=="V") or (band=="g"):
               value = 20.7
           if band=="C":
               value = 20.5
           if (band=="R") or (band=="r"):
               value = 20.3
           if (band=="I") or (band=="i"):
               value = 19.5
        else:
           if band=="U":
               value = 17.0
           if band=="B":
               value = 19.5
           if (band=="V") or (band=="g"):
               value = 20.0
           if band=="C":
               value = 19.9
           if (band=="R") or (band=="r"):
               value = 19.9
           if (band=="I") or (band=="i"):
               value = 19.2
           if band=="z":
               value = 17.0
           if band=="J":
               value = 15.7
           if band=="H":
               value = 14.1
           if band=="K":
               value = 13.0
        self._etc["param"]["local"]["msky"] = [value, "Sky brightness in the $band band at moon age $moon_age day (mag/arcsec2)", "SKYMAG"]
        return ""

    def _params_filter(self, band="V"):
        #---------------------------------------------------------
        #  brief set a filter from band
        #  remark called by _modify_band
        #  param band (default) V
        #  private
        #
        valid = True
        if band=="C":
            l = 0.6
            Dl = 0.3
            Fm0 = 3100
        elif band=="U":
            l = 0.36
            Dl = 0.15*l
            Fm0 = 1810
        elif band=="B":
            l = 0.44
            Dl = 0.22*l
            Fm0 = 4260
        elif band=="V":
            l = 0.55
            Dl = 0.16*l
            Fm0 = 3640
        elif band=="R":
            l = 0.64
            Dl = 0.23*l
            Fm0 = 3080
        elif band=="I":
            l = 0.79
            Dl = 0.19*l
            Fm0 = 2550
        elif band=="J":
            l = 1.26
            Dl = 0.16*l
            Fm0 = 1600
        elif band=="H":
            l = 1.60
            Dl = 0.23*l
            Fm0 = 1080
        elif band=="K":
            l = 2.22
            Dl = 0.23*l
            Fm0 = 670
        elif band=="g":
            l = 0.52
            Dl = 0.14*l
            Fm0 = 3730
        elif band=="r":
            l = 0.67
            Dl = 0.14*l
            Fm0 = 4490
        elif band=="i":
            l = 0.79
            Dl = 0.16*l
            Fm0 = 4760
        elif band=="z":
            l = 0.91
            Dl = 0.13*l
            Fm0 = 4810
        elif band=="CAB":
            lmin = 380e-9
            lmax = 860e-9
            fmin = 3e8/lmax
            fmax = 3e8/lmin
            fmean = math.sqrt(fmin*fmax)
            lmean = 3e8/fmean
            dl = lmax-lmin
            l = lmean
            Dl = dl
            Fm0 = 3631
            band = "C"
        else:
            valid = False
            #error "Filter $band not found"
        if valid == True:
            l_comment = "Central wavelength of the filter (micrometers)"
            Dl_comment = "Wavelength bandpass of the filter (micrometers)"
            Fm0_comment = "Flux for magnitude zero for the filter (Jy)"
            self._etc["param"]["filter"]["l"] = [l, l_comment, "FILT_WV"]
            self._etc["param"]["filter"]["Dl"] = [Dl, Dl_comment, "FILT_DWV"]
            self._etc["param"]["filter"]["Fm0"] = [Fm0, Fm0_comment, "FILT_JY0"]

    def inputs(self, key: str="", val: Any="?") -> Any:
        """Get/Set an input according its key.

        Settings of the ETC are "parameters" and "inputs".
        **Parameters** are settings of hardware and siteobs conditions (diameter of the optics, seeing, etc.).
        **Inputs** are settings related to the target or to observing conditions (exposure time, distance of the source, signal to noise ration, etc.)

        Args:
            key: Key of the input. If key is not given (or equal "") it returns the list of available keys

        Returns:
            Value of the input or list available keys if key is not given (or equal "").

        ::

            etc = ExposureTimeCalculator()
            key_list = etc.inputs()
        """
        #---------------------------------------------------------
        ## @brief get input list
        #  @code Exemples :
        #  - inputs m=12.5
        #  - inputs t=120
        #  @endcode
        #  @param args list of {key value}
        #
        keyfound = False
        allkeys = []
        jkeys = self._etc["input"].keys()
        for jkey in jkeys:
            ikeys = self._etc["input"][jkey].keys()
            allkeys.extend(ikeys)
            if key in ikeys:
                #print(f"INPUT {jkey}/{key} val={val}")
                if val == "?":
                    return self._etc["input"][jkey][key]
                keyfound = True
                if jkey=="object":
                    self._etc["input"][jkey][key][0] = val
                    delta = 5.0 * math.log10( self._etc["input"][jkey]["DL_pc"][0] / 10.)
                    if key == "M":
                        self._etc["input"][jkey]["m"][0] = self._etc["input"][jkey]["M"][0] + delta
                    elif key == "m":
                        self._etc["input"][jkey]["M"][0] = self._etc["input"][jkey]["m"][0] - delta
                    elif key == "DL_pc":
                        self._etc["input"][jkey]["m"][0] = self._etc["input"][jkey]["M"][0] + delta
                else:
                    self._etc["input"][jkey][key][0] = val
                res = self._etc["input"][jkey][key]
        if keyfound == False:
            if key=="":
                return allkeys
            msg = f"The asked key {key} is not found amongst {allkeys}"
            raise ExposureTimeCalculatorException(ExposureTimeCalculatorException.KEY_NOT_FOUND, msg)
        return res


    def _inputs_defaults(self):
        #---------------------------------------------------------
        #  brief set default parameters
        #  remark called by init
        #  private
        #
        band = self._etc["param"]["object"]["band"][0]
        self._etc["input"] = {}
        self._etc["input"]["object"] = {}
        self._etc["input"]["ccd"] = {}
        self._etc["input"]["constraint"] = {}

        self._etc["input"]["object"]["M"] = [-21, f"Absolute stellar magnitude in the {band} band", "MAG_ABS"]

        self._etc["input"]["object"]["DL_pc"] = [40e6, "Distance luminosity (pc)", "DIST_PC"]

        m = self._etc["input"]["object"]["M"][0] + 5. * math.log10 (self._etc["input"]["object"]["DL_pc"][0] / 10.)
        self._etc["input"]["object"]["m"] = [m, f"Apparent stellar magnitude in the {band} band", "MAG_APP"]

        self._etc["input"]["ccd"]["t"] = [30, "Exposure time (sec)", "EXPTIME"]

        self._etc["input"]["constraint"]["snr"] = [5, "SNR constrained", "ETC_SNR"]

    def params(self, key: str="", val: Any="?") -> Any:
        """Get/Set a parameter according its key.

        Settings of the ETC are "parameters" and "inputs".
        **Parameters** are settings of hardware and siteobs conditions (diameter of the optics, seeing, etc.).
        **Inputs** are settings related to the target or to observing conditions (exposure time, distance of the source, signal to noise ration, etc.)

        Args:
            key: Key of the parameter. If key is not given (or equal "") it returns the list of available keys

        Returns:
            Value of the parameter or list available keys if key is not given (or equal "").

        ::

            etc = ExposureTimeCalculator()
            key_list = etc.params()
        """
        band = self._etc["param"]["object"]["band"][0]

        keyfound = False
        allkeys = []
        jkeys = self._etc["param"].keys()
        for jkey in jkeys:
            ikeys = self._etc["param"][jkey].keys()
            allkeys.extend(ikeys)
            if key in ikeys:
                #print(f"PARAMS {jkey}/{key} val={val}")
                if val == "?":
                    return self._etc["param"][jkey][key]
                keyfound = True
                if jkey=="object" or jkey=="local":
                    self._etc["param"][jkey][key][0] = val
                    if key == "moon_age":
                        self._etc["param"][jkey][key][0] = val
                        self._params_msky(band, val)
                    elif key == "band":
                        bands = ["B", "C", "H", "I", "J", "K", "R", "U", "V", "g", "r", "i", "z", "CAB"]
                        if val in bands:
                            self._etc["param"][jkey][key][0] = val
                            self._modify_band(val)
                else:
                    self._etc["param"][jkey][key][0] = val
                res = self._etc["param"][jkey][key]
        if keyfound == False:
            if key=="":
                return allkeys
            msg = f"The asked key {key} is not found amongst {allkeys}"
            raise ExposureTimeCalculatorException(ExposureTimeCalculatorException.KEY_NOT_FOUND, msg)
        return res


    def _modify_band(self, band: str="V"):
        """Modify the band pass (filter) of the optics.

        Args:
            band: Symbol of the filter.

        """
        #---  calcul des cofficients L, Dl et Fm0
        self._params_filter(band)

        #--- calcul de Tatm0
        Tatm0 = self._params_Tatm0(300, 0.10)
        self._etc["param"]["local"]["Tatm0"] = [Tatm0 , "Elevation above horizon (deg)", "ELEV"]

        #---  calcul de msky
        self._params_msky(band, self._etc["param"]["local"]["moon_age"][0])

        #---  modification des commentaires
        try:
            self._etc["input"]["object"]["M"][1] = f"Absolute stellar magnitude in the {band} band"
            self._etc["input"]["object"]["m"][1] = f"Apparent stellar magnitude in the {band} band"
        except:
            pass

    def _spec(self, xy, **kwargs):
        # x, y = etc._spec([[400, 407, 550, 900, 935], [0, 20, 95, 30, 0]], method="splinefit", dy=5, s=0.1, limits = [0,100])
        # input
        # y
        # [ [x1, y1], [x2, y2], ...]
        # [ [x1, x2, ...], [y1, y2, ...]]
        # output
        config = {}
        config["method"] = "lininterp" # "splinefit"
        config["limits"] = [] # [0, 1]
        config["dy"] = 0 # len(xs)
        config["s"] = 0 # len(xs)
        lambd1= 350 ; lambd2= 3000; dlamb = 10 ; n = 1+int((lambd2-lambd1)/dlamb)
        config["xx"] = np.linspace(lambd1, lambd2, n)
        for key, val in kwargs.items():
            if key in config.keys():
                config[key] = val
        xxs = config["xx"]
        # ---
        if isinstance(xy, (int,float)) == True:
            yys = xxs*0 + xy
        else:
            xya = np.array(xy)
            s = xya.shape
            ns = len(s)
            if ns==1:
                yys = xxs*0 + xy[0]
                return xxs, yys
            n1 = s[0]
            n2 = s[1]
            if n1 == 2:
                # xy = [ [1, 2, 3], [5, 4, 6] ] => x = [1, 2, 3] y = [5, 4, 6]
                pass
            elif n2 == 2:
                # xy = [ [1, 5], [2, 4], [3, 6] ] => x = [1, 2, 3] y = [5, 4, 6]
                xya = xya.T
            xs = xya[0]
            ys = xya[1]
            k = 0
            for x, y in zip(xs, ys):
                if k == 0:
                    mini = [x, y]
                    maxi = [x, y]
                else:
                    if x < mini[0]:
                        mini = [x, y]
                    if y < maxi[0]:
                        maxi = [x, y]
                k += 1
            xx_mini = np.min(xxs)
            xx_maxi = np.max(xxs)
            if mini[0] < xx_mini :
                xs = np.append(xs,mini[0])
                ys = np.append(xs,mini[1])
            if mini[0] > xx_maxi :
                xs = np.append(xs,maxi[0])
                ys = np.append(xs,maxi[1])
            # ---
            inds = np.argsort(xs)
            xs = np.array([xs[i] for i in inds])
            ys = np.array([ys[i] for i in inds])
            # ---
            n = len(xs)
            nn = len(xxs)
            yys = np.zeros(nn)
            if config["method"] == "lininterp":
                kk = -1
                kfast = 0
                for xx in xxs:
                    kk += 1
                    for k in range(kfast, n-1):
                        x1 = xs[k]
                        x2 = xs[k+1]
                        if xx >= x1 and xx <= x2:
                            y1 = ys[k]
                            y2 = ys[k+1]
                            dy = y2-y1
                            dx = x2-x1
                            if dy == 0 or dx == 0:
                                yy = y1
                            else:
                                yy = y1 + dy * (xx-x1)/dx
                            yys[kk] = yy
                            kfast = k
                            break
            elif config["method"] == "splinefit":
                fiter = Splinefit()
                dy = config["dy"]
                s = config["s"]
                if s==0:
                    s = 1*len(xs)
                yys = fiter.fit(xs, ys, dy, s, xxs)
            # ---
            limits = config["limits"]
            nl = len(limits)
            if nl==2:
                y1 = limits[0]
                y2 = limits[1]
                for k in range(nn):
                    yy = yys[k]
                    if yy < y1:
                        yys[k] = y1
                    if yy > y2:
                        yys[k] = y2
        return xxs, yys

    def _preliminary_computations (self):
        """Preliminary common computations before resolving equations for SNR, exposure time or magnitude.

        This method is called automatically when needed.
        """
        self._etc["comp1"] = {}

        # --- Optics
        FonD = self._etc["param"]["optic"]["FonD"][0]
        D = self._etc["param"]["optic"]["D"][0]
        Foclen = FonD * D
        self._etc["comp1"]["Foclen"] = [Foclen, "Focal length (m)", "FOCLEN"]

        photocell1 = self._etc["param"]["ccd"]["photocell1"][0]
        bin1 = self._etc["param"]["ccd"]["bin1"][0]
        pixsize1 = photocell1 * bin1
        self._etc["comp1"]["pixsize1"] = [pixsize1, "Pixel length on axis1 (m)", "PIXSIZE1"]

        photocell2 = self._etc["param"]["ccd"]["photocell2"][0]
        bin2 = self._etc["param"]["ccd"]["bin2"][0]
        pixsize2 = photocell2 * bin2
        self._etc["comp1"]["pixsize2"] = [pixsize2, "Pixel length on axis2 (m)", "PIXSIZE2"]

        photocell1 = self._etc["param"]["ccd"]["photocell1"][0]
        bin1 = self._etc["param"]["ccd"]["bin1"][0]
        pixsize1 = photocell1 * bin1
        self._etc["comp1"]["pixsize1"] = [pixsize1, "Pixel length on axis1 (m)"]

        cdelt1 = 2 * math.atan ( pixsize1 / Foclen / 2.) * 180. / math.pi * 3600.
        self._etc["comp1"]["cdelt1"] = [cdelt1, "Pixel spatial sampling on axis1 (arcsec/pix)", "CDELT1"]

        cdelt2 = 2 * math.atan ( pixsize2 / Foclen / 2.) * 180. / math.pi * 3600.
        self._etc["comp1"]["cdelt2"] = [cdelt2, "Pixel spatial sampling on axis2 (arcsec/pix)", "CDELT2"]

        W = cdelt1 * cdelt2
        self._etc["comp1"]["W"] = [W ,"Pixel solid angle (arcsec2/pix)", "PIX_ANG"]

        nbcell1 = self._etc["param"]["ccd"]["nbcell1"][0]
        naxis1 = int(nbcell1 / bin1)
        self._etc["comp1"]["naxis1"] = [naxis1, "Number of pixels along axis 1", "NAXIS1"]
        FoV1 = 2 * math.atan ( nbcell1 * photocell1 / Foclen / 2.) * 180. / math.pi
        self._etc["comp1"]["FoV1"] = [FoV1, "Field of view of the CCD image on axis1 (deg)", "FOV1"]

        nbcell2 = self._etc["param"]["ccd"]["nbcell2"][0]
        naxis2 = int(nbcell2 / bin2)
        self._etc["comp1"]["naxis2"] = [naxis2, "Number of pixels along axis 2", "NAXIS2"]
        FoV2 = 2 * math.atan ( nbcell2 * photocell2 / Foclen / 2.) * 180. / math.pi
        self._etc["comp1"]["FoV2"] = [FoV2, "Field of view of the CCD image on axis2 (deg)", "FOV2"]

        seeing = self._etc["param"]["local"]["seeing"][0]
        Fwhm_psf_seeing = seeing / 3600. * math.pi / 180 * Foclen
        self._etc["comp1"]["Fwhm_psf_seeing"] = [Fwhm_psf_seeing, "Fwhm of the seeing in the image plane (m)", "LSEEING"]

        Fwhm_psf_opt = self._etc["param"]["optic"]["Fwhm_psf_opt"][0]
        Fwhm_psf = math.sqrt ( Fwhm_psf_opt * Fwhm_psf_opt + Fwhm_psf_seeing * Fwhm_psf_seeing )
        self._etc["comp1"]["Fwhm_psf"] = [Fwhm_psf, "Fwhm of the PSF in the image plane (m)", "LPSF"]

        # --- Optics : computation of the gaussian fraction covered by the brightest pixel
        oversampling = 12 ; # must be even and >10 to ensure a good resolution
        if pixsize1 >= pixsize2:
            p = pixsize2
            P = pixsize1
        else:
            p = pixsize1
            P = pixsize2

        dp  = p/oversampling
        sigma = Fwhm_psf / (2*math.sqrt(2*math.log(2)))
        sigma2 = sigma*sigma
        a1d = 1 / sigma / math.sqrt(2*math.pi)
        a2d = a1d*a1d

        x1 = -p/2. ; x2 = x1 + p
        y1 = -P/2. ; y2 = y1 + P
        som = 0
        xs = []
        x = x1
        while x <= x2:
            xs.append(x)
            x += dp
        ys = []
        y = y1
        while y <= y2:
            ys.append(y)
            y += dp
        for x in xs:
            dx2 = x*x
            for y in ys:
                dy2 = y*y
                d2 = dx2 + dy2
                som += math.exp(-0.5*d2/sigma2)
        fpix1 = a2d*dp*dp*som
        self._etc["comp1"]["fpix1"] = [fpix1, "Flux fraction in the brightest pixel in the favorable case (max flux at the center of the pixel)", "FLUPIX00"]

        x1 = 0. ; x2 = x1 + p
        y1 = 0. ; y2 = y1 + P
        som = 0
        xs = []
        x = x1
        while x <= x2:
            xs.append(x)
            x += dp
        ys = []
        y = y1
        while y <= y2:
            ys.append(y)
            y += dp
        for x in xs:
            dx2 = x*x
            for y in ys:
                dy2 = y*y
                d2 = dx2 + dy2
                som += math.exp(-0.5*d2/sigma2)
        fpix3 = a2d*dp*dp*som
        self._etc["comp1"]["fpix3"] = [fpix3, "Flux fraction in the brightest pixel in the worst case (max flux at the corner of the pixel)", "FLUPIX11"]

        x1 = 0. ; x2 = x1 + p
        y1 = -P/2. ; y2 = y1 + P
        som = 0
        xs = []
        x = x1
        while x <= x2:
            xs.append(x)
            x += dp
        ys = []
        y = y1
        while y <= y2:
            ys.append(y)
            y += dp
        for x in xs:
            dx2 = x*x
            for y in ys:
                dy2 = y*y
                d2 = dx2 + dy2
                som += math.exp(-0.5*d2/sigma2)
        fpix2 = a2d*dp*dp*som
        self._etc["comp1"]["fpix2"] = [fpix2, "Flux fraction in the brightest pixel in the intermediate case", "FLUPIX01"]

        # --- Mean value of fraction
        fpix = (1.0*fpix1+2.0*fpix2+1.0*fpix3)/4.0
        self._etc["comp1"]["fpix"] = [fpix, "Flux fraction in the mean case", "FLUPIX"]

        # --- Object
        Fm0 = self._etc["param"]["filter"]["Fm0"][0]
        m = self._etc["input"]["object"]["m"][0]
        F_Jy = Fm0 * math.pow(10,-0.4 * m)
        self._etc["comp1"]["F_Jy"] = [F_Jy, "Total flux of the object outside atmosphere (Jy)", "FLUX_JY"]

        Dl = self._etc["param"]["filter"]["Dl"][0]
        l = self._etc["param"]["filter"]["l"][0]
        F_ph = F_Jy * 1.51e7 * Dl/l
        self._etc["comp1"]["F_ph"] = [F_ph, "Total flux of the object outside atmosphere (photons / sec /m2)", "FLUX_PH"]

        Tatm0 = self._etc["param"]["local"]["Tatm0"][0]
        Elev = self._etc["param"]["local"]["Elev"][0]
        Tatm = Tatm0 * math.sin(Elev*math.pi/180.)
        self._etc["comp1"]["Tatm"] = [Tatm, "Transmission of the atmosphere at elevation", "TRANSATM"]

        Topt = self._etc["param"]["optic"]["Topt"][0]
        t = self._etc["input"]["ccd"]["t"][0]
        Ftot_ph = F_ph * math.pi * D*D / 4. * Tatm * Topt * t
        self._etc["comp1"]["Ftot_ph"] = [Ftot_ph, "Total flux of the object after passed thru the optics (photons / object)", "FTOT_PH"]

        eta = self._etc["param"]["ccd"]["eta"][0]
        Ftot_el = Ftot_ph * eta
        self._etc["comp1"]["Ftot_el"] = [Ftot_el, "Total flux of the object after passed thru the optics (electrons / object)", "FTOT_EL"]

        Fpix_el = Ftot_el * fpix
        self._etc["comp1"]["Fpix_el"] = [Fpix_el, "Brightest pixel flux of the object after passed thru the optics (electrons / pixel)", "FPIX_EL"]

        # --- Sky brightness
        msky = self._etc["param"]["local"]["msky"][0]
        Sky_Jy = Fm0 * math.pow(10,-0.4 * msky)
        self._etc["comp1"]["Sky_Jy"] = [Sky_Jy, "Brightness of the sky (Jy/arsec2)", "BSKY_JY"]

        Sky_ph = Sky_Jy * 1.51e7 * Dl/l
        self._etc["comp1"]["Sky_ph"] = [Sky_ph, "Brightness of the sky (photons / sec /m2)", "BSKY_PH"]

        Skypix_ph = Sky_ph * math.pi * D*D / 4. * W * Topt * t
        self._etc["comp1"]["Skypix_ph"] = [Skypix_ph, "Brightness of the sky after passed thru the optics (photons / pixel)", "BPIX_JY"]

        Skypix_el = Skypix_ph * eta
        self._etc["comp1"]["Skypix_el"] = [Skypix_el, "Brightness of the sky after passed thru the optics (electrons / pixel)", "BPIX_EL"]

        # --- EMCCD
        Em = self._etc["param"]["ccd"]["Em"][0]
        fex = 1. + math.pow( 2./math.pi * math.atan( (Em-1)*3 ) ,3)
        self._etc["comp1"]["fex"] = [fex, "EMCCD excess noise factor (empirical formula derived from a figure of a paper)", "CCD_EM"]

    def t2snr_computations(self):
        """Compute the Signal to Noise Ratio (SNR) from a given exposure time

        Before using t2snr_computations, parameters must be set:

        ::

            etc = ExposureTimeCalculator()
            camera = "ProLine 16803"
            etc.camera(camera)
            optic = "Takahashi_180ED"
            etc.optics("Takahashi_180ED")
            etc.params("Fwhm_psf_opt", 15e-6)
            etc.params("band","C")
            etc.params("seeing", 2.5)
            t = 100
            etc.inputs("t", t)
            m = 16.5
            etc.inputs("m", m)
            snr = etc.t2snr_computations()

        """
        #---------------------------------------------------------
        ## @brief compute SNR from a given time
        #  @code What is the SNR for t=20s ?
        #  init
        #  inputs t 20
        #  disp
        #  t2snr_computations
        #  @endcode
        #  @return SNR ratio
        #
        self._preliminary_computations()
        C_th = self._etc["param"]["ccd"]["C_th"][0]
        bin1 = self._etc["param"]["ccd"]["bin1"][0]
        bin2 = self._etc["param"]["ccd"]["bin2"][0]
        t = self._etc["input"]["ccd"]["t"][0]
        Em = self._etc["param"]["ccd"]["Em"][0]
        S_th = C_th * bin1 * bin2 * t * Em
        self._etc["compsnr"]["S_th"] = [S_th, "Thermic signal (electrons/pixel)", "S_TH_E"]

        Skypix_el = self._etc["comp1"]["Skypix_el"][0]
        S_sk = Skypix_el * Em
        self._etc["compsnr"]["S_sk"] = [S_sk, "Sky signal (electrons/pixel)", "S_SKY_E"]

        Fpix_el = self._etc["comp1"]["Fpix_el"][0]
        S_ph = Fpix_el * Em
        self._etc["compsnr"]["S_ph"] = [S_ph, "Object signal (electrons/pixel)", "S_OBJ_E"]

        fex = self._etc["comp1"]["fex"][0]
        N_th = math.sqrt(C_th * bin1 * bin2 * t * fex) * Em
        self._etc["compsnr"]["N_th"] = [N_th, "Thermic noise (electrons/pixel)", "N_TH_E"]

        N_sk = math.sqrt(Skypix_el * fex) * Em
        self._etc["compsnr"]["N_sk"] = [N_sk, "Sky noise (electrons/pixel)", "N_SKY_E"]

        N_ph = math.sqrt(Fpix_el * fex) * Em
        self._etc["compsnr"]["N_ph"] = [N_ph, "Object noise (electrons/pixel) = shot noise", "N_OBJ_E"]

        N_ro = self._etc["param"]["ccd"]["N_ro"][0]
        N_tot = math.sqrt ( N_ro*N_ro + N_th*N_th + N_sk*N_sk + N_ph*N_ph )
        self._etc["compsnr"]["N_tot"] = [N_tot, "Total noise (electrons/pixels)", "N_TOT_E"]

        SNR_obj = S_ph / N_tot
        self._etc["compsnr"]["SNR_obj"] = [SNR_obj, "Object signal/noise at the brightest pixel", "SNRPIX00"]
        self._etc["compsnr"]["t"][0] = t # ?

        G = self._etc["param"]["ccd"]["G"][0]
        S_th_adu = S_th / G
        self._etc["compsnr"]["S_th_adu"] = [S_th_adu, "Thermic signal (ADU/pixel)", "S_TH_A"]

        S_sk_adu = S_sk / G
        self._etc["compsnr"]["S_sk_adu"] = [S_sk_adu, "Sky signal (ADU/pixel)", "S_SKY_A"]

        S_ph_adu = S_ph / G
        self._etc["compsnr"]["S_ph_adu"] = [S_ph_adu, "Object signal (ADU/pixel)", "S_OBJ_A"]

        N_th_adu = N_th / G
        self._etc["compsnr"]["N_th_adu"] = [N_th_adu, "Thermic noise (ADU/pixel)", "N_TH_A"]

        N_sk_adu = N_sk / G
        self._etc["compsnr"]["N_sk_adu"] = [N_sk_adu, "Sky noise (ADU/pixel)", "N_SKY_A"]

        N_ph_adu = N_ph / G
        self._etc["compsnr"]["N_ph_adu"] = [N_ph_adu, "Object noise (ADU/pixel) = shot noise", "N_OBJ_A"]

        N_tot_adu = N_tot / G
        self._etc["compsnr"]["N_tot_adu"] = [N_tot_adu, "Total noise (ADU/pixels)", "N_TOT_A"]

        return SNR_obj


    def snr2t_computations(self):
        """compute exposure time from a given Signal to Noise Ratio (SNR)

        Before using snr2t_computations, parameters must be set:

        ::

            etc = ExposureTimeCalculator()
            camera = "ProLine 16803"
            etc.camera(camera)
            optic = "Takahashi_180ED"
            etc.optics("Takahashi_180ED")
            etc.params("Fwhm_psf_opt", 15e-6)
            etc.params("band","C")
            etc.params("seeing", 2.5)
            m = 16.5
            etc.inputs("m", m)
            snr = 5
            etc.inputs("snr", snr)
            t = etc.snr2t_computations()

        """
        #---------------------------------------------------------
        ## @brief compute time from a given SNR
        #  @code What is the exposure time for SNR = 5 and magnitude = 18 ?
        #  init
        #  inputs snr 5
        #  inputs m 18
        #  disp
        #  snr2t_computations
        #  @endcode
        #  @return exposure time
        #
        self._preliminary_computations()
        snr = self._etc["input"]["constraint"]["snr"][0]
        N_ro = self._etc["param"]["ccd"]["N_ro"][0]
        C_th = self._etc["param"]["ccd"]["C_th"][0]
        bin1 = self._etc["param"]["ccd"]["bin1"][0]
        bin2 = self._etc["param"]["ccd"]["bin2"][0]
        Em = self._etc["param"]["ccd"]["Em"][0]
        fex = self._etc["comp1"]["fex"][0]
        Sky_ph = self._etc["comp1"]["Sky_ph"][0]
        W = self._etc["comp1"]["W"][0]
        D = self._etc["param"]["optic"]["D"][0]
        Topt = self._etc["param"]["optic"]["Topt"][0]
        eta = self._etc["param"]["ccd"]["eta"][0]
        F_ph = self._etc["comp1"]["F_ph"][0]
        Tatm = self._etc["comp1"]["Tatm"][0]
        fpix = self._etc["comp1"]["fpix"][0]

        C = snr*snr * N_ro*N_ro
        B = snr*snr * ( (C_th * bin1 * bin2 * fex * Em*Em) + (Sky_ph * math.pi * D*D / 4. * W * Topt * eta * fex * Em*Em) + (F_ph * math.pi * D*D / 4. * Tatm * Topt * eta * fpix * fex * Em*Em) )
        A = -math.pow( F_ph * math.pi * D*D / 4. * Tatm * Topt * eta *fpix * Em , 2)
        # We have A<0, B>0 and C >0. From the equation A*t^2 + B*t + C = 0, we can find the t value:
        D = B*B - 4*A*C # (always positive)
        t = (-B - math.sqrt(D) ) / (2.*A)
        self._etc["compsnr"]["t"] = [t, "Exposure time computer from a SNR value constrained"]
        return t


    def snr2m_computations(self) -> float:
        """Compute apparent magnitude from a given Signal to Noise Ratio (SNR)

        Before using snr2m_computations, parameters must be set:

        ::

            etc = ExposureTimeCalculator()
            camera = "ProLine 16803"
            etc.camera(camera)
            optic = "Takahashi_180ED"
            etc.optics("Takahashi_180ED")
            etc.params("Fwhm_psf_opt", 15e-6)
            etc.params("band","C")
            etc.params("seeing", 2.5)
            t = 100
            etc.inputs("t", t)
            snr = 5
            etc.inputs("snr", snr)
            m = etc.snr2m_computations()

        """
        self._preliminary_computations()
        snr = self._etc["input"]["constraint"]["snr"][0]
        N_ro = self._etc["param"]["ccd"]["N_ro"][0]
        C_th = self._etc["param"]["ccd"]["C_th"][0]
        bin1 = self._etc["param"]["ccd"]["bin1"][0]
        bin2 = self._etc["param"]["ccd"]["bin2"][0]
        Em = self._etc["param"]["ccd"]["Em"][0]
        fex = self._etc["comp1"]["fex"][0]
        Sky_ph = self._etc["comp1"]["Sky_ph"][0]
        W = self._etc["comp1"]["W"][0]
        D = self._etc["param"]["optic"]["D"][0]
        Topt = self._etc["param"]["optic"]["Topt"][0]
        eta = self._etc["param"]["ccd"]["eta"][0]
        F_ph = self._etc["comp1"]["F_ph"][0]
        Tatm = self._etc["comp1"]["Tatm"][0]
        fpix = self._etc["comp1"]["fpix"][0]
        Dl = self._etc["param"]["filter"]["Dl"][0]
        l = self._etc["param"]["filter"]["l"][0]
        Fm0 = self._etc["param"]["filter"]["Fm0"][0]

        t = self._etc["input"]["ccd"]["t"][0]

        C = snr*snr * ( N_ro*N_ro +  C_th * bin1 * bin2 * fex * Em*Em*t + Sky_ph * math.pi * D*D / 4. * W * Topt * eta * fex * Em*Em*t)
        B = snr*snr * ( math.pi * D*D / 4. * Tatm * Topt * eta * fpix * fex * Em*Em*t)
        A = -math.pow(math.pi * D*D / 4. * Tatm * Topt * eta * fpix * Em * t , 2)
        # We have A<0, B>0 and C >0. From the equation A*t^2 + B*t + C = 0, we can find the t value:
        D = B*B - 4*A*C # (always positive)
        F_ph = ( -B - math.sqrt(D) ) / (2.*A)
        F_Jy = F_ph / (1.51e7 * Dl/l)
        m = -2.5 *math.log10( F_Jy / Fm0)
        self._etc["compsnr"]["m"] = [m, "Apparent magnitude computed from a SNR and exposure value constrained"]
        return m


    def init (self, band="V", moon_age=0):
        #---------------------------------------------------------
        #  brief initializations
        #  param band (default) V
        #  param moon_age (devault) 0
        #
        self._etc = {}
        self._etc["compsnr"] = {}
        self._etc["compsnr"]["t"] = [1, "Exposure time computer from a SNR value constrained", "EXPTIME"]
        self._etc["compsnr"]["m"] = [10, "Apparent magnitude computed from a SNR and exposure value constrained", "MAG_APP"]
        self._etc["compsnr"]["SNR_obj"] = [5, "SNR computed from a exposure and a apparent magnitude value constrained", "SNR_OBJ"]

        self._params_defaults(band, moon_age)
        self._inputs_defaults()


    def _superkeys(self):
        superkeys = {}
        k1s = self._etc.keys()
        for k1 in k1s:
            k2s = self._etc[k1].keys()
            ik2 = 0
            for k2 in k2s:
                ik2 += 1
                if type(self._etc[k1][k2]) is list:
                    val = self._etc[k1][k2][0]
                    com = self._etc[k1][k2][1]
                    superkeys[k2] = [val, com, k1]
                else:
                    k3s = self._etc[k1][k2].keys()
                    ik3 = 0
                    for k3 in k3s:
                        ik3 += 1
                        val = self._etc[k1][k2][k3][0]
                        com = self._etc[k1][k2][k3][1]
                        superkeys[k3] = [val, com, (k1, k2)]
        return superkeys

    def __str__(self):
        msg = ""
        k1s = self._etc.keys()
        for k1 in k1s:
            k2s = self._etc[k1].keys()
            ik2 = 0
            for k2 in k2s:
                ik2 += 1
                if type(self._etc[k1][k2]) is list:
                    tk2 = 0
                    if ik2 == 1:
                        msg += f"=== {k1} ===\n"
                    val = self._etc[k1][k2][0]
                    com = self._etc[k1][k2][1]
                    msg += f"{k2} = {val} // {com}\n"
                else:
                    tk2 = 1
                    k3s = self._etc[k1][k2].keys()
                    ik3 = 0
                    for k3 in k3s:
                        ik3 += 1
                        if ik3 == 1:
                            msg += f"=== {k1} {k2} ===\n"
                        val = self._etc[k1][k2][k3][0]
                        com = self._etc[k1][k2][k3][1]
                        msg += f"{k3} = {val} // {com}\n"
                    msg += "\n"
            if tk2 == 0:
                msg += "\n"
        return msg

    def print(self, forkey=""):
        superkeys = self._superkeys()
        for k, v in superkeys.items():
            val, com, prefix = v
            if forkey == "" or forkey in prefix:
                print(f"{prefix}/{k} = {val} // {com}")

    def simu_star_params(self, unit: str="adu") -> dict:
        """Return parameters to simulate a Gaussian image of a star

        Args:
            unit: Choice for output (unit), "adu" or "electron".

        Returns:
            Dict of Gaussian parameters. Each key of the dictionary is a list of value, comment. The keys are:

                * 'max_sig': Signal in the central pixel (unit)
                * 'fwhmx': Full Width at Half Maximum along x axis (pixels)
                * 'fwhmy': Full Width at Half Maximum along y axis (pixels)
                * 'tot_sig': Integral signal of the star (unit)
                * 'sky_sig': Sky signal (unit)
                * 'sky_std': Sky standard deviation (unit)

        """
        Skypix_el = self._etc["comp1"]["Skypix_el"][0]
        G = self._etc["param"]["ccd"]["G"][0] # e/adu
        Ftot_el = self._etc["comp1"]["Ftot_el"][0]
        fpix = self._etc["comp1"]["fpix"][0]
        Fwhm_psf = self._etc["comp1"]["Fwhm_psf"][0] # m
        Foclen = self._etc["comp1"]["Foclen"][0] # m
        cdelt1 = self._etc["comp1"]["cdelt1"][0] / 3600.0  # deg/pix
        cdelt2 = self._etc["comp1"]["cdelt2"][0] / 3600.0  # deg/pix
        # --- electrons
        sky_sig = Skypix_el
        sky_std = math.sqrt(Skypix_el)
        tot_sig = Ftot_el # electrons in the integral
        max_sig = Ftot_el*fpix # electrons in the central pixel
        if unit == "adu":
            # --- adu
            sky_sig /= G
            sky_std /= G
            max_sig /= G # adu in the central pixel
            tot_sig /= G # adu in the integral
        # --- pixels
        Fwhm_psf = math.degrees(Fwhm_psf/Foclen) # deg
        fwhmx = Fwhm_psf/cdelt1
        fwhmy = Fwhm_psf/cdelt2
        res = {}
        res["max_sig"] = [max_sig, f"Signal in the central pixel ({unit})"]
        res["fwhmx"] = [fwhmx, "Full Width at Half Maximum along x axis (pixels)"]
        res["fwhmy"] = [fwhmy, "Full Width at Half Maximum along y axis (pixels)"]
        res["tot_sig"] = [tot_sig, f"Integral signal of the star ({unit})"]
        res["sky_sig"] = [sky_sig, f"Sky signal ({unit})"]
        res["sky_std"] = [sky_std, f"Sky standard deviation ({unit})"]
        return res


# #####################################################################
# #####################################################################
# #####################################################################
# Main
# #####################################################################
# #####################################################################
# #####################################################################

if __name__ == "__main__":

    default = 4
    example = input(f"Select the example (0 to 4) ({default}) ")
    try:
        example = int(example)
    except:
        example = default

    print("Example       = {}".format(example))

    if example == 1:
        """
        List all cameras and optics
        """
        etc = ExposureTimeCalculator()
        print("List of cameras ===")
        print(etc.camera())
        print("List of optics ===")
        print(etc.optics())

    if example == 2:
        """
        TRE: Compute SNR and limiting magnitude
        """
        etc = ExposureTimeCalculator()
        camera = "ProLine 16803"
        etc.camera(camera)
        optic = "Takahashi_180ED"
        etc.optics("Takahashi_180ED")
        # ---
        etc.params("Fwhm_psf_opt", 15e-6)
        band = etc.params("band","C")[0]
        etc.params("seeing", 2.5)
        # --- 1st calculation
        t = 90
        etc.inputs("t", t)
        m= 17
        etc.inputs("m", m)
        snr = etc.t2snr_computations()
        print(f"Setup: Camera={camera} Optic={optic}")
        print(f"For an exposure of {t}s a star of mag({band})={m} has a SNR = {snr:.1f}")
        # --- 2nd calculation
        snr = 5
        etc.inputs("snr", snr)
        m = etc.snr2m_computations()
        print(f"For an exposure of {t}s a star of SNR={snr} has a mag({band}) = {m:.1f}")
        print("=== Params for star simulation ===")
        print(etc.simu_star_params())

    if example == 3:
        """
        T1M: Compute SNR and limiting magnitude
        """
        etc = ExposureTimeCalculator()
        optic = "T1M"
        etc.params("D", 1.05)
        etc.params("FonD", 11.0) # 11.0 17.5
        etc.params("Topt", 0.6)
        etc.params("Fwhm_psf_opt", 15e-6)
        etc.params("diamFull", 30e-3)
        etc.params("optPrice", 1e5)
        camera = "Andor DW936 BV" # Andor DW936 BV Lytid SIRIS
        etc.camera(camera)
        # ---
        band = etc.params("band","J")[0]
        etc.params("seeing", 1.2)
        etc.params("moon_age", 14)
        # --- 1st calculation
        etc.params("bin1", 2)
        etc.params("bin2", 2)
        t = 100
        etc.inputs("t", t)
        m= 17
        etc.inputs("m", m)
        snr = etc.t2snr_computations()
        print(f"Setup: Camera={camera} Optic={optic}")
        print(f"For an exposure of {t}s a star of mag({band})={m} has a SNR = {snr:.1f}")
        # --- 2nd calculation
        snr = 5
        etc.inputs("snr", snr)
        m = etc.snr2m_computations()
        print(f"For an exposure of {t}s a star of SNR={snr} has a mag({band}) = {m:.1f}")
        print("=== Params for star simulation ===")
        print(etc.simu_star_params())

    if example == 4:
        """
        Zadko: Compute SNR and limiting magnitude
        """
        etc = ExposureTimeCalculator()
        optic = "Zadko"
        etc.params("D", 1.0)
        etc.params("FonD", 4.0) # 11.0 17.5
        etc.params("Topt", 0.6)
        etc.params("Fwhm_psf_opt", 60e-6)
        etc.params("diamFull", 0.1)
        etc.params("optPrice", 1e5)
        camera = "Kepler 4040"
        etc.camera(camera)
        # ---
        band = etc.params("band","R")[0]
        etc.params("seeing", 1.5)
        etc.params("moon_age", 0)
        # --- 1st calculation
        etc.params("bin1", 1)
        etc.params("bin2", 1)
        t = 100
        etc.inputs("t", t)
        m= 17
        etc.inputs("m", m)
        snr = etc.t2snr_computations()
        print(f"Setup: Camera={camera} Optic={optic}")
        print(f"For an exposure of {t}s a star of mag({band})={m} has a SNR = {snr:.1f}")
        # --- 2nd calculation
        snr = 5
        etc.inputs("snr", snr)
        m = etc.snr2m_computations()
        print(f"For an exposure of {t}s a star of SNR={snr} has a mag({band}) = {m:.1f}")
        print("=== Params for star simulation ===")
        print(etc.simu_star_params())