#! /usr/bin/python3

# coding : utf-8

 

# Programme calculant la magnitude limite stellaire en CCD.

 

# Source : "La détection des objets faibles du Système Solaire",
# Christian Buil, CCD et Télescope, n° 7, Automne 1996.
# « Voici un petit formulaire que j'utilise couramment et qui
# donne une précision correcte (±0,5 magnitude typiquement). »

 

#----------------------------------------------------------------------#
# Importations

 

from sys  import argv
from math import pi
from math import sqrt
from math import atan
from math import log10

 

#----------------------------------------------------------------------#
# Paramètres de la ligne de commande :
# argv[1] : diamètre du télescope en millimètres
# argv[2] : focale du télescope en millimètres
# argv[3] : taille d'un pixel en microns
# argv[4] : nombre d'images unitaires
# argv[5] : temps de pose unitaire en secondes

 

if len(argv) <= 5 :
    print('Erreur : argument(s) manquant(s)')
    print('Syntaxe : ')
    print('> arguments: Diametre(mm) Focale(mm) pixel_size(µm) Nb_poses Tps_pose_unitaire(s)')
    exit(1)
argv1 = float(argv[1])
argv2 = float(argv[2])
argv3 = float(argv[3])
argv4 = int  (argv[4])
argv5 = float(argv[5])
print()

 

#----------------------------------------------------------------------#
# Constantes

 

# E0 = « éclairement d'une étoile de magnitude 0 et de type solaire en
# dehors de l'atmosphère en photons/cm²/s/micron. Autour de 0,65 microns
# nous avons E0 = 0,8.10⁷ ph/cm²/s/microns. »
E0 = 0.8*10**7

 

# Tau = « transmission totale du système de détection. Ce paramètre inclut
# la transmission de l'atmosphère, le coefficient de réflexion des miroirs,
# les pertes par réflexions sur les dioptres rencontrés et les pertes dues
# à l'obstruction centrale. Dans un bon site de montagne et avec un
# télescope de type Newton d'obstruction centrale 0,33 équipée d'un
# correcteur de champ et pour une visée à une distance zénithale de 30°
# on trouve Tau = 0,53. La même instrumentation utilisée en plaine donnera
# Tau = 0,43.
Tau = 0.4

 

# QE = « rendement quantique du CCD. Nous adoptons ici la valeur de 0,33
# qui est une moyenne avec un CCD du type Kodak-0400 utilisé sans filtre. »
# (Attention : c'est donc la moyenne sur la bande spectrale et non
# le maximum.)
# Default value : 0.33
QE = 0.80 #Adapté pour un capteur CMOS

 

# Dl = « largeur de la bande spectrale en microns. Adopter 0,4 microns
# si on désire se mettre dans la situation d'une caméra CCD utilisée
# sans filtre et utiliser 0,1 micron pour simuler l'utilisation d'un
# filtre R (dans ce cas faire QE = 0,37).
Dl = 0.4

 

# Msky = « magnitude du fond du ciel par seconde d'arc carrée. La valeur
# typique dans un bon site amateur est msky = 20,4 dans le rouge. »
# Bortle 8 = 18.00
# Bortle 7 = 18.30
# Bortle 6 = 18.80
# Bortle 5 = 19.80
# Bortle 4 = 20.90
# Exemple calculés avec l'outil de sharpcap -> http://www.tools.sharpcap.co.uk/
Msky = 18.80

 

# FWHM = « largeur à mi-hauteur des étoiles en secondes d'arc (ce
# paramètre tient compte aussi des défauts de suivi et des aberrations
# optiques). »
FWHM = 3.5

 

# Dark = « taux de production de charges thermiques par seconde.
# Un CCD MPP comme le Kodak-0400 refroidi vers -10°C produit environ
# 0,1 électron/seconde/pixel. »
Dark = 0.0622	# Valeur trouvée concernant l'imx290 à -10°C

 

# Sigma = « bruit de lecture de la caméra en électrons. Nous adopterons
# ici Sigma = 15 électrons qui est une bonne valeur typique pour une
# caméra de qualité aujourd'hui. »
SigmaCam = 1	# Modifié pour collé à un capteur CMOS

 

# SB = « rapport signal sur bruit à partir duquel on estime que la
# détection de l'objet est effective. En théorie S/B = .3, mais la
# pratique montre qu'il est opportun d'adopter une marge. Un S/B=5
# est généralement un bon choix.
SB = 4.0

 

#----------------------------------------------------------------------#
# Calculs 

 

# diam = diamètre du miroir primaire en cm
diam = 0.1 * argv1

 

# pix = « taille du pixel en secondes d'arc. »
pix = 7200.0 * atan(argv3 / (2000.0 * argv2)) * 180.0 / pi
print('    Échantillonnage    = {:4.2f}"'.format(pix))

 

# tpose = « temps de pose effectif en secondes »
tpose = argv4 * argv5
hpose = int(tpose / 3600.0)
treste = tpose - 3600.0 * hpose
mpose = int(treste / 60.0)
spose = treste - 60.0 * mpose
print('    Tps de pose total  = {:01d}h{:02d}m{:02.0f}s'.format(hpose, mpose, spose))

 

# σ² bruit du fond du ciel
vciel = 0.7 * E0 * Tau * diam**2 * Dl * pix**2 * 10**(-0.4*Msky) * tpose
sigmaciel = sqrt(vciel)
print('    σ bruit fond ciel  = {:6.2f}'.format(sigmaciel))

 

# σ² bruit thermique 
vdark = Dark * tpose
sigmadark = sqrt(vdark)
print('    σ bruit thermique  = {:6.2f}'.format(sigmadark))

 

# σ² bruit de lecture
vlect = argv4 * SigmaCam**2
sigmalect = sqrt(vlect)
print('    σ bruit de lecture = {:6.2f}'.format(sigmalect))

 

# σ² total
vtot =  vciel + vdark + vlect
sigmatot = sqrt(vtot)
print('    σ bruit total img  = {:6.2f}'.format(sigmatot))

 

# Fraction de signal de l'étoile tombant dans le pixel central de
# la tache image
fract = (4.0 * pix**2) / (pi * FWHM**2)
print('    Fraction ét./pixel = {:6.3f}'.format(fract), end="")
if fract < 0.5 :
    print()
else :
    print(' --> sous-échantillonnage, le calcul sera faussé')

 

# Magnitude limite
x1 = 1.0 + sqrt(1.0 + (4.0 * vtot) / SB**2)
x2 = SB**2 / (1.4 * E0 * QE * Tau * fract * Dl * diam**2 * tpose)
magn = -2.5 * log10(x1 * x2)
print()
print('--> Magnitude limite   = {:4.1f}'.format(magn))
print()